Nonlinear Model Predictive Control

Dynamic control is also known as Nonlinear Model Predictive Control (NMPC) or simply as Nonlinear Control (NLC). NLC with predictive models is a dynamic optimization approach that seeks to follow a trajectory or drive certain values to maximum or minimum levels.


Objective: Design a controller to maintain temperature of a chemical reactor. Develop 3 separate controllers (PID, Linear MPC, Nonlinear MPC) in Python, MATLAB, or Simulink. Demonstrate controller performance with steps in the set point and disturbance changes. Estimated time: 3 hours.

A reactor is used to convert a hazardous chemical A to an acceptable chemical B in waste stream before entering a nearby lake. This particular reactor is dynamically modeled as a Continuously Stirred Tank Reactor (CSTR) with a simplified kinetic mechanism that describes the conversion of reactant A to product B with an irreversible and exothermic reaction. It is desired to maintain the temperature at a constant setpoint that maximizes the destruction of A (highest possible temperature). Adjust the jacket temperature (Tc) to maintain a desired reactor temperature and minimize the concentration of A. The reactor temperature should never exceed 400 K. The cooling jacket temperature can be adjusted between 250 K and 350 K.

Step Test (Doublet)

A step test is required to obtain a process model for the PID controller and the linear model predictive controller. It is a first step in developing a controller. The following code implements a doublet test. A doublet test starts with the system at steady state. Three moves of Manipulated Variable (MV) are made with sufficient time to nearly reach steady state conditions at two other operating points. The steps are above and below the nominal operating conditions. In this case, the cooling jacket temperature is raised, lowered, and brought back to 300 K (nominal operating condition.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# define CSTR model
def cstr(x,t,u,Tf,Caf):
    # Inputs (3):
    # Temperature of cooling jacket (K)
    Tc = u
    # Tf = Feed Temperature (K)
    # Caf = Feed Concentration (mol/m^3)

    # States (2):
    # Concentration of A in CSTR (mol/m^3)
    Ca = x[0]
    # Temperature in CSTR (K)
    T = x[1]

    # Parameters:
    # Volumetric Flowrate (m^3/sec)
    q = 100
    # Volume of CSTR (m^3)
    V = 100
    # Density of A-B Mixture (kg/m^3)
    rho = 1000
    # Heat capacity of A-B Mixture (J/kg-K)
    Cp = 0.239
    # Heat of reaction for A->B (J/mol)
    mdelH = 5e4
    # E - Activation energy in the Arrhenius Equation (J/mol)
    # R - Universal Gas Constant = 8.31451 J/mol-K
    EoverR = 8750
    # Pre-exponential factor (1/sec)
    k0 = 7.2e10
    # U - Overall Heat Transfer Coefficient (W/m^2-K)
    # A - Area - this value is specific for the U calculation (m^2)
    UA = 5e4
    # reaction rate
    rA = k0*np.exp(-EoverR/T)*Ca

    # Calculate concentration derivative
    dCadt = q/V*(Caf - Ca) - rA
    # Calculate temperature derivative
    dTdt = q/V*(Tf - T) \
            + mdelH/(rho*Cp)*rA \
            + UA/V/rho/Cp*(Tc-T)

    # Return xdot:
    xdot = np.zeros(2)
    xdot[0] = dCadt
    xdot[1] = dTdt
    return xdot

# Steady State Initial Conditions for the States
Ca_ss = 0.87725294608097
T_ss = 324.475443431599
x0 = np.empty(2)
x0[0] = Ca_ss
x0[1] = T_ss

# Steady State Initial Condition
u_ss = 300.0
# Feed Temperature (K)
Tf = 350
# Feed Concentration (mol/m^3)
Caf = 1

# Time Interval (min)
t = np.linspace(0,25,251)

# Store results for plotting
Ca = np.ones(len(t)) * Ca_ss
T = np.ones(len(t)) * T_ss
u = np.ones(len(t)) * u_ss

# Step cooling temperature to 295
u[10:100] = 303.0
u[100:190] = 297.0
u[190:] = 300.0

# Simulate CSTR
for i in range(len(t)-1):
    ts = [t[i],t[i+1]]
    y = odeint(cstr,x0,ts,args=(u[i+1],Tf,Caf))
    Ca[i+1] = y[-1][0]
    T[i+1] = y[-1][1]
    x0[0] = Ca[i+1]
    x0[1] = T[i+1]

# Construct results and save data file
# Column 1 = time
# Column 2 = cooling temperature
# Column 3 = reactor temperature
data = np.vstack((t,u,T)) # vertical stack
data = data.T             # transpose data

# Plot the results
plt.ylabel('Cooling T (K)')
plt.legend(['Jacket Temperature'],loc='best')

plt.ylabel('Ca (mol/L)')
plt.legend(['Reactor Concentration'],loc='best')

plt.ylabel('T (K)')
plt.xlabel('Time (min)')
plt.legend(['Reactor Temperature'],loc='best')

See additional information on this application on the Process Control Class Web-page.

Solution in Python

Solution in Simulink