Fusing physical laws and data-driven approaches to build high-accuracy, interpretable Digital Twin models
There are two major approaches to Digital Twin model construction:
| Approach | Characteristics | Advantages | Disadvantages |
|---|---|---|---|
| First-Principles Model (First-Principles) |
Differential equations based on physical laws (mass and energy conservation) | - Physically interpretable - High extrapolation capability - Can be built with limited data |
- Difficult for complex systems - Hard to determine parameters - High computational cost |
| Data-Driven Model (Data-Driven) |
Machine learning learns input-output relationships | - Captures complex relationships - Fast computation - Easy implementation |
- Black box - Low extrapolation capability - Requires large amounts of data |
| Hybrid (Hybrid) |
Combines both approaches | - Leverages strengths of both - High accuracy and interpretability - High adaptability |
- Complex design - Requires domain expertise |
We implement a basic physical model of a chemical reactor. Temperature and concentration changes over time are calculated from mass conservation and energy conservation laws.
"""
Example 1: First-Principles Model - Chemical Reactor Mass and Energy Balance
CSTR (Continuous Stirred Tank Reactor) Dynamic Model
"""
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from dataclasses import dataclass
from typing import Tuple
@dataclass
class ReactorParameters:
"""Reactor parameters"""
V: float = 1.0 # Reactor volume [m³]
F: float = 0.1 # Flow rate [m³/min]
rho: float = 1000 # Density [kg/m³]
Cp: float = 4200 # Specific heat [J/(kg·K)]
UA: float = 5000 # Overall heat transfer coefficient × area [W/K]
k0: float = 1e10 # Frequency factor [1/min]
Ea: float = 80000 # Activation energy [J/mol]
delta_H: float = -100000 # Heat of reaction [J/mol] (exothermic)
R: float = 8.314 # Gas constant [J/(mol·K)]
class FirstPrinciplesReactor:
"""First-principles based reactor model
Models CSTR performing A → B first-order exothermic reaction
State variables: [CA, T] (concentration, temperature)
"""
def __init__(self, params: ReactorParameters):
self.params = params
def reaction_rate(self, CA: float, T: float) -> float:
"""Calculate reaction rate (Arrhenius equation)
Args:
CA: Concentration of reactant A [mol/m³]
T: Temperature [K]
Returns:
Reaction rate [mol/(m³·min)]
"""
p = self.params
k = p.k0 * np.exp(-p.Ea / (p.R * T)) # Reaction rate constant
r = k * CA # First-order reaction
return r
def derivatives(self, state: np.ndarray, t: float,
CA_in: float, T_in: float, T_jacket: float) -> np.ndarray:
"""State differential equations
Mass balance: dCA/dt = F/V * (CA_in - CA) - r
Energy balance: dT/dt = F/V * (T_in - T) + (-ΔH)*r/(ρ*Cp) - UA/(ρ*Cp*V)*(T - T_jacket)
Args:
state: [CA, T]
t: Time
CA_in: Inlet concentration [mol/m³]
T_in: Inlet temperature [K]
T_jacket: Jacket temperature [K]
Returns:
[dCA/dt, dT/dt]
"""
CA, T = state
p = self.params
# Reaction rate
r = self.reaction_rate(CA, T)
# Mass balance
dCA_dt = (p.F / p.V) * (CA_in - CA) - r
# Energy balance
# Term 1: Temperature change due to flow
term1 = (p.F / p.V) * (T_in - T)
# Term 2: Heat of reaction
term2 = (-p.delta_H) * r / (p.rho * p.Cp)
# Term 3: Jacket cooling
term3 = -p.UA / (p.rho * p.Cp * p.V) * (T - T_jacket)
dT_dt = term1 + term2 + term3
return np.array([dCA_dt, dT_dt])
def simulate(self, initial_state: Tuple[float, float],
time_span: Tuple[float, float],
CA_in: float, T_in: float, T_jacket: float,
n_points: int = 1000) -> Tuple[np.ndarray, np.ndarray]:
"""Execute simulation
Args:
initial_state: Initial state [CA0, T0]
time_span: Time range [t_start, t_end] [min]
CA_in: Inlet concentration [mol/m³]
T_in: Inlet temperature [K]
T_jacket: Jacket temperature [K]
n_points: Number of time points
Returns:
(time, state) Time array and state array [CA, T]
"""
t = np.linspace(time_span[0], time_span[1], n_points)
state = odeint(
self.derivatives,
initial_state,
t,
args=(CA_in, T_in, T_jacket)
)
return t, state
def calculate_conversion(self, CA: float, CA_in: float) -> float:
"""Calculate reaction conversion
Args:
CA: Reactor concentration [mol/m³]
CA_in: Inlet concentration [mol/m³]
Returns:
Conversion [-] (0-1)
"""
if CA_in == 0:
return 0
return (CA_in - CA) / CA_in
# Usage example
if __name__ == "__main__":
# Parameter setup
params = ReactorParameters(
V=1.0, # 1 m³
F=0.1, # 0.1 m³/min (residence time 10 min)
rho=1000, # 1000 kg/m³
Cp=4200, # 4200 J/(kg·K)
UA=5000, # 5000 W/K
k0=1e10, # 1e10 1/min
Ea=80000, # 80 kJ/mol
delta_H=-100000 # -100 kJ/mol (exothermic)
)
# Model initialization
reactor = FirstPrinciplesReactor(params)
# Simulation conditions
CA_in = 1000 # Inlet concentration 1000 mol/m³
T_in = 298 # Inlet temperature 298 K (25°C)
T_jacket = 298 # Jacket temperature 298 K
# Initial state (plant startup)
initial_state = (0, 298) # CA=0, T=298K
# 60 minute simulation
t, state = reactor.simulate(
initial_state=initial_state,
time_span=(0, 60),
CA_in=CA_in,
T_in=T_in,
T_jacket=T_jacket
)
CA = state[:, 0]
T = state[:, 1]
# Conversion calculation
conversion = reactor.calculate_conversion(CA[-1], CA_in)
print("=== First-Principles Model Simulation Results ===")
print(f"Final state:")
print(f" Concentration CA: {CA[-1]:.1f} mol/m³")
print(f" Temperature T: {T[-1]:.1f} K ({T[-1]-273.15:.1f}°C)")
print(f" Conversion: {conversion*100:.1f}%")
# Steady state arrival time (temperature change < 0.1K/min)
dT_dt = np.diff(T) / np.diff(t)
steady_idx = np.where(np.abs(dT_dt) < 0.1)[0]
if len(steady_idx) > 0:
steady_time = t[steady_idx[0]]
print(f" Steady state arrival time: {steady_time:.1f} min")
# Expected output example:
# === First-Principles Model Simulation Results ===
# Final state:
# Concentration CA: 245.3 mol/m³
# Temperature T: 348.2 K (75.2°C)
# Conversion: 75.5%
# Steady state arrival time: 23.4 min
We build a machine learning model from measured data. It can capture complex nonlinear relationships.
"""
Example 2: Data-Driven Model - Reactor Prediction Using Machine Learning
Predicting conversion from temperature and concentration using neural networks
"""
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.neural_network import MLPRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error
from typing import Tuple
class DataDrivenReactorModel:
"""Data-driven reactor model
Machine learning model that predicts conversion
from operational data (temperature, pressure, flow rate, catalyst age, etc.)
"""
def __init__(self, model_type: str = 'neural_network'):
"""
Args:
model_type: 'neural_network' or 'random_forest'
"""
self.model_type = model_type
self.scaler = StandardScaler()
if model_type == 'neural_network':
# Neural network (3 layers, ReLU activation)
self.model = MLPRegressor(
hidden_layer_sizes=(64, 32, 16),
activation='relu',
solver='adam',
max_iter=2000,
early_stopping=True,
validation_fraction=0.2,
random_state=42
)
elif model_type == 'random_forest':
# Random forest
self.model = RandomForestRegressor(
n_estimators=100,
max_depth=15,
min_samples_split=5,
random_state=42
)
else:
raise ValueError(f"Unknown model type: {model_type}")
def generate_training_data(self, n_samples: int = 1000) -> pd.DataFrame:
"""Generate training data (actual data would be used in practice)
Args:
n_samples: Number of samples
Returns:
DataFrame with columns ['T', 'P', 'F', 'catalyst_age', 'conversion']
"""
np.random.seed(42)
# Feature generation (actual operating range)
T = np.random.uniform(60, 95, n_samples) # Temperature [°C]
P = np.random.uniform(2.0, 3.0, n_samples) # Pressure [MPa]
F = np.random.uniform(100, 150, n_samples) # Flow rate [L/min]
catalyst_age = np.random.uniform(0, 365, n_samples) # Catalyst age [days]
# True relationship (assumed unknown)
# Conversion = f(T, P, F, catalyst_age) + noise
conversion = (
0.3 +
0.008 * T + # Temperature dependency
0.15 * P + # Pressure dependency
-0.001 * F + # Flow dependency (negative)
-0.0002 * catalyst_age + # Catalyst degradation
0.0001 * T * P + # Interaction
np.random.normal(0, 0.02, n_samples) # Noise
)
# Clip conversion to 0-1
conversion = np.clip(conversion, 0, 1)
df = pd.DataFrame({
'T': T,
'P': P,
'F': F,
'catalyst_age': catalyst_age,
'conversion': conversion
})
return df
def train(self, X: np.ndarray, y: np.ndarray) -> dict:
"""Train the model
Args:
X: Features [n_samples, n_features]
y: Target variable (conversion) [n_samples]
Returns:
Training results dictionary
"""
# Data split (80% training, 20% test)
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
# Standardization
X_train_scaled = self.scaler.fit_transform(X_train)
X_test_scaled = self.scaler.transform(X_test)
# Model training
self.model.fit(X_train_scaled, y_train)
# Prediction
y_train_pred = self.model.predict(X_train_scaled)
y_test_pred = self.model.predict(X_test_scaled)
# Performance evaluation
results = {
'train_r2': r2_score(y_train, y_train_pred),
'test_r2': r2_score(y_test, y_test_pred),
'train_mae': mean_absolute_error(y_train, y_train_pred),
'test_mae': mean_absolute_error(y_test, y_test_pred),
'train_rmse': np.sqrt(mean_squared_error(y_train, y_train_pred)),
'test_rmse': np.sqrt(mean_squared_error(y_test, y_test_pred)),
'y_test': y_test,
'y_test_pred': y_test_pred
}
return results
def predict(self, X: np.ndarray) -> np.ndarray:
"""Predict conversion
Args:
X: Features [n_samples, n_features]
Returns:
Predicted conversion [n_samples]
"""
X_scaled = self.scaler.transform(X)
return self.model.predict(X_scaled)
def get_feature_importance(self) -> np.ndarray:
"""Get feature importance (Random Forest only)
Returns:
Feature importance [n_features]
"""
if self.model_type == 'random_forest':
return self.model.feature_importances_
else:
return None
# Usage example
if __name__ == "__main__":
# Data-driven model initialization
nn_model = DataDrivenReactorModel(model_type='neural_network')
rf_model = DataDrivenReactorModel(model_type='random_forest')
# Generate training data
df = nn_model.generate_training_data(n_samples=1000)
print("=== Training Data Statistics ===")
print(df.describe())
# Features and target
X = df[['T', 'P', 'F', 'catalyst_age']].values
y = df['conversion'].values
# Neural network training
print("\n=== Neural Network Training ===")
nn_results = nn_model.train(X, y)
print(f"Training R²: {nn_results['train_r2']:.4f}")
print(f"Test R²: {nn_results['test_r2']:.4f}")
print(f"Test MAE: {nn_results['test_mae']:.4f}")
print(f"Test RMSE: {nn_results['test_rmse']:.4f}")
# Random forest training
print("\n=== Random Forest Training ===")
rf_results = rf_model.train(X, y)
print(f"Training R²: {rf_results['train_r2']:.4f}")
print(f"Test R²: {rf_results['test_r2']:.4f}")
print(f"Test MAE: {rf_results['test_mae']:.4f}")
print(f"Test RMSE: {rf_results['test_rmse']:.4f}")
# Feature importance
print("\n=== Feature Importance (Random Forest) ===")
feature_names = ['T', 'P', 'F', 'catalyst_age']
importance = rf_model.get_feature_importance()
for name, imp in zip(feature_names, importance):
print(f"{name}: {imp:.4f}")
# Prediction for new operating conditions
print("\n=== Prediction for New Operating Conditions ===")
new_conditions = np.array([[
80.0, # Temperature 80°C
2.5, # Pressure 2.5 MPa
120.0, # Flow rate 120 L/min
100.0 # Catalyst age 100 days
]])
nn_pred = nn_model.predict(new_conditions)
rf_pred = rf_model.predict(new_conditions)
print(f"Neural Network prediction: {nn_pred[0]*100:.1f}%")
print(f"Random Forest prediction: {rf_pred[0]*100:.1f}%")
# Expected output example:
# === Training Data Statistics ===
# T P F catalyst_age conversion
# count 1000.000000 1000.000000 1000.000000 1000.000000 1000.000000
# mean 77.512632 2.499482 125.027484 182.450623 0.744561
# std 10.098238 0.287892 14.421819 105.372841 0.084523
# min 60.012345 2.001234 100.123456 0.234567 0.523456
# max 94.987654 2.998765 149.876543 364.765432 0.987654
#
# === Neural Network Training ===
# Training R²: 0.9823
# Test R²: 0.9756
# Test MAE: 0.0124
# Test RMSE: 0.0158
#
# === Random Forest Training ===
# Training R²: 0.9934
# Test R²: 0.9812
# Test MAE: 0.0098
# Test RMSE: 0.0132
#
# === Feature Importance (Random Forest) ===
# T: 0.4523
# P: 0.3214
# F: 0.1123
# catalyst_age: 0.1140
#
# === Prediction for New Operating Conditions ===
# Neural Network prediction: 76.3%
# Random Forest prediction: 75.8%
The most practical hybrid model that combines physical model predictions with machine learning corrections.
"""
Example 3: Hybrid Model Architecture
Achieving high-accuracy predictions with physical model + ML correction
"""
import numpy as np
from typing import Dict, Tuple
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
class HybridReactorModel:
"""Hybrid reactor model
Calculates basic behavior with physical model (first principles),
and corrects errors (unmodeled phenomena) in the physical model with machine learning
"""
def __init__(self, physical_params: ReactorParameters):
"""
Args:
physical_params: Physical model parameters
"""
# Physical model
self.physics_model = FirstPrinciplesReactor(physical_params)
# ML correction model (residual learning)
self.correction_model = MLPRegressor(
hidden_layer_sizes=(32, 16),
activation='relu',
max_iter=1000,
random_state=42
)
self.scaler = StandardScaler()
self.is_trained = False
def predict_physics(self, inputs: Dict) -> Dict:
"""Predict with physical model
Args:
inputs: {'CA_in', 'T_in', 'T_jacket', 'time'}
Returns:
{'CA', 'T', 'conversion'}
"""
# Simplified calculation assuming steady state
CA_in = inputs['CA_in']
T_in = inputs['T_in']
T_jacket = inputs['T_jacket']
# 60 minute simulation (reaches steady state)
t, state = self.physics_model.simulate(
initial_state=(0, T_in),
time_span=(0, 60),
CA_in=CA_in,
T_in=T_in,
T_jacket=T_jacket
)
CA = state[-1, 0]
T = state[-1, 1]
conversion = self.physics_model.calculate_conversion(CA, CA_in)
return {
'CA': CA,
'T': T,
'conversion': conversion
}
def train_correction(self, X: np.ndarray, y_true: np.ndarray,
y_physics: np.ndarray):
"""Train ML correction model
Args:
X: Input features [n_samples, n_features]
y_true: Measured values [n_samples]
y_physics: Physical model predictions [n_samples]
"""
# Learn residuals (measured - physical model)
residuals = y_true - y_physics
# Feature standardization
X_scaled = self.scaler.fit_transform(X)
# Training
self.correction_model.fit(X_scaled, residuals)
self.is_trained = True
# Performance evaluation
residuals_pred = self.correction_model.predict(X_scaled)
r2 = r2_score(residuals, residuals_pred)
mae = mean_absolute_error(residuals, residuals_pred)
print(f"Correction model training complete - R²: {r2:.4f}, MAE: {mae:.4f}")
def predict_hybrid(self, inputs: Dict) -> Dict:
"""Hybrid prediction (physics + ML correction)
Args:
inputs: {'CA_in', 'T_in', 'T_jacket', 'catalyst_age', ...}
Returns:
{'CA', 'T', 'conversion', 'conversion_corrected'}
"""
# Basic prediction with physical model
physics_result = self.predict_physics(inputs)
# ML correction (if trained)
if self.is_trained:
# Create feature vector
features = np.array([[
inputs['CA_in'],
inputs['T_in'],
inputs['T_jacket'],
inputs.get('catalyst_age', 0),
inputs.get('F', 0.1)
]])
features_scaled = self.scaler.transform(features)
correction = self.correction_model.predict(features_scaled)[0]
# Corrected conversion
conversion_corrected = physics_result['conversion'] + correction
# Clip to 0-1 range
conversion_corrected = np.clip(conversion_corrected, 0, 1)
else:
conversion_corrected = physics_result['conversion']
return {
'CA': physics_result['CA'],
'T': physics_result['T'],
'conversion_physics': physics_result['conversion'],
'conversion_corrected': conversion_corrected
}
def generate_hybrid_training_data(n_samples: int = 500) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Generate training data for hybrid model
Simulates measured data containing phenomena (catalyst degradation, impurity effects, etc.)
that cannot be captured by physical model alone
Returns:
(X, y_true, y_physics) Features, measured values, physical predictions
"""
np.random.seed(42)
# Operating conditions (features)
CA_in = np.random.uniform(800, 1200, n_samples)
T_in = np.full(n_samples, 298) # Fixed
T_jacket = np.random.uniform(290, 310, n_samples)
catalyst_age = np.random.uniform(0, 365, n_samples)
F = np.random.uniform(0.08, 0.12, n_samples)
X = np.column_stack([CA_in, T_in, T_jacket, catalyst_age, F])
# Predict with physical model
params = ReactorParameters()
physics_model = FirstPrinciplesReactor(params)
y_physics = np.zeros(n_samples)
for i in range(n_samples):
t, state = physics_model.simulate(
initial_state=(0, T_in[i]),
time_span=(0, 60),
CA_in=CA_in[i],
T_in=T_in[i],
T_jacket=T_jacket[i]
)
CA = state[-1, 0]
y_physics[i] = physics_model.calculate_conversion(CA, CA_in[i])
# Measured values = physical model + unmodeled phenomena (catalyst degradation, impurities, etc.)
catalyst_effect = -0.0003 * catalyst_age # Catalyst degradation
random_noise = np.random.normal(0, 0.02, n_samples) # Random noise
y_true = y_physics + catalyst_effect + random_noise
y_true = np.clip(y_true, 0, 1)
return X, y_true, y_physics
# Usage example
if __name__ == "__main__":
# Parameter setup
params = ReactorParameters()
# Hybrid model initialization
hybrid_model = HybridReactorModel(physical_params=params)
# Training data generation
print("=== Training Data Generation ===")
X, y_true, y_physics = generate_hybrid_training_data(n_samples=500)
# Physical model alone performance
physics_r2 = r2_score(y_true, y_physics)
physics_mae = mean_absolute_error(y_true, y_physics)
print(f"Physical model alone - R²: {physics_r2:.4f}, MAE: {physics_mae:.4f}")
# ML correction model training
print("\n=== ML Correction Model Training ===")
hybrid_model.train_correction(X, y_true, y_physics)
# Hybrid model performance
y_hybrid = y_physics.copy()
for i in range(len(X)):
inputs = {
'CA_in': X[i, 0],
'T_in': X[i, 1],
'T_jacket': X[i, 2],
'catalyst_age': X[i, 3],
'F': X[i, 4]
}
result = hybrid_model.predict_hybrid(inputs)
y_hybrid[i] = result['conversion_corrected']
hybrid_r2 = r2_score(y_true, y_hybrid)
hybrid_mae = mean_absolute_error(y_true, y_hybrid)
print(f"\nHybrid model - R²: {hybrid_r2:.4f}, MAE: {hybrid_mae:.4f}")
# Improvement rate
improvement = (physics_mae - hybrid_mae) / physics_mae * 100
print(f"MAE improvement rate: {improvement:.1f}%")
# Prediction comparison for new conditions
print("\n=== Prediction for New Conditions ===")
test_inputs = {
'CA_in': 1000,
'T_in': 298,
'T_jacket': 300,
'catalyst_age': 180, # 6 months old catalyst
'F': 0.1
}
result = hybrid_model.predict_hybrid(test_inputs)
print(f"Physical model prediction: {result['conversion_physics']*100:.1f}%")
print(f"Hybrid prediction: {result['conversion_corrected']*100:.1f}%")
print(f"Correction amount: {(result['conversion_corrected'] - result['conversion_physics'])*100:.1f}%")
# Expected output example:
# === Training Data Generation ===
# Physical model alone - R²: 0.8234, MAE: 0.0412
#
# === ML Correction Model Training ===
# Correction model training complete - R²: 0.8923, MAE: 0.0098
#
# Hybrid model - R²: 0.9756, MAE: 0.0123
# MAE improvement rate: 70.1%
#
# === Prediction for New Conditions ===
# Physical model prediction: 75.5%
# Hybrid prediction: 70.1%
# Correction amount: -5.4%
We optimize (calibrate) physical model parameters using measured data.
"""
Example 4: Model Calibration
Optimizing physical model parameters from measured data
"""
from scipy.optimize import minimize, differential_evolution
import numpy as np
from typing import List, Tuple
class ModelCalibrator:
"""Physical model parameter calibration
Search for parameters that minimize error with measured data
"""
def __init__(self, reactor_model: FirstPrinciplesReactor):
self.reactor_model = reactor_model
self.original_params = reactor_model.params
def objective_function(self, param_values: np.ndarray,
measured_data: List[Dict]) -> float:
"""Objective function for optimization (error with measurements)
Args:
param_values: Parameters to adjust [k0, Ea, UA]
measured_data: List of measured data
[{'CA_in': 1000, 'T_in': 298, 'T_jacket': 300,
'measured_conversion': 0.75, 'measured_T': 348}, ...]
Returns:
Mean squared error (MSE)
"""
# Update parameters
k0, Ea, UA = param_values
self.reactor_model.params.k0 = k0
self.reactor_model.params.Ea = Ea
self.reactor_model.params.UA = UA
total_error = 0
for data in measured_data:
# Predict with model
t, state = self.reactor_model.simulate(
initial_state=(0, data['T_in']),
time_span=(0, 60),
CA_in=data['CA_in'],
T_in=data['T_in'],
T_jacket=data['T_jacket']
)
CA_pred = state[-1, 0]
T_pred = state[-1, 1]
# Conversion error
conversion_pred = self.reactor_model.calculate_conversion(
CA_pred, data['CA_in']
)
error_conversion = (conversion_pred - data['measured_conversion']) ** 2
# Temperature error (weighted: 0.1)
error_T = 0.1 * ((T_pred - data['measured_T']) / 100) ** 2
total_error += error_conversion + error_T
mse = total_error / len(measured_data)
return mse
def calibrate(self, measured_data: List[Dict],
method: str = 'differential_evolution') -> Dict:
"""Execute parameter calibration
Args:
measured_data: List of measured data
method: Optimization method ('L-BFGS-B' or 'differential_evolution')
Returns:
Optimization results {'k0', 'Ea', 'UA', 'mse'}
"""
# Initial guess and search range
initial_guess = [
self.original_params.k0,
self.original_params.Ea,
self.original_params.UA
]
bounds = [
(1e8, 1e12), # k0 range
(60000, 100000), # Ea range [J/mol]
(3000, 7000) # UA range [W/K]
]
if method == 'differential_evolution':
# Global optimization (recommended)
result = differential_evolution(
lambda x: self.objective_function(x, measured_data),
bounds=bounds,
seed=42,
maxiter=100,
popsize=15,
tol=1e-6
)
optimal_params = result.x
mse = result.fun
elif method == 'L-BFGS-B':
# Local optimization
result = minimize(
lambda x: self.objective_function(x, measured_data),
x0=initial_guess,
bounds=bounds,
method='L-BFGS-B'
)
optimal_params = result.x
mse = result.fun
else:
raise ValueError(f"Unknown method: {method}")
return {
'k0': optimal_params[0],
'Ea': optimal_params[1],
'UA': optimal_params[2],
'mse': mse
}
def generate_measured_data(n_points: int = 20) -> List[Dict]:
"""Generate measured data (simulated)
Assumes data obtained from actual plant
Returns:
List of measured data
"""
np.random.seed(42)
# True parameters (assumed unknown)
true_params = ReactorParameters(
k0=7.5e9, # Smaller than initial guess of 1e10
Ea=85000, # Larger than initial guess of 80000
UA=4500 # Smaller than initial guess of 5000
)
true_model = FirstPrinciplesReactor(true_params)
measured_data = []
for _ in range(n_points):
# Randomly vary operating conditions
CA_in = np.random.uniform(800, 1200)
T_in = 298
T_jacket = np.random.uniform(290, 310)
# Simulate with true model
t, state = true_model.simulate(
initial_state=(0, T_in),
time_span=(0, 60),
CA_in=CA_in,
T_in=T_in,
T_jacket=T_jacket
)
CA = state[-1, 0]
T = state[-1, 1]
conversion = true_model.calculate_conversion(CA, CA_in)
# Add measurement noise
measured_data.append({
'CA_in': CA_in,
'T_in': T_in,
'T_jacket': T_jacket,
'measured_conversion': conversion + np.random.normal(0, 0.01),
'measured_T': T + np.random.normal(0, 1.0)
})
return measured_data
# Usage example
if __name__ == "__main__":
# Initial parameters (inaccurate)
initial_params = ReactorParameters(
k0=1e10, # Larger than true value
Ea=80000, # Smaller than true value
UA=5000 # Larger than true value
)
reactor = FirstPrinciplesReactor(initial_params)
# Generate measured data
print("=== Generate Measured Data ===")
measured_data = generate_measured_data(n_points=20)
print(f"Number of data points: {len(measured_data)}")
# Performance evaluation before calibration
print("\n=== Before Calibration ===")
calibrator = ModelCalibrator(reactor)
mse_before = calibrator.objective_function(
[initial_params.k0, initial_params.Ea, initial_params.UA],
measured_data
)
print(f"MSE: {mse_before:.6f}")
# Execute calibration
print("\n=== Executing Calibration... ===")
optimal_result = calibrator.calibrate(
measured_data,
method='differential_evolution'
)
print("\n=== Optimization Results ===")
print(f"Optimal k0: {optimal_result['k0']:.2e} (Initial: {initial_params.k0:.2e})")
print(f"Optimal Ea: {optimal_result['Ea']:.0f} J/mol (Initial: {initial_params.Ea:.0f})")
print(f"Optimal UA: {optimal_result['UA']:.0f} W/K (Initial: {initial_params.UA:.0f})")
print(f"Minimum MSE: {optimal_result['mse']:.6f}")
# Improvement rate
improvement = (mse_before - optimal_result['mse']) / mse_before * 100
print(f"\nMSE improvement rate: {improvement:.1f}%")
# Update model with optimal parameters
reactor.params.k0 = optimal_result['k0']
reactor.params.Ea = optimal_result['Ea']
reactor.params.UA = optimal_result['UA']
# Validation with test data
print("\n=== Validation with Test Data ===")
test_data = measured_data[0]
t, state = reactor.simulate(
initial_state=(0, test_data['T_in']),
time_span=(0, 60),
CA_in=test_data['CA_in'],
T_in=test_data['T_in'],
T_jacket=test_data['T_jacket']
)
CA_pred = state[-1, 0]
conversion_pred = reactor.calculate_conversion(CA_pred, test_data['CA_in'])
print(f"Predicted conversion: {conversion_pred*100:.1f}%")
print(f"Measured conversion: {test_data['measured_conversion']*100:.1f}%")
print(f"Error: {abs(conversion_pred - test_data['measured_conversion'])*100:.1f}%")
# Expected output example:
# === Generate Measured Data ===
# Number of data points: 20
#
# === Before Calibration ===
# MSE: 0.012345
#
# === Executing Calibration... ===
#
# === Optimization Results ===
# Optimal k0: 7.48e+09 (Initial: 1.00e+10)
# Optimal Ea: 85023 J/mol (Initial: 80000)
# Optimal UA: 4512 W/K (Initial: 5000)
# Minimum MSE: 0.000123
#
# MSE improvement rate: 99.0%
#
# === Validation with Test Data ===
# Predicted conversion: 74.3%
# Measured conversion: 74.5%
# Error: 0.2%
An efficient hybrid method that learns the prediction residuals (errors) of the physical model with machine learning.
"""
Example 5: Residual Learning Approach
Learning error patterns of physical models with machine learning
"""
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel
from typing import Dict, Tuple
class ResidualLearningModel:
"""Residual learning hybrid model
Learns the difference (residual) between physical model prediction and measured values
with Gaussian process regression, and also quantifies uncertainty
"""
def __init__(self, physics_model: FirstPrinciplesReactor):
self.physics_model = physics_model
# Gaussian process regression (also estimates uncertainty)
kernel = RBF(length_scale=1.0) + WhiteKernel(noise_level=1e-5)
self.gp_model = GaussianProcessRegressor(
kernel=kernel,
n_restarts_optimizer=10,
random_state=42
)
self.is_trained = False
def compute_physics_residuals(self, X: np.ndarray,
y_true: np.ndarray) -> np.ndarray:
"""Calculate physical model residuals
Args:
X: Features [n_samples, n_features]
columns: [CA_in, T_in, T_jacket, catalyst_age]
y_true: Measured conversion [n_samples]
Returns:
Residuals [n_samples]
"""
n_samples = X.shape[0]
residuals = np.zeros(n_samples)
for i in range(n_samples):
CA_in, T_in, T_jacket, _ = X[i]
# Predict with physical model
t, state = self.physics_model.simulate(
initial_state=(0, T_in),
time_span=(0, 60),
CA_in=CA_in,
T_in=T_in,
T_jacket=T_jacket
)
CA_pred = state[-1, 0]
conversion_pred = self.physics_model.calculate_conversion(
CA_pred, CA_in
)
# Residual = measured - physical prediction
residuals[i] = y_true[i] - conversion_pred
return residuals
def train(self, X: np.ndarray, y_true: np.ndarray):
"""Train residual learning model
Args:
X: Features [n_samples, n_features]
y_true: Measured conversion [n_samples]
"""
print("Calculating physical model residuals...")
residuals = self.compute_physics_residuals(X, y_true)
print(f"Residual statistics: Mean={np.mean(residuals):.4f}, Std={np.std(residuals):.4f}")
# Learn residuals with Gaussian process regression
print("Training Gaussian process regression...")
self.gp_model.fit(X, residuals)
self.is_trained = True
print("Training complete")
def predict(self, X: np.ndarray,
return_std: bool = False) -> Tuple[np.ndarray, np.ndarray]:
"""Execute prediction (physics + residual correction)
Args:
X: Features [n_samples, n_features]
return_std: Also return standard deviation (uncertainty)
Returns:
(predictions, std) Predicted values and standard deviation
"""
n_samples = X.shape[0]
physics_predictions = np.zeros(n_samples)
# Predict with physical model
for i in range(n_samples):
CA_in, T_in, T_jacket, _ = X[i]
t, state = self.physics_model.simulate(
initial_state=(0, T_in),
time_span=(0, 60),
CA_in=CA_in,
T_in=T_in,
T_jacket=T_jacket
)
CA_pred = state[-1, 0]
conversion_pred = self.physics_model.calculate_conversion(
CA_pred, CA_in
)
physics_predictions[i] = conversion_pred
# Predict residuals
if self.is_trained:
if return_std:
residual_pred, residual_std = self.gp_model.predict(
X, return_std=True
)
hybrid_predictions = physics_predictions + residual_pred
return hybrid_predictions, residual_std
else:
residual_pred = self.gp_model.predict(X)
hybrid_predictions = physics_predictions + residual_pred
return hybrid_predictions, None
else:
if return_std:
return physics_predictions, np.zeros(n_samples)
else:
return physics_predictions, None
# Usage example
if __name__ == "__main__":
# Physical model initialization
params = ReactorParameters()
physics_model = FirstPrinciplesReactor(params)
# Residual learning model initialization
residual_model = ResidualLearningModel(physics_model)
# Training data generation
print("=== Training Data Generation ===")
np.random.seed(42)
n_train = 50
X_train = np.column_stack([
np.random.uniform(800, 1200, n_train), # CA_in
np.full(n_train, 298), # T_in
np.random.uniform(290, 310, n_train), # T_jacket
np.random.uniform(0, 365, n_train) # catalyst_age
])
# True conversion (physical model + unknown effects)
y_train_physics = np.zeros(n_train)
for i in range(n_train):
t, state = physics_model.simulate(
initial_state=(0, X_train[i, 1]),
time_span=(0, 60),
CA_in=X_train[i, 0],
T_in=X_train[i, 1],
T_jacket=X_train[i, 2]
)
CA = state[-1, 0]
y_train_physics[i] = physics_model.calculate_conversion(CA, X_train[i, 0])
# Add unknown effects (catalyst degradation, etc.)
catalyst_effect = -0.0003 * X_train[:, 3]
y_train_true = y_train_physics + catalyst_effect + np.random.normal(0, 0.01, n_train)
# Training
print("\n=== Residual Learning Model Training ===")
residual_model.train(X_train, y_train_true)
# Evaluation with test data
print("\n=== Evaluation with Test Data ===")
n_test = 10
X_test = np.column_stack([
np.random.uniform(800, 1200, n_test),
np.full(n_test, 298),
np.random.uniform(290, 310, n_test),
np.random.uniform(0, 365, n_test)
])
# True values (test)
y_test_physics = np.zeros(n_test)
for i in range(n_test):
t, state = physics_model.simulate(
initial_state=(0, X_test[i, 1]),
time_span=(0, 60),
CA_in=X_test[i, 0],
T_in=X_test[i, 1],
T_jacket=X_test[i, 2]
)
CA = state[-1, 0]
y_test_physics[i] = physics_model.calculate_conversion(CA, X_test[i, 0])
catalyst_effect_test = -0.0003 * X_test[:, 3]
y_test_true = y_test_physics + catalyst_effect_test
# Prediction (with uncertainty)
y_pred, y_std = residual_model.predict(X_test, return_std=True)
# Display results
print("\nPrediction Results:")
print("Index\tTrue\tPhysics\tHybrid\tStd\tError")
print("-" * 60)
for i in range(n_test):
error = abs(y_test_true[i] - y_pred[i])
print(f"{i}\t{y_test_true[i]:.3f}\t{y_test_physics[i]:.3f}\t"
f"{y_pred[i]:.3f}\t{y_std[i]:.3f}\t{error:.3f}")
# Performance evaluation
mae_physics = np.mean(np.abs(y_test_true - y_test_physics))
mae_hybrid = np.mean(np.abs(y_test_true - y_pred))
print(f"\n=== Performance Comparison ===")
print(f"Physical model MAE: {mae_physics:.4f}")
print(f"Hybrid MAE: {mae_hybrid:.4f}")
print(f"Improvement rate: {(mae_physics - mae_hybrid) / mae_physics * 100:.1f}%")
# Expected output example:
# === Training Data Generation ===
#
# === Residual Learning Model Training ===
# Calculating physical model residuals...
# Residual statistics: Mean=-0.0534, Std=0.0456
# Training Gaussian process regression...
# Training complete
#
# === Evaluation with Test Data ===
#
# Prediction Results:
# Index True Physics Hybrid Std Error
# ------------------------------------------------------------
# 0 0.712 0.756 0.708 0.012 0.004
# 1 0.745 0.789 0.742 0.015 0.003
# 2 0.698 0.745 0.695 0.018 0.003
# 3 0.723 0.768 0.721 0.013 0.002
# 4 0.734 0.779 0.732 0.014 0.002
# 5 0.689 0.738 0.687 0.019 0.002
# 6 0.756 0.798 0.754 0.011 0.002
# 7 0.701 0.748 0.699 0.016 0.002
# 8 0.745 0.787 0.743 0.013 0.002
# 9 0.712 0.759 0.710 0.015 0.002
#
# === Performance Comparison ===
# Physical model MAE: 0.0456
# Hybrid MAE: 0.0025
# Improvement rate: 94.5%
Basic implementation of PINN, a neural network that incorporates physical laws into the loss function.
"""
Example 6: Physics-Informed Neural Networks (PINN) Basics
Neural networks incorporating physical laws as constraints
"""
import torch
import torch.nn as nn
import numpy as np
from typing import Tuple
class PhysicsInformedNN(nn.Module):
"""Physics-Informed Neural Network (PINN)
Incorporates residuals of physical laws (differential equations)
into the loss function of neural networks
"""
def __init__(self, input_dim: int = 4, hidden_dim: int = 64):
"""
Args:
input_dim: Input dimension (CA_in, T_in, T_jacket, t)
hidden_dim: Hidden layer dimension
"""
super(PhysicsInformedNN, self).__init__()
self.network = nn.Sequential(
nn.Linear(input_dim, hidden_dim),
nn.Tanh(),
nn.Linear(hidden_dim, hidden_dim),
nn.Tanh(),
nn.Linear(hidden_dim, hidden_dim),
nn.Tanh(),
nn.Linear(hidden_dim, 2) # Output [CA, T]
)
# Physical parameters (learnable)
self.k0 = nn.Parameter(torch.tensor(1e10, dtype=torch.float32))
self.Ea = nn.Parameter(torch.tensor(80000.0, dtype=torch.float32))
self.UA = nn.Parameter(torch.tensor(5000.0, dtype=torch.float32))
# Fixed parameters
self.V = 1.0
self.F = 0.1
self.rho = 1000.0
self.Cp = 4200.0
self.delta_H = -100000.0
self.R = 8.314
def forward(self, x: torch.Tensor) -> torch.Tensor:
"""Forward pass
Args:
x: [batch, 4] (CA_in, T_in, T_jacket, t)
Returns:
[batch, 2] (CA, T)
"""
return self.network(x)
def reaction_rate(self, CA: torch.Tensor, T: torch.Tensor) -> torch.Tensor:
"""Reaction rate (Arrhenius equation)
Args:
CA: Concentration [batch]
T: Temperature [batch]
Returns:
Reaction rate [batch]
"""
k = self.k0 * torch.exp(-self.Ea / (self.R * T))
r = k * CA
return r
def physics_loss(self, x: torch.Tensor,
CA: torch.Tensor, T: torch.Tensor) -> torch.Tensor:
"""Physical law residual (differential equation error)
Residual of dCA/dt = F/V * (CA_in - CA) - r
Residual of dT/dt = ...
Args:
x: Input [batch, 4]
CA: Concentration prediction [batch]
T: Temperature prediction [batch]
Returns:
Physical residual (scalar)
"""
CA_in = x[:, 0]
T_in = x[:, 1]
T_jacket = x[:, 2]
# Automatic differentiation to calculate dCA/dt, dT/dt
CA.requires_grad_(True)
T.requires_grad_(True)
# Reaction rate
r = self.reaction_rate(CA, T)
# Right-hand side of mass balance
dCA_dt_theory = (self.F / self.V) * (CA_in - CA) - r
# Right-hand side of energy balance
term1 = (self.F / self.V) * (T_in - T)
term2 = (-self.delta_H) * r / (self.rho * self.Cp)
term3 = -self.UA / (self.rho * self.Cp * self.V) * (T - T_jacket)
dT_dt_theory = term1 + term2 + term3
# Actual time derivative (differentiate NN output with respect to t)
# For simplicity, direct comparison with theoretical values here
# In practice, use torch.autograd.grad
physics_residual = torch.mean(dCA_dt_theory ** 2 + dT_dt_theory ** 2)
return physics_residual
def data_loss(self, CA_pred: torch.Tensor, T_pred: torch.Tensor,
CA_true: torch.Tensor, T_true: torch.Tensor) -> torch.Tensor:
"""Data error (difference from measured values)
Args:
CA_pred, T_pred: Predicted values
CA_true, T_true: Measured values
Returns:
Data loss (scalar)
"""
loss_CA = torch.mean((CA_pred - CA_true) ** 2)
loss_T = torch.mean((T_pred - T_true) ** 2) / 10000 # Temperature scale adjustment
return loss_CA + loss_T
def train_pinn(model: PhysicsInformedNN,
X_data: np.ndarray, CA_data: np.ndarray, T_data: np.ndarray,
epochs: int = 1000,
lambda_physics: float = 0.1) -> list:
"""PINN training
Args:
model: PINN model
X_data: Input data [n_samples, 4]
CA_data: Concentration measurements [n_samples]
T_data: Temperature measurements [n_samples]
epochs: Number of epochs
lambda_physics: Weight of physics loss
Returns:
Loss history
"""
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)
# Convert data to Tensor
X_tensor = torch.tensor(X_data, dtype=torch.float32)
CA_tensor = torch.tensor(CA_data, dtype=torch.float32)
T_tensor = torch.tensor(T_data, dtype=torch.float32)
loss_history = []
for epoch in range(epochs):
optimizer.zero_grad()
# Forward pass
output = model(X_tensor)
CA_pred = output[:, 0]
T_pred = output[:, 1]
# Data loss
loss_data = model.data_loss(CA_pred, T_pred, CA_tensor, T_tensor)
# Physics loss
loss_physics = model.physics_loss(X_tensor, CA_pred, T_pred)
# Total loss
total_loss = loss_data + lambda_physics * loss_physics
# Backward pass
total_loss.backward()
optimizer.step()
loss_history.append(total_loss.item())
if (epoch + 1) % 100 == 0:
print(f"Epoch {epoch+1}/{epochs} - "
f"Total Loss: {total_loss.item():.4f}, "
f"Data Loss: {loss_data.item():.4f}, "
f"Physics Loss: {loss_physics.item():.4f}")
return loss_history
# Usage example (simplified)
if __name__ == "__main__":
print("=== Physics-Informed Neural Network (PINN) ===")
# Training data generation (simplified)
np.random.seed(42)
n_samples = 100
X_data = np.column_stack([
np.random.uniform(800, 1200, n_samples), # CA_in
np.full(n_samples, 298), # T_in
np.random.uniform(290, 310, n_samples), # T_jacket
np.random.uniform(0, 60, n_samples) # t
])
# Simulated measurements
CA_data = np.random.uniform(200, 800, n_samples)
T_data = np.random.uniform(320, 360, n_samples)
# PINN model initialization
pinn_model = PhysicsInformedNN(input_dim=4, hidden_dim=64)
print(f"\nInitial parameters:")
print(f"k0: {pinn_model.k0.item():.2e}")
print(f"Ea: {pinn_model.Ea.item():.0f} J/mol")
print(f"UA: {pinn_model.UA.item():.0f} W/K")
# Training (simplified - more epochs needed in practice)
print(f"\nTraining started...")
loss_history = train_pinn(
pinn_model,
X_data, CA_data, T_data,
epochs=500,
lambda_physics=0.1
)
print(f"\nFinal parameters:")
print(f"k0: {pinn_model.k0.item():.2e}")
print(f"Ea: {pinn_model.Ea.item():.0f} J/mol")
print(f"UA: {pinn_model.UA.item():.0f} W/K")
print("\nPINN training complete (more detailed implementation needed for actual applications)")
# Expected output example:
# === Physics-Informed Neural Network (PINN) ===
#
# Initial parameters:
# k0: 1.00e+10
# Ea: 80000 J/mol
# UA: 5000 W/K
#
# Training started...
# Epoch 100/500 - Total Loss: 125432.5678, Data Loss: 123456.7890, Physics Loss: 1975.7788
# Epoch 200/500 - Total Loss: 45678.1234, Data Loss: 44567.2345, Physics Loss: 1110.8889
# Epoch 300/500 - Total Loss: 12345.6789, Data Loss: 11234.5678, Physics Loss: 1111.1111
# Epoch 400/500 - Total Loss: 3456.7890, Data Loss: 2345.6789, Physics Loss: 1111.1101
# Epoch 500/500 - Total Loss: 1234.5678, Data Loss: 123.4567, Physics Loss: 1111.1111
#
# Final parameters:
# k0: 8.23e+09
# Ea: 83215 J/mol
# UA: 4687 W/K
#
# PINN training complete (more detailed implementation needed for actual applications)
Automatically selects the optimal model (physical/data-driven/hybrid) according to operating conditions.
"""
Example 7: Model Selection and Switching Logic
Dynamically selecting optimal model based on operating conditions
"""
from enum import Enum
from typing import Dict, Optional
import numpy as np
class ModelType(Enum):
"""Model types"""
PHYSICS_ONLY = "physics"
DATA_DRIVEN = "data_driven"
HYBRID = "hybrid"
class AdaptiveModelSelector:
"""Adaptive model selector
Dynamically selects optimal model based on
operating conditions, data quality, and prediction accuracy
"""
def __init__(self,
physics_model: FirstPrinciplesReactor,
data_model: DataDrivenReactorModel,
hybrid_model: HybridReactorModel):
"""
Args:
physics_model: Physical model
data_model: Data-driven model
hybrid_model: Hybrid model
"""
self.physics_model = physics_model
self.data_model = data_model
self.hybrid_model = hybrid_model
# Application range for each model (training data range)
self.data_model_range = {
'T': (60, 95),
'P': (2.0, 3.0),
'F': (100, 150),
'catalyst_age': (0, 365)
}
# Model performance history
self.performance_history = {
ModelType.PHYSICS_ONLY: [],
ModelType.DATA_DRIVEN: [],
ModelType.HYBRID: []
}
def is_within_training_range(self, inputs: Dict) -> bool:
"""Check if inputs are within data-driven model training range
Args:
inputs: Input dictionary
Returns:
True if within range
"""
if 'T' in inputs:
T = inputs['T']
if not (self.data_model_range['T'][0] <= T <= self.data_model_range['T'][1]):
return False
if 'P' in inputs:
P = inputs['P']
if not (self.data_model_range['P'][0] <= P <= self.data_model_range['P'][1]):
return False
# Check other variables similarly
return True
def assess_data_quality(self, sensor_data: Dict) -> float:
"""Assess sensor data quality
Args:
sensor_data: Sensor data dictionary
Returns:
Quality score (0-1, 1 is highest)
"""
quality_score = 1.0
# Check data quality flag
if sensor_data.get('quality') == 'Bad':
quality_score *= 0.3
elif sensor_data.get('quality') == 'Uncertain':
quality_score *= 0.7
# Check data completeness
required_fields = ['temperature', 'pressure', 'flow_rate']
missing_count = sum(1 for field in required_fields if field not in sensor_data)
quality_score *= (1.0 - 0.2 * missing_count)
# Check noise level (if standard deviation is large)
if 'noise_level' in sensor_data:
quality_score *= max(0.5, 1.0 - sensor_data['noise_level'])
return quality_score
def select_model(self, inputs: Dict,
sensor_data: Dict,
data_quality: Optional[float] = None) -> ModelType:
"""Select optimal model based on operating conditions
Args:
inputs: Operating conditions
sensor_data: Sensor data
data_quality: Data quality score (auto-calculated if None)
Returns:
Selected model type
"""
if data_quality is None:
data_quality = self.assess_data_quality(sensor_data)
in_training_range = self.is_within_training_range(inputs)
# Decision logic
if data_quality < 0.5:
# Low data quality → Use physical model
return ModelType.PHYSICS_ONLY
elif in_training_range and data_quality >= 0.8:
# Within training range and high quality data → Hybrid model (highest accuracy)
return ModelType.HYBRID
elif in_training_range:
# Within training range but medium quality → Data-driven model
return ModelType.DATA_DRIVEN
else:
# Outside training range → Physical model (high extrapolation capability)
return ModelType.PHYSICS_ONLY
def predict(self, inputs: Dict, sensor_data: Dict) -> Dict:
"""Execute prediction (automatically select optimal model)
Args:
inputs: Operating conditions
sensor_data: Sensor data
Returns:
Prediction result dictionary {'conversion', 'model_used', 'confidence'}
"""
# Model selection
selected_model = self.select_model(inputs, sensor_data)
# Data quality evaluation
data_quality = self.assess_data_quality(sensor_data)
# Predict with selected model
if selected_model == ModelType.PHYSICS_ONLY:
result = self.physics_model.predict_physics(inputs)
confidence = 0.7 # Physical model confidence
elif selected_model == ModelType.DATA_DRIVEN:
X = np.array([[
inputs.get('T', 80),
inputs.get('P', 2.5),
inputs.get('F', 120),
inputs.get('catalyst_age', 0)
]])
conversion_pred = self.data_model.predict(X)[0]
result = {'conversion': conversion_pred}
confidence = 0.9 * data_quality # Depends on data quality
elif selected_model == ModelType.HYBRID:
result = self.hybrid_model.predict_hybrid(inputs)
confidence = 0.95 * data_quality # Highest accuracy
else:
raise ValueError(f"Unknown model type: {selected_model}")
# Add information to result
result['model_used'] = selected_model.value
result['confidence'] = confidence
result['data_quality'] = data_quality
return result
def update_performance(self, model_type: ModelType,
prediction: float, actual: float):
"""Update model performance history (used for online learning)
Args:
model_type: Model type used
prediction: Predicted value
actual: Measured value
"""
error = abs(prediction - actual)
self.performance_history[model_type].append(error)
# Keep only latest 100
if len(self.performance_history[model_type]) > 100:
self.performance_history[model_type] = \
self.performance_history[model_type][-100:]
def get_model_performance_summary(self) -> Dict:
"""Get performance summary for each model
Returns:
{'physics': {'mae': 0.05, 'count': 120}, ...}
"""
summary = {}
for model_type, errors in self.performance_history.items():
if len(errors) > 0:
summary[model_type.value] = {
'mae': np.mean(errors),
'std': np.std(errors),
'count': len(errors)
}
return summary
# Usage example
if __name__ == "__main__":
# Initialize each model (from previous examples)
params = ReactorParameters()
physics_model = FirstPrinciplesReactor(params)
data_model = DataDrivenReactorModel(model_type='random_forest')
# (training omitted)
hybrid_model = HybridReactorModel(physical_params=params)
# (training omitted)
# Adaptive model selector initialization
model_selector = AdaptiveModelSelector(
physics_model=physics_model,
data_model=data_model,
hybrid_model=hybrid_model
)
# Scenario 1: Normal operation (high quality data, within training range)
print("=== Scenario 1: Normal Operation ===")
inputs1 = {
'CA_in': 1000,
'T_in': 298,
'T_jacket': 300,
'T': 80,
'P': 2.5,
'F': 120,
'catalyst_age': 100
}
sensor_data1 = {
'temperature': 80.0,
'pressure': 2.5,
'flow_rate': 120.0,
'quality': 'Good'
}
result1 = model_selector.predict(inputs1, sensor_data1)
print(f"Model used: {result1['model_used']}")
print(f"Confidence: {result1['confidence']:.2f}")
print(f"Data quality: {result1['data_quality']:.2f}")
# Scenario 2: Abnormal operation (low quality data)
print("\n=== Scenario 2: Abnormal Operation (Low Quality Data) ===")
sensor_data2 = {
'temperature': 80.0,
'quality': 'Bad', # Bad data
'noise_level': 0.8
}
result2 = model_selector.predict(inputs1, sensor_data2)
print(f"Model used: {result2['model_used']}")
print(f"Confidence: {result2['confidence']:.2f}")
print(f"Data quality: {result2['data_quality']:.2f}")
# Scenario 3: Outside training range (extrapolation)
print("\n=== Scenario 3: Outside Training Range ===")
inputs3 = {
'CA_in': 1000,
'T_in': 298,
'T_jacket': 300,
'T': 105, # Outside training range (max 95°C)
'P': 2.5,
'F': 120,
'catalyst_age': 100
}
sensor_data3 = {
'temperature': 105.0,
'pressure': 2.5,
'flow_rate': 120.0,
'quality': 'Good'
}
result3 = model_selector.predict(inputs3, sensor_data3)
print(f"Model used: {result3['model_used']}")
print(f"Confidence: {result3['confidence']:.2f}")
print(f"Reason: Physical model used due to being outside training range")
# Update performance history (simulated)
print("\n=== Model Performance Update ===")
model_selector.update_performance(ModelType.HYBRID, 0.75, 0.74)
model_selector.update_performance(ModelType.HYBRID, 0.76, 0.75)
model_selector.update_performance(ModelType.PHYSICS_ONLY, 0.70, 0.75)
summary = model_selector.get_model_performance_summary()
for model, stats in summary.items():
print(f"{model}: MAE={stats['mae']:.4f}, Times used={stats['count']}")
# Expected output example:
# === Scenario 1: Normal Operation ===
# Model used: hybrid
# Confidence: 0.95
# Data quality: 1.00
#
# === Scenario 2: Abnormal Operation (Low Quality Data) ===
# Model used: physics
# Confidence: 0.70
# Data quality: 0.06
#
# === Scenario 3: Outside Training Range ===
# Model used: physics
# Confidence: 0.70
# Reason: Physical model used due to being outside training range
#
# === Model Performance Update ===
# hybrid: MAE=0.0100, Times used=2
# physics: MAE=0.0500, Times used=1
Upon completing this chapter, you will be able to implement the following:
Q1: Which of the following is NOT an advantage of hybrid models?
a) Can achieve both physical interpretability and high accuracy
b) High performance even with limited training data
c) Easiest to implement
d) Achieves both extrapolation capability and data fitting
Correct answer: c) Easiest to implement
Explanation:
Hybrid models are the most complex to implement because they require integration of physical and data-driven models. However, in exchange for this complexity, the following advantages are obtained:
Q2: In physical model calibration, the parameter optimization objective function (MSE) improved from initial value 0.012345 to 0.000123. How much did RMSE improve? Also calculate the improvement rate.
Calculation:
Correct answer: RMSE improvement rate = 90.0%
Important points:
Q3: You have been assigned to develop a Digital Twin for a new chemical plant. Given the following conditions, propose which of first-principles model, data-driven model, or hybrid model should be adopted, along with your reasoning.
Conditions:
Recommendation: Phased Approach (First-Principles → Hybrid)
Phase 1 (First month): First-Principles Model Development
Phase 2 (Months 2-3): Hybridization
Reasons for not choosing data-driven alone:
Implementation Schedule (3 months):
| Month | Task | Milestone |
|---|---|---|
| 1 | First-principles model development, calibration with pilot data | ±5% accuracy achieved |
| 2 | Full-scale operation start, data collection, ML correction model development | 500 hours data accumulated |
| 3 | Hybrid model integration, validation, deployment | ±2% accuracy achieved |
Risk Mitigation Measures:
Having mastered the basics of hybrid modeling, next we move on to practical applications of Digital Twin. We will learn about real business use cases such as process optimization, predictive maintenance, and what-if analysis.