"""
Parameter Estimation for SPROCLIB
This module provides parameter estimation tools for process control
systems using experimental data and optimization techniques.
Author: Thorsten Gressling <gressling@paramus.ai>
License: MIT License
"""
import numpy as np
from typing import Optional, Tuple, Dict, Any, List, Callable
from scipy.optimize import minimize, differential_evolution
import logging
logger = logging.getLogger(__name__)
[docs]
class ParameterEstimation:
"""Parameter estimation for process models."""
[docs]
def __init__(self, name: str = "Parameter Estimation"):
"""
Initialize parameter estimation.
Args:
name: Estimation name
"""
self.name = name
self.results = {}
logger.info(f"Parameter estimation '{name}' initialized")
[docs]
def estimate_parameters(
self,
model_func: Callable,
t_data: np.ndarray,
y_data: np.ndarray,
param_bounds: List[Tuple[float, float]],
initial_guess: Optional[np.ndarray] = None,
method: str = 'least_squares'
) -> Dict[str, Any]:
"""
Estimate model parameters from experimental data.
Args:
model_func: Model function that takes (t, params) and returns y
t_data: Time data
y_data: Output data
param_bounds: Parameter bounds [(min, max), ...]
initial_guess: Initial parameter guess
method: Estimation method
Returns:
Parameter estimation results
"""
def objective(params):
"""Objective function for parameter estimation."""
try:
y_pred = model_func(t_data, params)
residuals = y_data - y_pred
return np.sum(residuals**2) # Sum of squared errors
except Exception as e:
logger.warning(f"Model evaluation error: {e}")
return 1e6 # Large penalty for invalid parameters
try:
if method.lower() == 'least_squares':
if initial_guess is None:
# Use midpoint of bounds as initial guess
initial_guess = np.array([(b[0] + b[1]) / 2 for b in param_bounds])
result = minimize(
objective, initial_guess, method='L-BFGS-B',
bounds=param_bounds
)
elif method.lower() == 'differential_evolution':
result = differential_evolution(
objective, param_bounds, seed=42
)
else:
raise ValueError(f"Unknown estimation method: {method}")
if result.success:
# Calculate model predictions and statistics
y_pred = model_func(t_data, result.x)
residuals = y_data - y_pred
mse = np.mean(residuals**2)
rmse = np.sqrt(mse)
# R-squared
ss_res = np.sum(residuals**2)
ss_tot = np.sum((y_data - np.mean(y_data))**2)
r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0
self.results = {
'success': True,
'parameters': result.x,
'objective_value': result.fun,
'y_predicted': y_pred,
'residuals': residuals,
'mse': mse,
'rmse': rmse,
'r_squared': r_squared,
'method': method,
'message': result.message
}
logger.info(f"Parameter estimation successful: RMSE = {rmse:.6f}, R² = {r_squared:.4f}")
else:
self.results = {
'success': False,
'message': result.message,
'method': method
}
logger.error(f"Parameter estimation failed: {result.message}")
return self.results
except Exception as e:
logger.error(f"Parameter estimation error: {e}")
return {
'success': False,
'error': str(e),
'method': method
}
[docs]
def validate_model(
self,
model_func: Callable,
parameters: np.ndarray,
t_validation: np.ndarray,
y_validation: np.ndarray
) -> Dict[str, Any]:
"""
Validate model with independent validation data.
Args:
model_func: Model function
parameters: Estimated parameters
t_validation: Validation time data
y_validation: Validation output data
Returns:
Validation results
"""
try:
# Generate predictions
y_pred = model_func(t_validation, parameters)
# Calculate validation metrics
residuals = y_validation - y_pred
mse = np.mean(residuals**2)
rmse = np.sqrt(mse)
mae = np.mean(np.abs(residuals))
# R-squared
ss_res = np.sum(residuals**2)
ss_tot = np.sum((y_validation - np.mean(y_validation))**2)
r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0
validation_results = {
'y_predicted': y_pred,
'residuals': residuals,
'mse': mse,
'rmse': rmse,
'mae': mae,
'r_squared': r_squared,
'validation_points': len(y_validation)
}
logger.info(f"Model validation: RMSE = {rmse:.6f}, R² = {r_squared:.4f}")
return validation_results
except Exception as e:
logger.error(f"Model validation error: {e}")
return {'error': str(e)}
[docs]
def describe(self) -> Dict[str, Any]:
"""
Provide comprehensive information about ParameterEstimation capabilities.
Returns:
Dictionary containing parameter estimation algorithms, applications,
parameters, and technical specifications for chemical engineering models.
"""
return {
'class': 'ParameterEstimation',
'purpose': 'Parameter estimation for chemical process models from experimental data',
'version': '1.0.0',
'algorithms': [
'Least Squares Optimization',
'Differential Evolution (Global)',
'L-BFGS-B (Bound Constrained)',
'Trust Region Methods',
'Genetic Algorithm',
'Particle Swarm Optimization'
],
'applications': [
'FOPDT (First Order Plus Dead Time) model identification',
'Reaction kinetics parameter estimation',
'Heat transfer coefficient determination',
'Mass transfer parameter identification',
'Process control model identification',
'Thermodynamic property estimation',
'Equipment performance parameter fitting',
'Catalyst activity parameter estimation'
],
'model_types': [
'First Order Plus Dead Time (FOPDT)',
'Second Order Plus Dead Time (SOPDT)',
'Arrhenius kinetics models',
'Heat transfer models',
'Mass transfer models',
'Adsorption isotherms',
'PID controller models',
'State-space models'
],
'parameters': {
'model_func': 'Mathematical model function f(t, params)',
't_data': 'Time vector for experimental data',
'y_data': 'Output measurements corresponding to time points',
'param_bounds': 'Lower and upper bounds for each parameter',
'initial_guess': 'Starting point for optimization (optional)',
'method': 'Optimization algorithm selection'
},
'estimation_methods': {
'least_squares': 'Minimize sum of squared residuals',
'differential_evolution': 'Global optimization for multimodal problems',
'maximum_likelihood': 'Statistical parameter estimation',
'bayesian': 'Probabilistic parameter estimation',
'robust': 'Outlier-resistant parameter estimation'
},
'objective_functions': [
'Sum of Squared Errors (SSE)',
'Mean Squared Error (MSE)',
'Root Mean Squared Error (RMSE)',
'Mean Absolute Error (MAE)',
'Maximum Likelihood',
'Weighted Least Squares'
],
'statistical_metrics': [
'R-squared (coefficient of determination)',
'Adjusted R-squared',
'Root Mean Squared Error (RMSE)',
'Mean Absolute Error (MAE)',
'Confidence intervals',
'Parameter covariance matrix',
'Residual analysis',
'F-statistic'
],
'data_requirements': {
'minimum_points': 'At least 3x number of parameters',
'time_span': 'Cover full dynamic response',
'sampling_rate': 'Adequate for system dynamics',
'signal_noise_ratio': 'SNR > 10 for reliable estimation',
'excitation': 'Rich input signals for identifiability'
},
'chemical_engineering_models': {
'reaction_kinetics': {
'arrhenius': 'k = A * exp(-Ea/(R*T))',
'power_law': 'r = k * C_A^n * C_B^m',
'langmuir_hinshelwood': 'Complex catalytic reactions'
},
'heat_transfer': {
'convection': 'h = Nu * k / L',
'radiation': 'q = σ * ε * A * (T1^4 - T2^4)',
'conduction': 'q = k * A * ΔT / L'
},
'mass_transfer': {
'film_theory': 'N = k_L * (C_i - C_bulk)',
'penetration_theory': 'k_L = 2 * sqrt(D / (π * t))',
'surface_renewal': 'k_L = sqrt(D * s)'
},
'thermodynamics': {
'antoine_equation': 'log(P) = A - B/(C + T)',
'riedel_equation': 'Extended vapor pressure correlation',
'wilson_equation': 'Activity coefficient model'
}
},
'validation_techniques': [
'Cross-validation',
'Hold-out validation',
'Bootstrap resampling',
'Residual analysis',
'Prediction confidence intervals',
'Parameter sensitivity analysis'
],
'optimization_features': {
'global_search': 'Find global minimum in parameter space',
'constraint_handling': 'Physical bounds and equality constraints',
'multi_objective': 'Balance fit quality vs model complexity',
'robustness': 'Handle outliers and measurement noise',
'convergence_criteria': 'Adaptive tolerance settings'
},
'typical_applications': {
'reactor_modeling': {
'parameters': ['reaction rate constants', 'activation energy', 'pre-exponential factor'],
'data_needed': 'Concentration vs time profiles',
'typical_ranges': 'k: 1e-6 to 1e6 1/s, Ea: 20-200 kJ/mol'
},
'heat_exchanger': {
'parameters': ['overall heat transfer coefficient', 'heat capacity', 'flow rates'],
'data_needed': 'Temperature profiles vs time',
'typical_ranges': 'U: 100-2000 W/m²/K'
},
'distillation_column': {
'parameters': ['tray efficiency', 'pressure drop', 'mass transfer coefficients'],
'data_needed': 'Composition and temperature profiles',
'typical_ranges': 'Efficiency: 0.3-0.9'
},
'process_control': {
'parameters': ['process gain', 'time constant', 'dead time'],
'data_needed': 'Step response or impulse response',
'typical_ranges': 'K: 0.1-10, τ: 1-1000 s, θ: 0-100 s'
}
},
'error_analysis': [
'Parameter uncertainty quantification',
'Confidence interval calculation',
'Sensitivity coefficient analysis',
'Correlation matrix computation',
'Model discrimination tests'
],
'computational_aspects': {
'convergence_criteria': 'Gradient norm < tolerance',
'numerical_derivatives': 'Finite difference approximations',
'scaling': 'Parameter normalization for better conditioning',
'initial_guess_strategies': 'Multiple starting points',
'regularization': 'Ridge regression for ill-conditioned problems'
},
'limitations': [
'Requires good initial parameter estimates for local methods',
'Model structure must be specified a priori',
'Identifiability issues with correlated parameters',
'Sensitive to measurement noise and outliers',
'May converge to local minima in complex parameter spaces',
'Requires sufficient data excitation for all parameters'
],
'performance_metrics': {
'small_problems': 'Sub-second estimation (< 10 parameters)',
'medium_problems': '1-10 seconds (10-50 parameters)',
'large_problems': '10-300 seconds (> 50 parameters)',
'accuracy': 'Parameter uncertainty typically 1-10%',
'robustness': 'Handles 5-20% measurement noise'
},
'integration_capabilities': [
'Process simulation software',
'Laboratory data management systems',
'Real-time optimization systems',
'Model predictive control',
'Statistical analysis packages'
]
}
# Standalone functions for backward compatibility
[docs]
def estimate_fopdt_parameters(
t_data: np.ndarray,
y_data: np.ndarray,
step_magnitude: float = 1.0
) -> Dict[str, Any]:
"""
Estimate FOPDT parameters from step response data.
Args:
t_data: Time data
y_data: Step response data
step_magnitude: Magnitude of step input
Returns:
Estimated FOPDT parameters
"""
def fopdt_model(t, params):
"""FOPDT model: y = K * (1 - exp(-(t-theta)/tau)) for t >= theta"""
K, tau, theta = params
y = np.zeros_like(t)
mask = t >= theta
y[mask] = K * (1 - np.exp(-(t[mask] - theta) / tau))
return y * step_magnitude
# Parameter bounds: K, tau, theta
param_bounds = [
(0.1, 10.0), # K: process gain
(0.01, max(t_data)), # tau: time constant
(0.0, max(t_data) / 2) # theta: dead time
]
estimator = ParameterEstimation("FOPDT Estimation")
result = estimator.estimate_parameters(
fopdt_model, t_data, y_data, param_bounds
)
if result.get('success', False):
K, tau, theta = result['parameters']
result.update({
'K': K,
'tau': tau,
'theta': theta
})
return result
__all__ = [
'ParameterEstimation',
'estimate_fopdt_parameters'
]