Implementing Economic Optimization and Model Predictive Control in Chemical Plants
What You Will Learn in This Chapter:
In chemical plant operations, Real-Time Optimization (RTO) and Advanced Process Control (APC) are essential technologies for achieving both economic performance and safety. This chapter covers optimization with SciPy and Pyomo, the implementation of Model Predictive Control (MPC), and next-generation process control using deep reinforcement learning (DQN, PPO), all at the implementation level.
In modern chemical plants, control systems are organized hierarchically. Each layer operates on a different time scale, with lower layers realizing the optimization goals set by the layers above them.
In a petroleum refinery (FCC unit):
This hierarchical structure has been shown to improve plant profitability by several hundred million yen annually in real-world deployments.
Real-time optimization takes the current plant state as input and computes optimal operating conditions that satisfy an economic objective function (e.g., profit maximization). We implement this using the operating optimization of a Continuous Stirred Tank Reactor (CSTR) as an example.
"""
===================================
Example 1: Online Optimization (SciPy)
===================================
Optimization of temperature and flow rate in a Continuous Stirred Tank Reactor (CSTR).
Maximizes profit per unit time by considering the trade-off between reaction rate and selectivity.
Objective: Maximize product value while minimizing energy and raw material costs
"""
import numpy as np
from scipy.optimize import minimize, NonlinearConstraint
from typing import Dict, Tuple
import pandas as pd
class CSTROptimizer:
"""Real-time optimization for a Continuous Stirred Tank Reactor"""
def __init__(self):
# Process parameters
self.volume = 10.0 # Reactor volume [m3]
self.heat_capacity = 4.18 # Heat capacity [kJ/kg*K]
self.density = 1000.0 # Density [kg/m3]
# Reaction rate constants (Arrhenius equation)
self.A1 = 1.2e10 # Pre-exponential factor (main reaction) [1/h]
self.E1 = 75000.0 # Activation energy (main reaction) [J/mol]
self.A2 = 3.5e9 # Pre-exponential factor (side reaction) [1/h]
self.E2 = 68000.0 # Activation energy (side reaction) [J/mol]
# Economic parameters
self.product_price = 150.0 # Product price [$/kg]
self.byproduct_price = 40.0 # Byproduct price [$/kg]
self.feed_cost = 50.0 # Feed cost [$/kg]
self.energy_cost = 0.08 # Energy cost [$/kWh]
# Physical constraints
self.T_min, self.T_max = 320.0, 380.0 # Temperature range [K]
self.F_min, self.F_max = 0.5, 5.0 # Flow rate range [m3/h]
self.T_feed = 298.0 # Feed temperature [K]
def reaction_rates(self, T: float) -> Tuple[float, float]:
"""Compute reaction rate constants (Arrhenius equation)
Args:
T: Reaction temperature [K]
Returns:
(main reaction rate constant, side reaction rate constant) [1/h]
"""
R = 8.314 # Gas constant [J/mol*K]
k1 = self.A1 * np.exp(-self.E1 / (R * T))
k2 = self.A2 * np.exp(-self.E2 / (R * T))
return k1, k2
def conversion_selectivity(self, T: float, tau: float) -> Tuple[float, float]:
"""Compute conversion and selectivity
Args:
T: Reaction temperature [K]
tau: Residence time [h]
Returns:
(conversion, selectivity)
"""
k1, k2 = self.reaction_rates(T)
# First-order reaction conversion
conversion = 1.0 - np.exp(-(k1 + k2) * tau)
# Selectivity (main product / total product)
selectivity = k1 / (k1 + k2)
return conversion, selectivity
def heating_power(self, T: float, F: float) -> float:
"""Compute the power required for heating
Args:
T: Reaction temperature [K]
F: Flow rate [m3/h]
Returns:
Heating power [kW]
"""
delta_T = T - self.T_feed
mass_flow = F * self.density # [kg/h]
heat_duty = mass_flow * self.heat_capacity * delta_T # [kJ/h]
return heat_duty / 3600.0 # [kW]
def objective(self, x: np.ndarray) -> float:
"""Objective function: negative profit (converted to a minimization problem)
Args:
x: [temperature [K], flow rate [m3/h]]
Returns:
-profit [$/h]
"""
T, F = x
tau = self.volume / F # Residence time [h]
# Conversion and selectivity
conversion, selectivity = self.conversion_selectivity(T, tau)
# Product formation rate
feed_mass = F * self.density # [kg/h]
product_mass = feed_mass * conversion * selectivity
byproduct_mass = feed_mass * conversion * (1 - selectivity)
# Revenue
revenue = (product_mass * self.product_price +
byproduct_mass * self.byproduct_price)
# Cost
feed_cost_total = feed_mass * self.feed_cost
energy_cost_total = self.heating_power(T, F) * self.energy_cost
profit = revenue - feed_cost_total - energy_cost_total
return -profit # Negated for minimization
def optimize(self, initial_guess: np.ndarray = None) -> Dict:
"""Run the optimization
Args:
initial_guess: Initial guess [temperature, flow rate]
Returns:
Dictionary of optimization results
"""
if initial_guess is None:
initial_guess = np.array([350.0, 2.0]) # [K, m3/h]
# Bound constraints
bounds = [(self.T_min, self.T_max),
(self.F_min, self.F_max)]
# Nonlinear constraint: conversion must be at least 0.85 (safety/quality requirement)
def conversion_constraint(x):
T, F = x
tau = self.volume / F
conversion, _ = self.conversion_selectivity(T, tau)
return conversion - 0.85
nlc = NonlinearConstraint(conversion_constraint, 0, np.inf)
# Run optimization
result = minimize(
self.objective,
initial_guess,
method='SLSQP',
bounds=bounds,
constraints=[nlc],
options={'ftol': 1e-6, 'disp': False}
)
# Organize the results
T_opt, F_opt = result.x
tau_opt = self.volume / F_opt
conversion, selectivity = self.conversion_selectivity(T_opt, tau_opt)
return {
'success': result.success,
'temperature': T_opt,
'flow_rate': F_opt,
'residence_time': tau_opt,
'conversion': conversion,
'selectivity': selectivity,
'profit_per_hour': -result.fun,
'heating_power': self.heating_power(T_opt, F_opt)
}
# ===================================
# Example run
# ===================================
if __name__ == "__main__":
optimizer = CSTROptimizer()
print("="*70)
print("CSTR Real-Time Optimization")
print("="*70)
# Run optimization
result = optimizer.optimize()
if result['success']:
print("\n[Optimal Operating Conditions]")
print(f" Reaction temperature: {result['temperature']:.1f} K ({result['temperature']-273.15:.1f} deg C)")
print(f" Flow rate: {result['flow_rate']:.2f} m3/h")
print(f" Residence time: {result['residence_time']:.2f} h")
print(f"\n[Process Performance]")
print(f" Conversion: {result['conversion']:.1%}")
print(f" Selectivity: {result['selectivity']:.1%}")
print(f" Heating power: {result['heating_power']:.1f} kW")
print(f"\n[Economics]")
print(f" Profit: ${result['profit_per_hour']:.2f}/h")
print(f" Annual profit: ${result['profit_per_hour'] * 8760:.0f}/year")
else:
print("Optimization failed.")
# Sensitivity analysis: impact of feed cost variation on the optimal conditions
print("\n" + "="*70)
print("Sensitivity Analysis: Effect of Feed Cost")
print("="*70)
feed_costs = [40, 50, 60, 70]
results = []
for cost in feed_costs:
optimizer.feed_cost = cost
res = optimizer.optimize()
results.append({
'Feed cost [$/kg]': cost,
'Optimal temperature [K]': res['temperature'],
'Optimal flow rate [m3/h]': res['flow_rate'],
'Hourly profit [$/h]': res['profit_per_hour']
})
df = pd.DataFrame(results)
print(df.to_string(index=False))
print("\n[OK] As feed cost rises, the optimum shifts to higher temperature and lower flow rate (favoring selectivity)")
For more complex optimization problems, we use Pyomo, an algebraic modeling language. We implement a formulation that simultaneously considers maximizing product value and minimizing utility costs.
"""
===================================
Example 2: Economic Optimization (Pyomo)
===================================
Economic optimization of a chemical plant using Pyomo.
For a process with multiple products, maximizes product value while
minimizing utility costs (steam, electricity, cooling water).
Pyomo requires a separate installation in practice; here we show a conceptual implementation.
"""
import numpy as np
from scipy.optimize import minimize
from typing import Dict, List
import pandas as pd
class EconomicOptimizer:
"""Economic optimization system for a chemical plant"""
def __init__(self):
# Product prices [$/ton]
self.product_prices = {
'ProductA': 800.0, # High value-added product
'ProductB': 500.0, # Intermediate product
'ProductC': 300.0 # Commodity product
}
# Utility costs
self.steam_cost = 25.0 # Steam [$/ton]
self.power_cost = 0.10 # Electricity [$/kWh]
self.cooling_cost = 0.5 # Cooling water [$/ton]
# Process constraint parameters
self.max_capacity = 100.0 # Maximum processing capacity [ton/h]
self.min_turndown = 0.4 # Minimum turndown ratio
def production_model(self, feed_rate: float, temperature: float,
pressure: float) -> Dict[str, float]:
"""Production model: compute product yields from operating conditions
Args:
feed_rate: Feed flow rate [ton/h]
temperature: Reaction temperature [K]
pressure: Reaction pressure [bar]
Returns:
Production rate of each product [ton/h]
"""
# Simplified yield model (a detailed reaction model would be used in practice)
T_ref = 400.0 # Reference temperature [K]
P_ref = 20.0 # Reference pressure [bar]
# Temperature and pressure effect factors
temp_factor = np.exp(-0.005 * (temperature - T_ref)**2)
press_factor = 1.0 + 0.02 * (pressure - P_ref)
# Base yield of each product
yield_A_base = 0.35 * temp_factor * press_factor
yield_B_base = 0.45 * (2.0 - temp_factor)
yield_C_base = 0.20
# Yield constraint (total <= 1.0)
total_yield = yield_A_base + yield_B_base + yield_C_base
if total_yield > 1.0:
scale = 1.0 / total_yield
yield_A_base *= scale
yield_B_base *= scale
yield_C_base *= scale
return {
'ProductA': feed_rate * yield_A_base,
'ProductB': feed_rate * yield_B_base,
'ProductC': feed_rate * yield_C_base
}
def utility_consumption(self, feed_rate: float, temperature: float,
pressure: float) -> Dict[str, float]:
"""Compute utility consumption
Args:
feed_rate: Feed flow rate [ton/h]
temperature: Reaction temperature [K]
pressure: Reaction pressure [bar]
Returns:
Utility consumption
"""
# Steam consumption (proportional to heating load)
T_feed = 298.0 # Feed temperature [K]
heating_load = feed_rate * 2.5 * (temperature - T_feed) # Simplified formula
steam = heating_load / 2000.0 # [ton/h]
# Electricity consumption (compressor, pumps, agitator)
compressor_power = 50.0 * (pressure / 20.0)**0.8 # [kW]
pump_power = 10.0 * feed_rate
agitator_power = 15.0
power = compressor_power + pump_power + agitator_power
# Cooling water (heat of reaction removal)
exothermic_heat = feed_rate * 500.0 # [kW] (assumed)
cooling_water = exothermic_heat / 40.0 # [ton/h]
return {
'steam': steam,
'power': power,
'cooling': cooling_water
}
def objective_function(self, x: np.ndarray) -> float:
"""Objective function: negative profit (minimization problem)
Args:
x: [feed_rate, temperature, pressure]
Returns:
-profit [$/h]
"""
feed_rate, temperature, pressure = x
# Product production rates
products = self.production_model(feed_rate, temperature, pressure)
# Revenue
revenue = sum(
products[prod] * self.product_prices[prod]
for prod in products
)
# Utility consumption
utilities = self.utility_consumption(feed_rate, temperature, pressure)
# Cost
utility_cost = (
utilities['steam'] * self.steam_cost +
utilities['power'] * self.power_cost +
utilities['cooling'] * self.cooling_cost
)
# Feed cost (assumed: $200/ton)
feed_cost = feed_rate * 200.0
profit = revenue - utility_cost - feed_cost
return -profit
def optimize_economics(self) -> Dict:
"""Run economic optimization
Returns:
Optimization results
"""
# Initial guess: [feed_rate, temperature, pressure]
x0 = np.array([60.0, 400.0, 20.0])
# Bound constraints
bounds = [
(self.max_capacity * self.min_turndown, self.max_capacity), # feed_rate
(350.0, 450.0), # temperature [K]
(10.0, 40.0) # pressure [bar]
]
# Optimization
result = minimize(
self.objective_function,
x0,
method='L-BFGS-B',
bounds=bounds,
options={'ftol': 1e-6}
)
feed_opt, temp_opt, press_opt = result.x
# Organize the results
products = self.production_model(feed_opt, temp_opt, press_opt)
utilities = self.utility_consumption(feed_opt, temp_opt, press_opt)
return {
'success': result.success,
'feed_rate': feed_opt,
'temperature': temp_opt,
'pressure': press_opt,
'products': products,
'utilities': utilities,
'profit_per_hour': -result.fun
}
# ===================================
# Example run
# ===================================
if __name__ == "__main__":
optimizer = EconomicOptimizer()
print("="*70)
print("Chemical Plant Economic Optimization")
print("="*70)
result = optimizer.optimize_economics()
if result['success']:
print("\n[Optimal Operating Conditions]")
print(f" Feed rate: {result['feed_rate']:.1f} ton/h")
print(f" Reaction temperature: {result['temperature']:.1f} K ({result['temperature']-273.15:.1f} deg C)")
print(f" Reaction pressure: {result['pressure']:.1f} bar")
print("\n[Product Production]")
for prod, amount in result['products'].items():
price = optimizer.product_prices[prod]
value = amount * price
print(f" {prod}: {amount:.2f} ton/h (value: ${value:.2f}/h)")
print("\n[Utility Consumption]")
util = result['utilities']
print(f" Steam: {util['steam']:.2f} ton/h (${util['steam'] * optimizer.steam_cost:.2f}/h)")
print(f" Power: {util['power']:.1f} kW (${util['power'] * optimizer.power_cost:.2f}/h)")
print(f" Cooling water: {util['cooling']:.2f} ton/h (${util['cooling'] * optimizer.cooling_cost:.2f}/h)")
print("\n[Economics]")
print(f" Hourly profit: ${result['profit_per_hour']:.2f}/h")
print(f" Daily profit: ${result['profit_per_hour'] * 24:.2f}/day")
print(f" Annual profit: ${result['profit_per_hour'] * 8760:.0f}/year")
# Sensitivity analysis: impact of product price variation
print("\n" + "="*70)
print("Sensitivity Analysis: Effect of Product A Price")
print("="*70)
original_price = optimizer.product_prices['ProductA']
price_scenarios = [600, 700, 800, 900, 1000]
results_table = []
for price in price_scenarios:
optimizer.product_prices['ProductA'] = price
res = optimizer.optimize_economics()
results_table.append({
'Product A price [$/ton]': price,
'Optimal feed rate [ton/h]': res['feed_rate'],
'Optimal temperature [K]': res['temperature'],
'Product A output [ton/h]': res['products']['ProductA'],
'Hourly profit [$/h]': res['profit_per_hour']
})
df = pd.DataFrame(results_table)
print(df.to_string(index=False))
print("\n[OK] As the price of the high value-added product rises, higher-temperature operation maximizes yield")
Model Predictive Control (MPC) is an advanced control technique that predicts future process behavior using a dynamic model and optimizes the manipulated variables while satisfying constraints. We show a basic implementation using distillation column temperature control as an example.
"""
===================================
Example 3: Basic MPC Implementation
===================================
Applying Model Predictive Control (MPC) to temperature control in a distillation column.
Optimizes reflux ratio and reflux flow rate while predicting future behavior.
Advantages of MPC:
- Multivariable control (multiple manipulated and controlled variables)
- Explicit handling of constraints
- Prediction and compensation of disturbances
"""
import numpy as np
from scipy.optimize import minimize
from typing import List, Tuple
import matplotlib.pyplot as plt
class DistillationMPC:
"""Model predictive control for a distillation column"""
def __init__(self, prediction_horizon: int = 10, control_horizon: int = 5):
"""
Args:
prediction_horizon: Prediction horizon (number of steps)
control_horizon: Control horizon (number of steps)
"""
self.Np = prediction_horizon
self.Nc = control_horizon
self.dt = 1.0 # Sampling time [min]
# Process model (state-space model)
# x[k+1] = A*x[k] + B*u[k]
# y[k] = C*x[k]
# State: x = [top temperature deviation, bottom temperature deviation]
# Input: u = [reflux ratio change, reflux flow rate change]
# Output: y = [top temperature, bottom temperature]
self.A = np.array([
[0.85, 0.10],
[0.05, 0.90]
])
self.B = np.array([
[0.5, 0.1],
[0.1, 0.4]
])
self.C = np.eye(2)
# Constraints
self.u_min = np.array([-2.0, -5.0]) # [reflux ratio, reflux flow rate kg/min]
self.u_max = np.array([2.0, 5.0])
self.delta_u_max = np.array([0.5, 1.0]) # Rate-of-change constraint
# Weighting matrices
self.Q = np.diag([10.0, 8.0]) # Output tracking weight
self.R = np.diag([1.0, 1.0]) # Input change weight
def predict(self, x0: np.ndarray, u_sequence: np.ndarray) -> np.ndarray:
"""Predict the output over the prediction horizon
Args:
x0: Current state
u_sequence: Input sequence (Nc x 2)
Returns:
Predicted output sequence (Np x 2)
"""
x = x0.copy()
y_pred = np.zeros((self.Np, 2))
for k in range(self.Np):
# Within the control horizon these are optimization variables; beyond it, hold the last value
if k < self.Nc:
u = u_sequence[k]
else:
u = u_sequence[-1]
# State update
x = self.A @ x + self.B @ u
# Output computation
y_pred[k] = self.C @ x
return y_pred
def mpc_objective(self, u_flat: np.ndarray, x0: np.ndarray,
r: np.ndarray, u_prev: np.ndarray) -> float:
"""MPC objective function
Args:
u_flat: Flattened input sequence (Nc*2,)
x0: Current state
r: Setpoint sequence (Np x 2)
u_prev: Input at the previous step
Returns:
Value of the cost function
"""
# Reconstruct the input sequence
u_sequence = u_flat.reshape(self.Nc, 2)
# Prediction
y_pred = self.predict(x0, u_sequence)
# Tracking error
tracking_error = 0.0
for k in range(self.Np):
e = y_pred[k] - r[k]
tracking_error += e.T @ self.Q @ e
# Input change penalty
control_effort = 0.0
for k in range(self.Nc):
if k == 0:
du = u_sequence[k] - u_prev
else:
du = u_sequence[k] - u_sequence[k-1]
control_effort += du.T @ self.R @ du
return tracking_error + control_effort
def solve(self, x0: np.ndarray, r: np.ndarray,
u_prev: np.ndarray) -> np.ndarray:
"""Solve the MPC optimization problem
Args:
x0: Current state
r: Setpoint sequence (Np x 2)
u_prev: Input at the previous step
Returns:
First value of the optimal input sequence (2,)
"""
# Initial guess
u0_flat = np.zeros(self.Nc * 2)
# Bound constraints
bounds = []
for _ in range(self.Nc):
bounds.extend([
(self.u_min[0], self.u_max[0]),
(self.u_min[1], self.u_max[1])
])
# Rate-of-change constraint
def delta_u_constraint(u_flat):
u_seq = u_flat.reshape(self.Nc, 2)
violations = []
for k in range(self.Nc):
if k == 0:
du = np.abs(u_seq[k] - u_prev)
else:
du = np.abs(u_seq[k] - u_seq[k-1])
violations.extend((self.delta_u_max - du).tolist())
return np.array(violations)
from scipy.optimize import NonlinearConstraint
nlc = NonlinearConstraint(delta_u_constraint, 0, np.inf)
# Optimization
result = minimize(
lambda u: self.mpc_objective(u, x0, r, u_prev),
u0_flat,
method='SLSQP',
bounds=bounds,
constraints=[nlc],
options={'ftol': 1e-4, 'disp': False}
)
u_opt = result.x.reshape(self.Nc, 2)
return u_opt[0] # Execute only the first step
# ===================================
# Simulation run
# ===================================
if __name__ == "__main__":
mpc = DistillationMPC(prediction_horizon=10, control_horizon=5)
# Simulation settings
T_sim = 50 # Simulation time [min]
x = np.zeros(2) # Initial state (deviation from setpoint)
u = np.zeros(2) # Initial input
# Setpoint changes (step changes)
r_top = np.zeros(T_sim)
r_bottom = np.zeros(T_sim)
r_top[10:] = -1.5 # Lower the top temperature by 1.5 deg C after 10 minutes
r_bottom[30:] = 1.0 # Raise the bottom temperature by 1.0 deg C after 30 minutes
# For recording
x_history = [x.copy()]
u_history = [u.copy()]
print("="*70)
print("Distillation Column MPC Control Simulation")
print("="*70)
for k in range(T_sim):
# Setpoint sequence (over the prediction horizon)
r_horizon = np.zeros((mpc.Np, 2))
for i in range(mpc.Np):
if k + i < T_sim:
r_horizon[i] = [r_top[k+i], r_bottom[k+i]]
else:
r_horizon[i] = [r_top[-1], r_bottom[-1]]
# Compute optimal input with MPC
u = mpc.solve(x, r_horizon, u)
# Process update (with disturbance)
disturbance = np.random.randn(2) * 0.05
x = mpc.A @ x + mpc.B @ u + disturbance
# Recording
x_history.append(x.copy())
u_history.append(u.copy())
if k % 10 == 0:
print(f"Time {k:2d} min: top deviation={x[0]:+.2f} deg C, bottom deviation={x[1]:+.2f} deg C, "
f"reflux ratio={u[0]:+.2f}, reflux flow={u[1]:+.2f} kg/min")
# Convert results to arrays
x_history = np.array(x_history)
u_history = np.array(u_history)
print("\n" + "="*70)
print("Control Performance Evaluation")
print("="*70)
# Steady-state error
steady_state_error_top = np.abs(x_history[-10:, 0] - r_top[-1]).mean()
steady_state_error_bottom = np.abs(x_history[-10:, 1] - r_bottom[-1]).mean()
print(f"Top temperature steady-state error: {steady_state_error_top:.3f} deg C")
print(f"Bottom temperature steady-state error: {steady_state_error_bottom:.3f} deg C")
print(f"\n[OK] MPC achieved high-precision temperature control under constraints")
print(f"[OK] Good tracking performance across multiple setpoint changes")
Nonlinear processes require nonlinear MPC. We show a nonlinear MPC implementation for a reactor using CasADi (conceptual implementation).
"""
===================================
Example 4: Nonlinear MPC (CasADi)
===================================
Nonlinear model predictive control using CasADi.
Simultaneously controls the concentration and temperature of a
Continuous Stirred Tank Reactor (CSTR).
Nonlinear process model:
- Mass balance: dC/dt = (C_in - C)/tau - k*C
- Energy balance: dT/dt = (T_in - T)/tau + (-dH)*k*C/(rho*Cp) + Q/(V*rho*Cp)
CasADi requires a separate installation in practice; here we substitute a numerical implementation.
"""
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
from typing import Tuple
class NonlinearCSTRMPC:
"""MPC for a nonlinear CSTR reactor"""
def __init__(self, prediction_horizon: int = 20, control_horizon: int = 10):
self.Np = prediction_horizon
self.Nc = control_horizon
self.dt = 0.5 # Sampling time [min]
# Process parameters
self.V = 1.0 # Reactor volume [m3]
self.rho = 1000.0 # Density [kg/m3]
self.Cp = 4.18 # Specific heat [kJ/kg*K]
self.delta_H = -50000.0 # Heat of reaction [kJ/kmol]
# Arrhenius reaction rate
self.A = 1.0e10 # Pre-exponential factor [1/min]
self.Ea = 75000.0 # Activation energy [J/mol]
self.R = 8.314 # Gas constant [J/mol*K]
# Manipulated variable constraints
self.F_min, self.F_max = 0.02, 0.20 # Flow rate [m3/min]
self.Q_min, self.Q_max = -5000, 5000 # Heating/cooling [kJ/min]
# Inputs
self.C_in = 10.0 # Inlet concentration [kmol/m3]
self.T_in = 300.0 # Inlet temperature [K]
def reaction_rate(self, C: float, T: float) -> float:
"""Reaction rate constant (Arrhenius equation)"""
return self.A * np.exp(-self.Ea / (self.R * T)) * C
def cstr_model(self, state: np.ndarray, t: float,
F: float, Q: float) -> np.ndarray:
"""CSTR reactor differential equations
Args:
state: [concentration C, temperature T]
t: Time
F: Flow rate [m3/min]
Q: Heating rate [kJ/min]
Returns:
[dC/dt, dT/dt]
"""
C, T = state
tau = self.V / F # Residence time
# Reaction rate
r = self.reaction_rate(C, T)
# Mass balance
dC_dt = (self.C_in - C) / tau - r
# Energy balance
heat_reaction = (-self.delta_H) * r / (self.rho * self.Cp)
heat_exchange = Q / (self.V * self.rho * self.Cp)
dT_dt = (self.T_in - T) / tau + heat_reaction + heat_exchange
return np.array([dC_dt, dT_dt])
def simulate_step(self, state: np.ndarray,
F: float, Q: float) -> np.ndarray:
"""Simulate a single step
Args:
state: [C, T]
F: Flow rate
Q: Heating rate
Returns:
State at the next step
"""
t_span = [0, self.dt]
result = odeint(self.cstr_model, state, t_span, args=(F, Q))
return result[-1]
def predict_trajectory(self, state0: np.ndarray,
u_sequence: np.ndarray) -> np.ndarray:
"""Compute the predicted trajectory
Args:
state0: Initial state [C, T]
u_sequence: Input sequence (Nc x 2) [F, Q]
Returns:
Predicted state sequence (Np x 2)
"""
state = state0.copy()
trajectory = np.zeros((self.Np, 2))
for k in range(self.Np):
if k < self.Nc:
F, Q = u_sequence[k]
else:
F, Q = u_sequence[-1]
state = self.simulate_step(state, F, Q)
trajectory[k] = state
return trajectory
def mpc_objective(self, u_flat: np.ndarray, state0: np.ndarray,
setpoint: np.ndarray) -> float:
"""Nonlinear MPC objective function
Args:
u_flat: Flattened input sequence
state0: Current state
setpoint: Setpoint [C_sp, T_sp]
Returns:
Value of the cost function
"""
u_sequence = u_flat.reshape(self.Nc, 2)
# Predicted trajectory
trajectory = self.predict_trajectory(state0, u_sequence)
# Tracking error (weighted sum of squares)
Q = np.diag([10.0, 5.0]) # Weight for concentration and temperature
tracking_cost = 0.0
for k in range(self.Np):
error = trajectory[k] - setpoint
tracking_cost += error.T @ Q @ error
# Input change penalty
R = np.diag([100.0, 0.01]) # Weight for flow rate and heating rate
control_cost = 0.0
for k in range(self.Nc):
if k > 0:
du = u_sequence[k] - u_sequence[k-1]
control_cost += du.T @ R @ du
return tracking_cost + control_cost
def solve(self, state0: np.ndarray, setpoint: np.ndarray,
u_prev: np.ndarray) -> np.ndarray:
"""Nonlinear MPC optimization
Args:
state0: Current state
setpoint: Setpoint
u_prev: Input at the previous step
Returns:
Optimal input
"""
# Initial guess (hold the previous input)
u0_flat = np.tile(u_prev, self.Nc)
# Bound constraints
bounds = []
for _ in range(self.Nc):
bounds.append((self.F_min, self.F_max))
bounds.append((self.Q_min, self.Q_max))
# Optimization
result = minimize(
lambda u: self.mpc_objective(u, state0, setpoint),
u0_flat,
method='L-BFGS-B',
bounds=bounds,
options={'ftol': 1e-3, 'maxiter': 50}
)
u_opt = result.x.reshape(self.Nc, 2)
return u_opt[0]
# ===================================
# Simulation run
# ===================================
if __name__ == "__main__":
mpc = NonlinearCSTRMPC(prediction_horizon=20, control_horizon=10)
# Initial state
state = np.array([2.0, 350.0]) # [C=2 kmol/m3, T=350 K]
u = np.array([0.1, 0.0]) # [F=0.1 m3/min, Q=0 kJ/min]
# Setpoint
C_setpoint = 5.0 # kmol/m3
T_setpoint = 360.0 # K
setpoint = np.array([C_setpoint, T_setpoint])
# Simulation
T_sim = 100 # Number of steps
state_history = [state.copy()]
u_history = [u.copy()]
print("="*70)
print("MPC Control of a Nonlinear CSTR Reactor")
print("="*70)
print(f"Setpoint: C={C_setpoint} kmol/m3, T={T_setpoint} K\n")
for k in range(T_sim):
# Disturbance (inlet concentration variation)
if k == 40:
mpc.C_in = 12.0 # Concentration increases
print(f"Time {k*mpc.dt:.1f} min: disturbance occurred (inlet concentration 10 -> 12 kmol/m3)\n")
# Compute optimal input with MPC
u = mpc.solve(state, setpoint, u)
# Process update
state = mpc.simulate_step(state, u[0], u[1])
# Add noise
state += np.random.randn(2) * [0.05, 0.5]
# Recording
state_history.append(state.copy())
u_history.append(u.copy())
if k % 20 == 0:
print(f"Time {k*mpc.dt:5.1f} min: C={state[0]:.2f} kmol/m3, T={state[1]:.1f} K, "
f"F={u[0]:.3f} m3/min, Q={u[1]:+6.1f} kJ/min")
state_history = np.array(state_history)
u_history = np.array(u_history)
print("\n" + "="*70)
print("Control Performance")
print("="*70)
# Steady-state evaluation (last 20 steps)
C_error = np.abs(state_history[-20:, 0] - C_setpoint).mean()
T_error = np.abs(state_history[-20:, 1] - T_setpoint).mean()
print(f"Concentration error (steady-state): {C_error:.3f} kmol/m3")
print(f"Temperature error (steady-state): {T_error:.3f} K")
print(f"\n[OK] High-precision control achieved even for a nonlinear process")
print(f"[OK] Rapid disturbance rejection")
Reinforcement learning learns an optimal control policy through trial and error. We implement Deep Q-Network (DQN)-based trajectory optimization for a batch reactor.
"""
===================================
Example 5: Batch Process Control with DQN
===================================
Temperature trajectory optimization for a batch reactor using a Deep Q-Network (DQN).
Objective: maximize product purity while minimizing batch time.
Reinforcement learning formulation:
- State: [concentration A, concentration B, temperature, elapsed time]
- Action: temperature increase/decrease (discretized)
- Reward: product yield - time penalty
"""
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from collections import deque
import random
from typing import Tuple, List
class BatchReactorEnv:
"""Batch reactor environment"""
def __init__(self):
self.dt = 1.0 # Time step [min]
self.max_time = 120.0 # Maximum batch time [min]
# Reaction rate parameters
self.A1 = 1.0e8 # A->B (desired reaction)
self.E1 = 70000.0
self.A2 = 5.0e7 # B->C (side reaction)
self.E2 = 65000.0
self.R = 8.314
# Temperature constraint
self.T_min, self.T_max = 320.0, 380.0 # [K]
self.reset()
def reset(self) -> np.ndarray:
"""Reset the environment"""
self.C_A = 10.0 # Initial concentration A [mol/L]
self.C_B = 0.0 # Initial concentration B [mol/L]
self.C_C = 0.0 # Initial concentration C [mol/L]
self.T = 340.0 # Initial temperature [K]
self.time = 0.0
return self._get_state()
def _get_state(self) -> np.ndarray:
"""Get the state vector"""
return np.array([
self.C_A / 10.0, # Normalized
self.C_B / 10.0,
(self.T - 350.0) / 30.0,
self.time / self.max_time
], dtype=np.float32)
def step(self, action: int) -> Tuple[np.ndarray, float, bool]:
"""Execute a single step
Args:
action: 0=cool, 1=hold, 2=heat
Returns:
(next_state, reward, done)
"""
# Temperature change
delta_T = [-2.0, 0.0, 2.0][action]
self.T = np.clip(self.T + delta_T, self.T_min, self.T_max)
# Reaction rate constants
k1 = self.A1 * np.exp(-self.E1 / (self.R * self.T))
k2 = self.A2 * np.exp(-self.E2 / (self.R * self.T))
# Concentration update (first-order reaction)
dC_A = -k1 * self.C_A * self.dt
dC_B = (k1 * self.C_A - k2 * self.C_B) * self.dt
dC_C = k2 * self.C_B * self.dt
self.C_A += dC_A
self.C_B += dC_B
self.C_C += dC_C
self.time += self.dt
# Reward computation
# Maximize product B concentration (target: 8 mol/L or more)
product_reward = self.C_B
# Penalty for byproduct C
byproduct_penalty = -0.5 * self.C_C
# Time penalty (finish as fast as possible)
time_penalty = -0.01 * self.time
reward = product_reward + byproduct_penalty + time_penalty
# Termination condition
done = (self.time >= self.max_time) or (self.C_A < 0.5)
# Bonus: when the target is achieved
if done and self.C_B >= 7.5:
reward += 50.0
return self._get_state(), reward, done
class DQN(nn.Module):
"""Deep Q-Network"""
def __init__(self, state_dim: int, action_dim: int):
super(DQN, self).__init__()
self.fc = nn.Sequential(
nn.Linear(state_dim, 64),
nn.ReLU(),
nn.Linear(64, 64),
nn.ReLU(),
nn.Linear(64, action_dim)
)
def forward(self, x):
return self.fc(x)
class DQNAgent:
"""DQN agent"""
def __init__(self, state_dim: int, action_dim: int):
self.state_dim = state_dim
self.action_dim = action_dim
self.q_network = DQN(state_dim, action_dim)
self.target_network = DQN(state_dim, action_dim)
self.target_network.load_state_dict(self.q_network.state_dict())
self.optimizer = optim.Adam(self.q_network.parameters(), lr=1e-3)
self.memory = deque(maxlen=10000)
self.batch_size = 64
self.gamma = 0.99
self.epsilon = 1.0
self.epsilon_min = 0.05
self.epsilon_decay = 0.995
def select_action(self, state: np.ndarray) -> int:
"""epsilon-greedy action selection"""
if random.random() < self.epsilon:
return random.randint(0, self.action_dim - 1)
with torch.no_grad():
state_t = torch.FloatTensor(state).unsqueeze(0)
q_values = self.q_network(state_t)
return q_values.argmax().item()
def store_transition(self, state, action, reward, next_state, done):
"""Store an experience"""
self.memory.append((state, action, reward, next_state, done))
def train(self):
"""Train the network"""
if len(self.memory) < self.batch_size:
return
# Batch sampling
batch = random.sample(self.memory, self.batch_size)
states, actions, rewards, next_states, dones = zip(*batch)
states = torch.FloatTensor(np.array(states))
actions = torch.LongTensor(actions)
rewards = torch.FloatTensor(rewards)
next_states = torch.FloatTensor(np.array(next_states))
dones = torch.FloatTensor(dones)
# Current Q values
current_q = self.q_network(states).gather(1, actions.unsqueeze(1))
# Target Q values
with torch.no_grad():
next_q = self.target_network(next_states).max(1)[0]
target_q = rewards + self.gamma * next_q * (1 - dones)
# Loss computation and update
loss = nn.MSELoss()(current_q.squeeze(), target_q)
self.optimizer.zero_grad()
loss.backward()
self.optimizer.step()
# Epsilon decay
self.epsilon = max(self.epsilon_min, self.epsilon * self.epsilon_decay)
def update_target_network(self):
"""Update the target network"""
self.target_network.load_state_dict(self.q_network.state_dict())
# ===================================
# Training run
# ===================================
if __name__ == "__main__":
env = BatchReactorEnv()
agent = DQNAgent(state_dim=4, action_dim=3)
num_episodes = 200
rewards_history = []
print("="*70)
print("Training a Batch Reactor Controller with DQN")
print("="*70)
for episode in range(num_episodes):
state = env.reset()
episode_reward = 0
done = False
while not done:
action = agent.select_action(state)
next_state, reward, done = env.step(action)
agent.store_transition(state, action, reward, next_state, done)
agent.train()
state = next_state
episode_reward += reward
rewards_history.append(episode_reward)
# Update the target network
if episode % 10 == 0:
agent.update_target_network()
if episode % 20 == 0:
avg_reward = np.mean(rewards_history[-20:])
print(f"Episode {episode:3d}: avg reward={avg_reward:6.2f}, "
f"epsilon={agent.epsilon:.3f}, final C_B={env.C_B:.2f} mol/L")
print("\n" + "="*70)
print("Training complete - Test run")
print("="*70)
# Test run (greedy policy)
agent.epsilon = 0.0
state = env.reset()
done = False
trajectory = []
while not done:
action = agent.select_action(state)
trajectory.append({
'time': env.time,
'C_A': env.C_A,
'C_B': env.C_B,
'C_C': env.C_C,
'T': env.T,
'action': ['cool', 'hold', 'heat'][action]
})
next_state, reward, done = env.step(action)
state = next_state
print(f"\nBatch completion time: {env.time:.1f} min")
print(f"Final product concentration C_B: {env.C_B:.2f} mol/L")
print(f"Byproduct C_C: {env.C_C:.2f} mol/L")
print(f"Selectivity: {env.C_B / (env.C_B + env.C_C):.1%}")
print("\nOptimal temperature trajectory (first 10 steps):")
for i in range(min(10, len(trajectory))):
t = trajectory[i]
print(f" {t['time']:5.1f} min: T={t['T']:.1f}K, C_B={t['C_B']:.2f}, action={t['action']}")
print("\n[OK] DQN learned an optimal trajectory that balances batch time and product purity")
Proximal Policy Optimization (PPO) is effective for processes with continuous action spaces. We show an example of applying it to CSTR control.
"""
===================================
Example 6: Continuous Process Control with PPO
===================================
Control of a Continuous Stirred Tank Reactor using Proximal Policy Optimization (PPO).
An Actor-Critic architecture is used to handle the continuous action space (temperature, flow rate).
Objective: maintain the product concentration at its target value while minimizing energy consumption
"""
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from torch.distributions import Normal
from typing import Tuple
class CSTREnv:
"""Continuous stirred tank reactor environment"""
def __init__(self):
self.dt = 0.5 # Sampling time [min]
self.V = 1.0 # Reactor volume [m3]
# Reaction rate parameters
self.A = 1.0e9
self.Ea = 72000.0
self.R = 8.314
# Setpoint
self.C_target = 3.0 # Target product concentration [mol/L]
self.reset()
def reset(self) -> np.ndarray:
"""Reset the environment"""
self.C = 2.0 + np.random.randn() * 0.5 # Concentration [mol/L]
self.T = 350.0 + np.random.randn() * 5.0 # Temperature [K]
self.C_in = 8.0 # Inlet concentration [mol/L]
return self._get_state()
def _get_state(self) -> np.ndarray:
"""State vector"""
return np.array([
(self.C - self.C_target) / self.C_target, # Concentration deviation (normalized)
(self.T - 350.0) / 30.0 # Temperature deviation (normalized)
], dtype=np.float32)
def step(self, action: np.ndarray) -> Tuple[np.ndarray, float, bool]:
"""Execute a single step
Args:
action: [flow rate change rate, temperature change] (normalized to -1..1)
Returns:
(next_state, reward, done)
"""
# Convert the action to actual physical quantities
F = 0.1 + 0.05 * action[0] # Flow rate 0.05-0.15 [m3/min]
delta_T = 5.0 * action[1] # Temperature change -5..+5 [K]
self.T = np.clip(self.T + delta_T, 320.0, 380.0)
# Reaction rate constant
k = self.A * np.exp(-self.Ea / (self.R * self.T))
# Concentration update (CSTR material balance)
tau = self.V / F
dC = ((self.C_in - self.C) / tau - k * self.C) * self.dt
self.C += dC
# Noise (disturbance)
self.C += np.random.randn() * 0.05
self.T += np.random.randn() * 1.0
# Reward computation
# 1. Concentration tracking (primary objective)
error = abs(self.C - self.C_target)
tracking_reward = -10.0 * error
# 2. Energy penalty (heating cost)
heating_cost = -0.01 * abs(delta_T)
# 3. Flow rate change penalty (smooth operation)
flow_penalty = -0.1 * abs(action[0])
reward = tracking_reward + heating_cost + flow_penalty
# Bonus: within the target range (+/-0.2 mol/L)
if error < 0.2:
reward += 5.0
done = False # No termination for a continuous process
return self._get_state(), reward, done
class ActorCritic(nn.Module):
"""Actor-Critic network"""
def __init__(self, state_dim: int, action_dim: int):
super(ActorCritic, self).__init__()
# Shared layers
self.shared = nn.Sequential(
nn.Linear(state_dim, 64),
nn.Tanh(),
nn.Linear(64, 64),
nn.Tanh()
)
# Actor (policy)
self.actor_mean = nn.Linear(64, action_dim)
self.actor_log_std = nn.Parameter(torch.zeros(action_dim))
# Critic (value function)
self.critic = nn.Linear(64, 1)
def forward(self, state):
"""Forward pass"""
shared_features = self.shared(state)
# Actor: parameters of the normal distribution
action_mean = torch.tanh(self.actor_mean(shared_features))
action_std = torch.exp(self.actor_log_std)
# Critic: state value
value = self.critic(shared_features)
return action_mean, action_std, value
def get_action(self, state):
"""Sample an action"""
action_mean, action_std, value = self.forward(state)
dist = Normal(action_mean, action_std)
action = dist.sample()
log_prob = dist.log_prob(action).sum(-1)
return action, log_prob, value
class PPOAgent:
"""PPO agent"""
def __init__(self, state_dim: int, action_dim: int):
self.ac = ActorCritic(state_dim, action_dim)
self.optimizer = optim.Adam(self.ac.parameters(), lr=3e-4)
self.gamma = 0.99
self.lam = 0.95 # GAE lambda
self.clip_epsilon = 0.2
self.epochs = 10
def compute_gae(self, rewards, values, dones):
"""Generalized Advantage Estimation"""
advantages = []
gae = 0
for t in reversed(range(len(rewards))):
if t == len(rewards) - 1:
next_value = 0
else:
next_value = values[t + 1]
delta = rewards[t] + self.gamma * next_value - values[t]
gae = delta + self.gamma * self.lam * gae
advantages.insert(0, gae)
return advantages
def update(self, states, actions, old_log_probs, returns, advantages):
"""PPO update"""
states = torch.FloatTensor(np.array(states))
actions = torch.FloatTensor(np.array(actions))
old_log_probs = torch.FloatTensor(old_log_probs)
returns = torch.FloatTensor(returns)
advantages = torch.FloatTensor(advantages)
# Normalization
advantages = (advantages - advantages.mean()) / (advantages.std() + 1e-8)
for _ in range(self.epochs):
action_mean, action_std, values = self.ac(states)
dist = Normal(action_mean, action_std)
new_log_probs = dist.log_prob(actions).sum(-1)
# Importance sampling ratio
ratio = torch.exp(new_log_probs - old_log_probs)
# Clipped surrogate objective
surr1 = ratio * advantages
surr2 = torch.clamp(ratio, 1 - self.clip_epsilon, 1 + self.clip_epsilon) * advantages
actor_loss = -torch.min(surr1, surr2).mean()
# Value loss
critic_loss = nn.MSELoss()(values.squeeze(), returns)
# Total loss
loss = actor_loss + 0.5 * critic_loss
self.optimizer.zero_grad()
loss.backward()
self.optimizer.step()
# ===================================
# Training run
# ===================================
if __name__ == "__main__":
env = CSTREnv()
agent = PPOAgent(state_dim=2, action_dim=2)
num_episodes = 300
steps_per_episode = 200
print("="*70)
print("Training Continuous CSTR Control with PPO")
print("="*70)
for episode in range(num_episodes):
state = env.reset()
states, actions, log_probs, rewards, values = [], [], [], [], []
for step in range(steps_per_episode):
state_t = torch.FloatTensor(state)
action, log_prob, value = agent.ac.get_action(state_t.unsqueeze(0))
action_np = action.squeeze().detach().numpy()
next_state, reward, done = env.step(action_np)
states.append(state)
actions.append(action_np)
log_probs.append(log_prob.item())
rewards.append(reward)
values.append(value.item())
state = next_state
# GAE computation
advantages = agent.compute_gae(rewards, values, [False] * len(rewards))
returns = [adv + val for adv, val in zip(advantages, values)]
# PPO update
agent.update(states, actions, log_probs, returns, advantages)
if episode % 30 == 0:
avg_reward = np.mean(rewards)
final_error = abs(env.C - env.C_target)
print(f"Episode {episode:3d}: avg reward={avg_reward:6.2f}, "
f"final error={final_error:.3f} mol/L")
print("\n" + "="*70)
print("Training complete - Test run")
print("="*70)
# Test run
state = env.reset()
trajectory = []
for step in range(100):
state_t = torch.FloatTensor(state)
with torch.no_grad():
action_mean, _, _ = agent.ac(state_t.unsqueeze(0))
action = action_mean.squeeze().numpy()
trajectory.append({
'step': step,
'C': env.C,
'T': env.T,
'action_flow': action[0],
'action_temp': action[1]
})
next_state, reward, done = env.step(action)
state = next_state
# Performance evaluation
concentrations = [t['C'] for t in trajectory]
mean_error = np.mean([abs(c - env.C_target) for c in concentrations])
std_error = np.std([abs(c - env.C_target) for c in concentrations])
print(f"\nMean tracking error: {mean_error:.3f} mol/L")
print(f"Error standard deviation: {std_error:.3f} mol/L")
print("\nControl trajectory (first 10 steps):")
for i in range(10):
t = trajectory[i]
print(f" Step {t['step']:3d}: C={t['C']:.2f} mol/L, T={t['T']:.1f}K, "
f"Flow action={t['action_flow']:+.2f}, Temp action={t['action_temp']:+.2f}")
print("\n[OK] PPO learned continuous-valued control and achieved setpoint tracking")
In chemical plants, multiple objectives such as yield and energy consumption are often in a trade-off relationship. We implement multi-objective optimization using NSGA-II.
"""
===================================
Example 7: Multi-Objective Optimization (NSGA-II)
===================================
Multi-objective optimization of a chemical process using NSGA-II
(Non-dominated Sorting Genetic Algorithm II).
Objectives:
1. Maximize product yield
2. Minimize energy consumption
These objectives are mutually conflicting; we seek the Pareto-optimal solution set.
pymoo requires a separate installation in practice; here we substitute a basic genetic algorithm implementation.
"""
import numpy as np
from typing import List, Tuple
import pandas as pd
class MultiObjectiveOptimizer:
"""Multi-objective optimizer (NSGA-II style)"""
def __init__(self, population_size: int = 50, generations: int = 100):
self.pop_size = population_size
self.n_gen = generations
# Range of decision variables
# x = [temperature [K], pressure [bar], flow rate [m3/h], catalyst amount [kg]]
self.bounds_lower = np.array([350.0, 10.0, 1.0, 50.0])
self.bounds_upper = np.array([420.0, 50.0, 5.0, 200.0])
def process_model(self, x: np.ndarray) -> Tuple[float, float]:
"""Process model: compute yield and energy
Args:
x: [temperature, pressure, flow rate, catalyst amount]
Returns:
(yield, energy consumption)
"""
T, P, F, Cat = x
# Yield model (simplified)
# Higher temperature, pressure, and catalyst amount improve yield, but with saturation
T_factor = 1.0 / (1.0 + np.exp(-(T - 380.0) / 10.0))
P_factor = np.log(P / 10.0) / np.log(5.0)
Cat_factor = np.sqrt(Cat / 50.0)
yield_fraction = 0.5 + 0.4 * T_factor + 0.2 * P_factor + 0.15 * Cat_factor
# Production rate proportional to flow rate
productivity = F * yield_fraction # [m3/h * fraction]
# Energy consumption model
# Heating energy (proportional to temperature)
heating_energy = F * 1000.0 * 4.18 * (T - 300.0) / 3600.0 # [kW]
# Compression energy (nonlinear in pressure)
compression_energy = 10.0 * F * (P / 10.0)**1.2 # [kW]
total_energy = heating_energy + compression_energy
return productivity, total_energy
def objectives(self, x: np.ndarray) -> np.ndarray:
"""Objective functions (unified into a minimization problem)
Returns:
[-yield (maximize -> minimize), energy consumption (minimize)]
"""
productivity, energy = self.process_model(x)
return np.array([-productivity, energy])
def dominates(self, obj1: np.ndarray, obj2: np.ndarray) -> bool:
"""Determine whether obj1 dominates obj2"""
return np.all(obj1 <= obj2) and np.any(obj1 < obj2)
def non_dominated_sort(self, objectives: np.ndarray) -> List[List[int]]:
"""Non-dominated sorting
Args:
objectives: (pop_size x 2) objective function values
Returns:
List of indices for each front
"""
pop_size = len(objectives)
domination_count = np.zeros(pop_size, dtype=int)
dominated_solutions = [[] for _ in range(pop_size)]
fronts = [[]]
for i in range(pop_size):
for j in range(i + 1, pop_size):
if self.dominates(objectives[i], objectives[j]):
dominated_solutions[i].append(j)
domination_count[j] += 1
elif self.dominates(objectives[j], objectives[i]):
dominated_solutions[j].append(i)
domination_count[i] += 1
for i in range(pop_size):
if domination_count[i] == 0:
fronts[0].append(i)
current_front = 0
while len(fronts[current_front]) > 0:
next_front = []
for i in fronts[current_front]:
for j in dominated_solutions[i]:
domination_count[j] -= 1
if domination_count[j] == 0:
next_front.append(j)
current_front += 1
fronts.append(next_front)
return fronts[:-1] # Exclude the trailing empty list
def crowding_distance(self, objectives: np.ndarray, front: List[int]) -> np.ndarray:
"""Compute the crowding distance"""
n = len(front)
if n <= 2:
return np.full(n, np.inf)
distances = np.zeros(n)
for m in range(objectives.shape[1]): # For each objective
sorted_idx = np.argsort(objectives[front, m])
distances[sorted_idx[0]] = np.inf
distances[sorted_idx[-1]] = np.inf
obj_range = objectives[front[sorted_idx[-1]], m] - objectives[front[sorted_idx[0]], m]
if obj_range == 0:
continue
for i in range(1, n - 1):
distances[sorted_idx[i]] += (
(objectives[front[sorted_idx[i + 1]], m] -
objectives[front[sorted_idx[i - 1]], m]) / obj_range
)
return distances
def optimize(self) -> Tuple[np.ndarray, np.ndarray]:
"""Run NSGA-II optimization
Returns:
(Pareto solution set, objective function values)
"""
# Generate the initial population
population = np.random.uniform(
self.bounds_lower,
self.bounds_upper,
(self.pop_size, len(self.bounds_lower))
)
for generation in range(self.n_gen):
# Evaluate the objective functions
objectives = np.array([self.objectives(ind) for ind in population])
# Non-dominated sorting
fronts = self.non_dominated_sort(objectives)
# Selection for the next generation
next_population = []
for front in fronts:
if len(next_population) + len(front) <= self.pop_size:
next_population.extend(front)
else:
# Sort by crowding distance
distances = self.crowding_distance(objectives, front)
sorted_idx = np.argsort(distances)[::-1]
remaining = self.pop_size - len(next_population)
next_population.extend([front[i] for i in sorted_idx[:remaining]])
break
# Crossover and mutation
selected = population[next_population]
offspring = []
for i in range(0, len(selected) - 1, 2):
# SBX crossover (simplified)
alpha = np.random.rand(len(self.bounds_lower))
child1 = alpha * selected[i] + (1 - alpha) * selected[i + 1]
child2 = (1 - alpha) * selected[i] + alpha * selected[i + 1]
# Polynomial mutation (simplified)
if np.random.rand() < 0.1:
child1 += np.random.randn(len(self.bounds_lower)) * 0.1 * (self.bounds_upper - self.bounds_lower)
if np.random.rand() < 0.1:
child2 += np.random.randn(len(self.bounds_lower)) * 0.1 * (self.bounds_upper - self.bounds_lower)
# Bound constraints
child1 = np.clip(child1, self.bounds_lower, self.bounds_upper)
child2 = np.clip(child2, self.bounds_lower, self.bounds_upper)
offspring.extend([child1, child2])
population = np.array(offspring[:self.pop_size])
if generation % 20 == 0:
print(f"Generation {generation}: Pareto front size={len(fronts[0])}")
# Final evaluation
objectives = np.array([self.objectives(ind) for ind in population])
fronts = self.non_dominated_sort(objectives)
pareto_front = fronts[0]
return population[pareto_front], objectives[pareto_front]
# ===================================
# Example run
# ===================================
if __name__ == "__main__":
optimizer = MultiObjectiveOptimizer(population_size=50, generations=100)
print("="*70)
print("Multi-Objective Optimization (NSGA-II)")
print("="*70)
print("Objective 1: Maximize product yield")
print("Objective 2: Minimize energy consumption\n")
# Run optimization
pareto_solutions, pareto_objectives = optimizer.optimize()
print("\n" + "="*70)
print(f"Pareto-optimal solutions: {len(pareto_solutions)}")
print("="*70)
# Organize the results
results = []
for i, (sol, obj) in enumerate(zip(pareto_solutions, pareto_objectives)):
T, P, F, Cat = sol
productivity = -obj[0] # Undo the minimization sign flip
energy = obj[1]
results.append({
'Solution #': i + 1,
'Temperature [K]': T,
'Pressure [bar]': P,
'Flow rate [m3/h]': F,
'Catalyst [kg]': Cat,
'Yield productivity [m3/h]': productivity,
'Energy [kW]': energy,
'Energy intensity [kW/(m3/h)]': energy / productivity
})
df = pd.DataFrame(results)
# Display representative solutions (yield-focused, balanced, energy-focused)
print("\n[Representative Pareto Solutions]")
print("\n1. Yield-focused (high energy consumption):")
idx_max_yield = df['Yield productivity [m3/h]'].idxmax()
print(df.loc[idx_max_yield].to_string())
print("\n2. Balanced:")
df['Balance score'] = (df['Yield productivity [m3/h]'].rank() + (1 / df['Energy [kW]']).rank()) / 2
idx_balanced = df['Balance score'].idxmax()
print(df.loc[idx_balanced].to_string())
print("\n3. Energy-focused (lower yield):")
idx_min_energy = df['Energy [kW]'].idxmin()
print(df.loc[idx_min_energy].to_string())
print("\n" + "="*70)
print("Trade-off Analysis")
print("="*70)
# Energy cost increase for a 10% yield improvement
sorted_df = df.sort_values('Yield productivity [m3/h]')
if len(sorted_df) > 1:
yield_range = sorted_df['Yield productivity [m3/h]'].max() - sorted_df['Yield productivity [m3/h]'].min()
energy_range = sorted_df['Energy [kW]'].max() - sorted_df['Energy [kW]'].min()
print(f"Energy increase for a 10% yield improvement: approx. {energy_range / yield_range * 0.1 * sorted_df['Yield productivity [m3/h]'].mean():.1f} kW")
print("\n[OK] The Pareto front lets decision-makers choose their preferred yield-energy trade-off")
Finally, we show a practical implementation example of a plant control system that integrates the RTO and APC layers.
"""
===================================
Example 8: Integrated APC + Optimization System
===================================
A hierarchical plant control system that integrates the Real-Time Optimization (RTO)
layer and the Advanced Process Control (APC) layer.
Hierarchical structure:
- RTO layer (upper): determines optimal operating conditions through economic optimization
- APC layer (lower): uses MPC to track the RTO setpoints with high precision
Application example: a chemical plant (distillation column + reactor)
"""
import numpy as np
from scipy.optimize import minimize
from typing import Dict, Tuple
import pandas as pd
class RTOLayer:
"""Real-time optimization layer"""
def __init__(self):
# Economic parameters
self.product_price = 120.0 # $/ton
self.feed_cost = 50.0 # $/ton
self.steam_cost = 25.0 # $/ton
self.power_cost = 0.10 # $/kWh
def steady_state_model(self, x: np.ndarray) -> Dict:
"""Steady-state process model
Args:
x: [reaction temperature, distillation reflux ratio]
Returns:
Dictionary of process outputs
"""
T_reactor, reflux_ratio = x
# Reaction yield (temperature-dependent)
yield_base = 0.75
temp_effect = 0.002 * (T_reactor - 370.0)
yield_fraction = yield_base + temp_effect
# Product purity (reflux ratio-dependent)
purity_base = 0.90
reflux_effect = 0.15 * (1.0 - np.exp(-0.5 * (reflux_ratio - 2.0)))
purity = min(0.99, purity_base + reflux_effect)
# Utility consumption
reactor_heat = 50.0 + 0.5 * (T_reactor - 350.0)**2 # kW
steam_consumption = 2.0 + 0.8 * reflux_ratio # ton/h
power = 30.0 + 5.0 * reflux_ratio # kW
return {
'yield': yield_fraction,
'purity': purity,
'reactor_heat': reactor_heat,
'steam': steam_consumption,
'power': power
}
def economic_objective(self, x: np.ndarray, feed_rate: float) -> float:
"""Economic objective function: negative profit
Args:
x: [reaction temperature, reflux ratio]
feed_rate: Feed flow rate [ton/h]
Returns:
-profit [$/h]
"""
outputs = self.steady_state_model(x)
# Product production rate
product_rate = feed_rate * outputs['yield'] * outputs['purity']
# Revenue
revenue = product_rate * self.product_price
# Cost
feed_cost = feed_rate * self.feed_cost
steam_cost = outputs['steam'] * self.steam_cost
power_cost = outputs['power'] * self.power_cost
profit = revenue - feed_cost - steam_cost - power_cost
return -profit # Negated for minimization
def optimize(self, feed_rate: float) -> Dict:
"""Run RTO
Args:
feed_rate: Current feed flow rate [ton/h]
Returns:
Optimal operating conditions
"""
# Initial guess
x0 = np.array([370.0, 3.0])
# Bound constraints
bounds = [
(350.0, 390.0), # Reaction temperature [K]
(2.0, 5.0) # Reflux ratio
]
# Constraint: product purity of at least 95%
def purity_constraint(x):
outputs = self.steady_state_model(x)
return outputs['purity'] - 0.95
from scipy.optimize import NonlinearConstraint
nlc = NonlinearConstraint(purity_constraint, 0, np.inf)
# Optimization
result = minimize(
lambda x: self.economic_objective(x, feed_rate),
x0,
method='SLSQP',
bounds=bounds,
constraints=[nlc]
)
T_opt, reflux_opt = result.x
outputs = self.steady_state_model(result.x)
return {
'T_reactor_sp': T_opt,
'reflux_ratio_sp': reflux_opt,
'predicted_yield': outputs['yield'],
'predicted_purity': outputs['purity'],
'predicted_profit': -result.fun
}
class APCLayer:
"""Advanced control layer (MPC)"""
def __init__(self):
# Process model (linearized approximation)
# State: [reaction temperature deviation, reflux ratio deviation]
self.A = np.array([
[0.90, 0.05],
[0.00, 0.85]
])
self.B = np.array([
[0.8, 0.0],
[0.0, 0.6]
])
# Control horizon
self.Np = 15
self.Nc = 8
def mpc_control(self, current_state: np.ndarray,
setpoint: np.ndarray) -> np.ndarray:
"""MPC control
Args:
current_state: Current state deviation [T deviation, reflux ratio deviation]
setpoint: Setpoint deviation (command from RTO)
Returns:
Optimal manipulated variable
"""
# Simplified MPC (a more detailed implementation would be needed in practice)
# Here we substitute proportional control
Kp = np.array([2.0, 1.5])
u = Kp * (setpoint - current_state)
# Manipulated variable constraints
u = np.clip(u, [-5.0, -0.5], [5.0, 0.5])
return u
class IntegratedControlSystem:
"""Integrated control system"""
def __init__(self):
self.rto = RTOLayer()
self.apc = APCLayer()
# RTO execution interval (longer than APC)
self.rto_interval = 30 # 30x the APC cycle
# Current state
self.T_reactor = 365.0
self.reflux_ratio = 3.2
def run_rto(self, feed_rate: float) -> Dict:
"""Execute the RTO layer"""
print("\n" + "="*70)
print("RTO Layer: Running economic optimization")
print("="*70)
result = self.rto.optimize(feed_rate)
print(f" Optimal reaction temperature: {result['T_reactor_sp']:.1f} K")
print(f" Optimal reflux ratio: {result['reflux_ratio_sp']:.2f}")
print(f" Predicted yield: {result['predicted_yield']:.1%}")
print(f" Predicted purity: {result['predicted_purity']:.1%}")
print(f" Predicted profit: ${result['predicted_profit']:.2f}/h")
return result
def run_apc(self, rto_setpoint: Dict) -> np.ndarray:
"""Execute the APC layer"""
# Current deviation
current_state = np.array([
self.T_reactor - rto_setpoint['T_reactor_sp'],
self.reflux_ratio - rto_setpoint['reflux_ratio_sp']
])
# Target deviation (zero)
setpoint = np.array([0.0, 0.0])
# MPC control
u = self.apc.mpc_control(current_state, setpoint)
return u
def simulate(self, feed_rate: float, simulation_steps: int = 100):
"""Simulate the integrated system
Args:
feed_rate: Feed flow rate [ton/h]
simulation_steps: Number of simulation steps
"""
print("="*70)
print("Integrated APC + Optimization System Simulation")
print("="*70)
# Initial RTO run
rto_result = self.run_rto(feed_rate)
history = []
for step in range(simulation_steps):
# RTO update (periodic)
if step % self.rto_interval == 0 and step > 0:
rto_result = self.run_rto(feed_rate)
# APC execution (every step)
u = self.run_apc(rto_result)
# Process update (simplified)
self.T_reactor += u[0] * 0.5 + np.random.randn() * 0.5
self.reflux_ratio += u[1] * 0.3 + np.random.randn() * 0.05
# Physical constraints
self.T_reactor = np.clip(self.T_reactor, 350.0, 390.0)
self.reflux_ratio = np.clip(self.reflux_ratio, 2.0, 5.0)
# Recording
outputs = self.rto.steady_state_model(
[self.T_reactor, self.reflux_ratio]
)
history.append({
'step': step,
'T_reactor': self.T_reactor,
'T_setpoint': rto_result['T_reactor_sp'],
'reflux': self.reflux_ratio,
'reflux_setpoint': rto_result['reflux_ratio_sp'],
'yield': outputs['yield'],
'purity': outputs['purity']
})
if step % 20 == 0:
print(f"\nStep {step:3d}:")
print(f" Reaction temperature: {self.T_reactor:.1f}K (SP: {rto_result['T_reactor_sp']:.1f}K)")
print(f" Reflux ratio: {self.reflux_ratio:.2f} (SP: {rto_result['reflux_ratio_sp']:.2f})")
print(f" Yield: {outputs['yield']:.1%}, Purity: {outputs['purity']:.1%}")
return pd.DataFrame(history)
# ===================================
# Example run
# ===================================
if __name__ == "__main__":
system = IntegratedControlSystem()
# Feed flow rate
feed_rate = 10.0 # ton/h
# Run simulation
df = system.simulate(feed_rate, simulation_steps=100)
print("\n" + "="*70)
print("Control Performance Evaluation")
print("="*70)
# Tracking performance
T_error = np.abs(df['T_reactor'] - df['T_setpoint']).mean()
reflux_error = np.abs(df['reflux'] - df['reflux_setpoint']).mean()
print(f"\nMean reaction temperature tracking error: {T_error:.2f} K")
print(f"Mean reflux ratio tracking error: {reflux_error:.3f}")
# Process performance
avg_yield = df['yield'].mean()
avg_purity = df['purity'].mean()
print(f"\nAverage yield: {avg_yield:.1%}")
print(f"Average purity: {avg_purity:.1%}")
print("\n" + "="*70)
print("System Highlights")
print("="*70)
print("[OK] The RTO layer determines economically optimal operating conditions")
print("[OK] The APC layer (MPC) achieves high-precision setpoint tracking")
print("[OK] Hierarchical structure distributes computational load and enables real-time control")
print("[OK] Maximizes overall plant economics while guaranteeing quality")
In this chapter, we learned implementation techniques for real-time optimization and Advanced Process Control in chemical plants. The key points are as follows:
Examples of implementation effects:
In the next chapter, we will learn about supply chain optimization and digital twins for the entire plant. We will cover demand forecasting, production planning, inventory optimization, and real-time simulation for plant operations support at the implementation level.