Nonlinear Model Predictive Control
Main.NonlinearControl History
Hide minor edits - Show changes to markup
m = GEKKO(remote=True) # remote=True for MacOS
m = GEKKO(remote=True)
plt.plot(t,u,'b--',linewidth=3)
plt.plot(t,u,'b--',lw=3)
plt.plot(t,Ca,'r-',linewidth=3)
plt.plot(t,Ca,'r-',lw=3)
plt.plot(t,T,'k.-',linewidth=3)
plt.plot(t,T,'k.-',lw=3)
plt.plot(t,Tc,'b--',linewidth=3)
plt.plot(t,Tc,'b--',lw=3)
plt.plot(t,Ca,'r-',linewidth=3)
plt.plot(t,Ca,'r-',lw=3)
plt.plot(t,T,'k.-',linewidth=3)
plt.plot(t,T,'k.-',lw=3)
plt.plot(t[0:i],u[0:i],'b--',linewidth=3)
plt.plot(t[0:i],u[0:i],'b--',lw=3)
plt.plot(t[0:i],Ca[0:i],'r-',linewidth=3)
plt.plot(t[0:i],Ca[0:i],'r-',lw=3)
plt.plot(t[0:i],Tsp[0:i],'k-',linewidth=3,label=r'$T_{sp}$') plt.plot(t[0:i],T[0:i],'b.-',linewidth=3,label=r'$T_{meas}$')
plt.plot(t[0:i],Tsp[0:i],'k-',lw=3,label=r'$T_{sp}$') plt.plot(t[0:i],T[0:i],'b.-',lw=3,label=r'$T_{meas}$')
plt.plot(t[0:i],Tc[0:i],'b--',linewidth=3)
plt.plot(t[0:i],Tc[0:i],'b--',lw=3)
plt.plot(t[0:i],Ca[0:i],'r-',linewidth=3)
plt.plot(t[0:i],Ca[0:i],'r-',lw=3)
plt.plot(t[0:i],Tsp[0:i],'k-',linewidth=3,label=r'$T_{sp}$') plt.plot(t[0:i],T[0:i],'b.-',linewidth=3,label=r'$T_{meas}$')
plt.plot(t[0:i],Tsp[0:i],'k-',lw=3,label=r'$T_{sp}$') plt.plot(t[0:i],T[0:i],'b.-',lw=3,label=r'$T_{meas}$')
plt.plot(t[0:i],u[0:i],'b--',linewidth=3)
plt.plot(t[0:i],u[0:i],'b--',lw=3)
plt.plot(t[0:i],Ca[0:i],'b.-',linewidth=3,label=r'$C_A$') plt.plot([0,t[i-1]],[0.2,0.2],'r--',linewidth=2,label='limit')
plt.plot(t[0:i],Ca[0:i],'b.-',lw=3,label=r'$C_A$') plt.plot([0,t[i-1]],[0.2,0.2],'r--',lw=2,label='limit')
plt.plot(t[0:i],Tsp[0:i],'k-',linewidth=3,label=r'$T_{sp}$') plt.plot(t[0:i],T[0:i],'b.-',linewidth=3,label=r'$T_{meas}$') plt.plot([0,t[i-1]],[400,400],'r--',linewidth=2,label='limit')
plt.plot(t[0:i],Tsp[0:i],'k-',lw=3,label=r'$T_{sp}$') plt.plot(t[0:i],T[0:i],'b.-',lw=3,label=r'$T_{meas}$') plt.plot([0,t[i-1]],[400,400],'r--',lw=2,label='limit')
m.time = np.linspace(0,5,51)
m.time = np.linspace(0,2,21)
m.T.TAU = 0.5
m.T.TAU = 1.2
t = np.linspace(0,10,101)
t = np.linspace(0,20,201)
Tsp[0:20] = 330.0 Tsp[20:40] = 350.0 Tsp[40:60] = 370.0 Tsp[60:] = 390.0
Tsp[0:40] = 330.0 Tsp[40:80] = 350.0 Tsp[80:120] = 370.0 Tsp[120:] = 390.0
plt.legend(['Reactor Temperature'],loc='best')
plt.legend(['Temperature SP','Reactor Temperature'],loc='best')
(:toggle hide gekko_mpc button show="Show GEKKO Linear MPC Code":)
(:toggle hide gekko_mpc button show="Show GEKKO Linear First-Order MPC Code":)
plt.ylabel('Cooling T (K)') plt.legend(['Jacket Temperature'],loc='best') plt.subplot(3,1,2) plt.plot(t[0:i],Ca[0:i],'r-',linewidth=3) plt.ylabel('Ca (mol/L)') plt.legend(['Reactor Concentration'],loc='best') plt.subplot(3,1,3) plt.plot(t[0:i],Tsp[0:i],'k-',linewidth=3,label=r'$T_{sp}$') plt.plot(t[0:i],T[0:i],'b.-',linewidth=3,label=r'$T_{meas}$') plt.ylabel('T (K)') plt.xlabel('Time (min)') plt.legend(['Reactor Temperature'],loc='best') plt.draw() plt.pause(0.01)
(:sourceend:) (:divend:)
(:toggle hide gekko_arx_mpc button show="Show GEKKO ARX MPC Code":) (:div id=gekko_arx_mpc:)
Linear ARX MPC: Reactor Runaway and Upper Limit Temperature Violation (>400K)
(:source lang=python:) import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint import pandas as pd from gekko import GEKKO
- load data and parse into columns
url = 'http://apmonitor.com/do/uploads/Main/cstr_step_tests.txt' data = pd.read_csv(url) print(data.head())
- generate time-series model
t = data['Time'] u = data['Tc'] y = data['T'] m = GEKKO(remote=True)
- system identification
na = 2 # output coefficients nb = 2 # input coefficients yp,p,K = m.sysid(t,u,y,na,nb,shift='init',scale=True,objf=100,diaglevel=1)
- plot results of fitting
plt.figure() plt.subplot(2,1,1) plt.plot(t,u) plt.legend([r'$T_c$']) plt.ylabel('MV') plt.subplot(2,1,2) plt.plot(t,y) plt.plot(t,yp) plt.legend([r'$T_{meas}$',r'$T_{pred}$']) plt.ylabel('CV') plt.xlabel('Time') plt.savefig('sysid.png') plt.show()
- step test model
yc,uc = m.arx(p)
- rename MV and CV
m.Tc = uc[0] m.T = yc[0]
- steady state initialization
m.options.IMODE = 1 m.Tc.value = 280 m.solve(disp=True)
- GEKKO linear MPC
m.time = np.linspace(0,5,51)
- MV tuning
m.Tc.STATUS = 1 m.Tc.FSTATUS = 0 m.Tc.DMAX = 100 m.Tc.DCOST = 0.1 m.Tc.DMAXHI = 5 # constrain movement up m.Tc.DMAXLO = -100 # quick action down m.Tc.UPPER = 350 m.Tc.LOWER = 250
- CV tuning
m.T.STATUS = 1 m.T.FSTATUS = 1 m.T.SP = 330 m.T.TR_INIT = 1 m.T.TAU = 0.5 m.options.CV_TYPE = 2 m.options.IMODE = 6 m.options.SOLVER = 3
- define CSTR (plant)
def cstr(x,t,Tc):
Ca,T = x Tf = 350; Caf = 1.0; q = 100; V = 100 rho = 1000; Cp = 0.239; mdelH = 5e4 EoverR = 8750; k0 = 7.2e10; UA = 5e4 rA = k0*np.exp(-EoverR/T)*Ca dCadt = q/V*(Caf - Ca) - rA dTdt = q/V*(Tf - T) + mdelH/(rho*Cp)*rA + UA/V/rho/Cp*(Tc-T) return [dCadt,dTdt]
- Time Interval (min)
t = np.linspace(0,10,101)
- Store results for plotting
Ca_ss = 1; T_ss = 304; Tc_ss = 280 Ca = np.ones(len(t)) * Ca_ss T = np.ones(len(t)) * T_ss Tsp = np.ones(len(t)) * T_ss Tc = np.ones(len(t)) * Tc_ss
- Set point steps
Tsp[0:20] = 330.0 Tsp[20:40] = 350.0 Tsp[40:60] = 370.0 Tsp[60:] = 390.0
- Create plot
plt.figure(figsize=(10,7)) plt.ion() plt.show()
- Simulate CSTR
x0 = [Ca_ss,T_ss] for i in range(len(t)-1):
y = odeint(cstr,x0,[0,0.05],args=(Tc[i],)) # retrieve measurements Ca[i+1] = y[-1][0] T[i+1] = y[-1][1] # insert measurement m.T.MEAS = T[i+1] # update setpoint m.T.SP = Tsp[i+1] # solve MPC m.solve(disp=True) # retrieve new Tc value Tc[i+1] = m.Tc.NEWVAL # update initial conditions x0[0] = Ca[i+1] x0[1] = T[i+1] # Plot the results plt.clf() plt.subplot(3,1,1) plt.plot(t[0:i],Tc[0:i],'b--',linewidth=3)
(:toggle hide arx_fit button show="Show Python ARX System ID":) (:div id=arx_fit:)
(:source lang=python:) from gekko import GEKKO import pandas as pd import matplotlib.pyplot as plt import numpy as np
- load data and parse into columns
url = 'http://apmonitor.com/do/uploads/Main/cstr_step_tests.txt' data = pd.read_csv(url) print(data.head())
- generate time-series model
t = data['Time'] u = data['Tc'] y = data['T'] m = GEKKO(remote=True) # remote=True for MacOS
- system identification
na = 2 # output coefficients nb = 2 # input coefficients yp,p,K = m.sysid(t,u,y,na,nb,shift='init',scale=True,objf=100,diaglevel=1)
- plot results of fitting
plt.figure() plt.subplot(2,1,1) plt.plot(t,u) plt.legend([r'$T_c$']) plt.ylabel('MV') plt.subplot(2,1,2) plt.plot(t,y) plt.plot(t,yp) plt.legend([r'$T_{meas}$',r'$T_{pred}$']) plt.ylabel('CV') plt.xlabel('Time') plt.savefig('sysid.png')
- step test model
yc,uc = m.arx(p)
- rename MV and CV
Tc = uc[0] T = yc[0]
- steady state initialization
m.options.IMODE = 1 Tc.value = 300 m.solve(disp=False)
- dynamic simulation (step test validation)
m.time = np.linspace(0,2,21) m.options.IMODE = 4 Tc.value = np.ones(21)*300 Tc.value[5:] = 305 m.solve(disp=False)
plt.figure() plt.subplot(2,1,1) plt.title('Step Test') plt.plot(m.time,Tc.value,'b-',label='Cooling Jacket') plt.ylabel(r'$T_c (K)$') plt.legend() plt.subplot(2,1,2) plt.plot(m.time,T.value,'r-',label='Reactor') plt.ylabel('T (K)') plt.xlabel('Time (min)') plt.legend()
plt.show() (:sourceend:) (:divend:)
See additional information on this application on the Process Control Class Web-page.
See additional information on this application with CSTR PID Control with MATLAB and CSTR PID Control with Python.
Linear MPC: Reactor Runaway and Upper Limit Temperature Violation (>400K)
Nonlinear MPC: Acceptable Temperature (<400K) and Concentration (<0.2) Control
Model Identification
There are many methods to develop a controller model. For a PID controller, an FOPDT model is one method to obtain IMC tuning parameters. For linear MPC, there are many options to obtain a controller model through identification methods. For nonlinear MPC, the nonlinear simulator equations can be used to develop the controller. This section demonstrates how to obtain a linear model for the MPC application using the step test data generated in the prior section.
Predictive Control
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.
Step Testing
Step testing 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 either a doublet test or multiple steps to different levels. 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.
np.savetxt('data_doublet.txt',data,delimiter=',')
np.savetxt('data_doublet.txt',data,delimiter=',', header='Time,Tc,T',comments='')
(:toggle hide step_tests button show="Show Python Step Test Generator":) (:div id=step_tests:)
Additional steps are preferred for systems that show a high degree a nonlinearity or when there is little additional expense to obtain the data. The following code generates data at multiple input levels and with varying different step time intervals. The cooling jacket temperature is not raised above 305 K to avoid reactor instability in open loop.
(:toggle hide step_tests button show="Show Python Step Test Generator":) (:div id=step_tests:)
np.savetxt('data_doublet.txt',data,delimiter=',',\
np.savetxt('cstr_step_tests.txt',data,delimiter=',',\
(:toggle hide python button show="Show Python Simulation Code":) (:div id=python:)
(:toggle hide doublet button show="Show Python Simulation Code":) (:div id=doublet:)
plt.ylabel('Cooling T (K)') plt.legend(['Jacket Temperature'],loc='best')
plt.subplot(3,1,2) plt.plot(t,Ca,'r-',linewidth=3) plt.ylabel('Ca (mol/L)') plt.legend(['Reactor Concentration'],loc='best')
plt.subplot(3,1,3) plt.plot(t,T,'k.-',linewidth=3) plt.ylabel('T (K)') plt.xlabel('Time (min)') plt.legend(['Reactor Temperature'],loc='best')
plt.show() (:sourceend:) (:divend:)
(:toggle hide step_tests button show="Show Python Step Test Generator":) (:div id=step_tests:)
(:source lang=python:) import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint
- define CSTR model
def cstr(x,t,Tc):
Ca = x[0] T = x[1] Tf = 350 Caf = 1.0 q = 100 V = 100 rho = 1000 Cp = 0.239 mdelH = 5e4 EoverR = 8750 k0 = 7.2e10 UA = 5e4 rA = k0*np.exp(-EoverR/T)*Ca dCadt = q/V*(Caf - Ca) - rA dTdt = q/V*(Tf - T) + mdelH/(rho*Cp)*rA + UA/V/rho/Cp*(Tc-T) 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
Tc_ss = 300.0
- Time Interval (min)
t = np.linspace(0,50,501)
- Store results for plotting
Ca = np.ones(len(t)) * Ca_ss T = np.ones(len(t)) * T_ss Tc = np.ones(len(t)) * Tc_ss
- Step cooling temperature
Tc[10:100] = 303.0 Tc[100:200] = 297.0 Tc[200:300] = 300.0 Tc[300:350] = 290.0 Tc[350:400] = 302.0 Tc[400:450] = 302.0 Tc[450:] = 299.0
- Simulate CSTR
for i in range(len(t)-1):
ts = [t[i],t[i+1]] y = odeint(cstr,x0,ts,args=(Tc[i+1],)) 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,Tc,T)) # vertical stack data = data.T # transpose data np.savetxt('data_doublet.txt',data,delimiter=',', header='Time,Tc,T',comments='')
- Plot the results
plt.figure() plt.subplot(3,1,1) plt.plot(t,Tc,'b--',linewidth=3)
cstr_linear_mpc.png
from gekko import gekko
from gekko import GEKKO
m = gekko()
m = GEKKO(remote=True)
m.solve(disp=True,remote=True) # remote=False for local solve
m.solve(disp=True)
(:toggle hide gekko_nmpc button show="Show GEKKO Non-Linear MPC Code":)
(:toggle hide gekko_nmpc button show="Show GEKKO Nonlinear MPC Code":)
from gekko import gekko
from gekko import GEKKO
m = gekko()
m = GEKKO(remote=True)
m.solve(disp=True,remote=True) # remote=False for local solve
m.solve(disp=True)
- m.server = 'https://127.0.0.1'
m.server='https://127.0.0.1'
plt.draw() plt.pause(0.01)
(:sourceend:) (:divend:)
(:toggle hide gekko_nmpc button show="Show GEKKO Non-Linear MPC Code":) (:div id=gekko_nmpc:) (:source lang=python:) import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint from gekko import gekko
- Steady State Initial Condition
u_ss = 280.0
- Feed Temperature (K)
Tf = 350
- Feed Concentration (mol/m^3)
Caf = 1
- Steady State Initial Conditions for the States
Ca_ss = 1 T_ss = 304 x0 = np.empty(2) x0[0] = Ca_ss x0[1] = T_ss
- GEKKO nonlinear MPC
m = gekko() m.server='https://127.0.0.1' m.time = [0,0.02,0.04,0.06,0.08,0.1,0.12,0.15,0.2]
- 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
- initial conditions
Tc0 = 280 T0 = 304 Ca0 = 1.0
tau = m.Const(value=0.5) Kp = m.Const(value=1)
m.Tc = m.MV(value=Tc0,lb=250,ub=350) m.T = m.CV(value=T_ss) m.rA = m.Var(value=0) m.Ca = m.CV(value=Ca_ss)
m.Equation(m.rA == k0*m.exp(-EoverR/m.T)*m.Ca)
m.Equation(m.T.dt() == q/V*(Tf - m.T) + mdelH/(rho*Cp)*m.rA + UA/V/rho/Cp*(m.Tc-m.T))
m.Equation(m.Ca.dt() == q/V*(Caf - m.Ca) - m.rA)
- MV tuning
m.Tc.STATUS = 1 m.Tc.FSTATUS = 0 m.Tc.DMAX = 100 m.Tc.DMAXHI = 20 # constrain movement up m.Tc.DMAXLO = -100 # quick action down
- CV tuning
m.T.STATUS = 1 m.T.FSTATUS = 1 m.T.TR_INIT = 1 m.T.TAU = 1.0 DT = 0.5 # deadband
m.Ca.STATUS = 0 m.Ca.FSTATUS = 0 # no measurement m.Ca.TR_INIT = 0
m.options.CV_TYPE = 1 m.options.IMODE = 6 m.options.SOLVER = 3
- 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
- Time Interval (min)
t = np.linspace(0,8,401)
- Store results for plotting
Ca = np.ones(len(t)) * Ca_ss T = np.ones(len(t)) * T_ss Tsp = np.ones(len(t)) * T_ss u = np.ones(len(t)) * u_ss
- Set point steps
Tsp[0:100] = 330.0 Tsp[100:200] = 350.0 Tsp[200:300] = 370.0 Tsp[300:] = 390.0
- Create plot
plt.figure(figsize=(10,7)) plt.ion() plt.show()
- Simulate CSTR
for i in range(len(t)-1):
# simulate one time period (0.05 sec each loop) ts = [t[i],t[i+1]] y = odeint(cstr,x0,ts,args=(u[i],Tf,Caf)) # retrieve measurements Ca[i+1] = y[-1][0] T[i+1] = y[-1][1] # insert measurement m.T.MEAS = T[i+1] # solve MPC m.solve(disp=True,remote=True) # remote=False for local solve m.T.SPHI = Tsp[i+1] + DT m.T.SPLO = Tsp[i+1] - DT # retrieve new Tc value u[i+1] = m.Tc.NEWVAL # update initial conditions x0[0] = Ca[i+1] x0[1] = T[i+1] # Plot the results plt.clf() plt.subplot(3,1,1) plt.plot(t[0:i],u[0:i],'b--',linewidth=3) plt.ylabel('Cooling T (K)') plt.legend(['Jacket Temperature'],loc='best') plt.subplot(3,1,2) plt.plot(t[0:i],Ca[0:i],'b.-',linewidth=3,label=r'$C_A$') plt.plot([0,t[i-1]],[0.2,0.2],'r--',linewidth=2,label='limit') plt.ylabel(r'$C_A$ (mol/L)') plt.legend(loc='best') plt.subplot(3,1,3) plt.plot(t[0:i],Tsp[0:i],'k-',linewidth=3,label=r'$T_{sp}$') plt.plot(t[0:i],T[0:i],'b.-',linewidth=3,label=r'$T_{meas}$') plt.plot([0,t[i-1]],[400,400],'r--',linewidth=2,label='limit') plt.ylabel('T (K)') plt.xlabel('Time (min)') plt.legend(loc='best')
(:toggle hide python button show="Show Python Simulation Code":) (:div id=python:)
(:toggle hide gekko_mpc button show="Show GEKKO Linear MPC Code":) (:div id=gekko_mpc:)
(:sourceend:) (:divend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/ZvvQ0_PdMPk" frameborder="0" allow="autoplay; encrypted-media" allowfullscreen></iframe> (:htmlend:)
(:toggle hide python button show="Show Python Simulation Code":) (:div id=python:) (:source lang=python:) import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint from gekko import gekko
- Steady State Initial Condition
u_ss = 280.0
- Feed Temperature (K)
Tf = 350
- Feed Concentration (mol/m^3)
Caf = 1
- Steady State Initial Conditions for the States
Ca_ss = 1 T_ss = 304 x0 = np.empty(2) x0[0] = Ca_ss x0[1] = T_ss
- GEKKO linear MPC
m = gekko()
- m.server = 'https://127.0.0.1'
m.time = [0,0.02,0.04,0.06,0.08,0.1,0.15,0.2,0.3,0.4,0.5]
- initial conditions
Tc0 = 280 T0 = 304 Ca0 = 1.0
tau = m.Const(value = 0.5) Kp = m.Const(value = 1)
m.Tc = m.MV(value = Tc0,lb=250,ub=350) m.T = m.CV(value = T_ss)
m.Equation(tau * m.T.dt() == -(m.T - T0) + Kp * (m.Tc - Tc0))
- MV tuning
m.Tc.STATUS = 1 m.Tc.FSTATUS = 0 m.Tc.DMAX = 100 m.Tc.DMAXHI = 5 # constrain movement up m.Tc.DMAXLO = -100 # quick action down
- CV tuning
m.T.STATUS = 1 m.T.FSTATUS = 1 m.T.SP = 330 m.T.TR_INIT = 2 m.T.TAU = 1.0 m.options.CV_TYPE = 2 m.options.IMODE = 6 m.options.SOLVER = 3
- 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
- Time Interval (min)
t = np.linspace(0,10,501)
- Store results for plotting
Ca = np.ones(len(t)) * Ca_ss T = np.ones(len(t)) * T_ss Tsp = np.ones(len(t)) * T_ss u = np.ones(len(t)) * u_ss
- Set point steps
Tsp[0:100] = 330.0 Tsp[100:200] = 350.0 Tsp[200:300] = 370.0 Tsp[300:] = 390.0
- Create plot
plt.figure(figsize=(10,7)) plt.ion() plt.show()
- Simulate CSTR
for i in range(len(t)-1):
# simulate one time period (0.05 sec each loop) ts = [t[i],t[i+1]] y = odeint(cstr,x0,ts,args=(u[i],Tf,Caf)) # retrieve measurements Ca[i+1] = y[-1][0] T[i+1] = y[-1][1] # insert measurement m.T.MEAS = T[i+1] # update setpoint m.T.SP = Tsp[i+1] # solve MPC m.solve(disp=True,remote=True) # remote=False for local solve # change to a fixed starting point for trajectory m.T.TR_INIT = 2 # retrieve new Tc value u[i+1] = m.Tc.NEWVAL # update initial conditions x0[0] = Ca[i+1] x0[1] = T[i+1] # Plot the results plt.clf() plt.subplot(3,1,1) plt.plot(t[0:i],u[0:i],'b--',linewidth=3) plt.ylabel('Cooling T (K)') plt.legend(['Jacket Temperature'],loc='best') plt.subplot(3,1,2) plt.plot(t[0:i],Ca[0:i],'r-',linewidth=3) plt.ylabel('Ca (mol/L)') plt.legend(['Reactor Concentration'],loc='best') plt.subplot(3,1,3) plt.plot(t[0:i],Tsp[0:i],'k-',linewidth=3,label=r'$T_{sp}$') plt.plot(t[0:i],T[0:i],'b.-',linewidth=3,label=r'$T_{meas}$') plt.ylabel('T (K)') plt.xlabel('Time (min)') plt.legend(['Reactor Temperature'],loc='best') plt.draw() plt.pause(0.01)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/tSOMSxGLzQo" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/nqv6jFeVUYA" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/Jxpk4-daDLI" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/sfhHcSF2i90" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/lBx10LvT8uA" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/PyrLMlht-PU" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:toggle hide python button show="Show Python Simulation Code":) (:div id=python:)
(:divend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/Jxpk4-daDLI" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/nqv6jFeVUYA" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/tSOMSxGLzQo" frameborder="0" allowfullscreen></iframe> (:htmlend:)
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.
Step Test (Doublet)
(:source lang=python:) 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 np.savetxt('data_doublet.txt',data,delimiter=',')
- Plot the results
plt.figure() plt.subplot(3,1,1) plt.plot(t,u,'b--',linewidth=3) plt.ylabel('Cooling T (K)') plt.legend(['Jacket Temperature'],loc='best')
plt.subplot(3,1,2) plt.plot(t,Ca,'r-',linewidth=3) plt.ylabel('Ca (mol/L)') plt.legend(['Reactor Concentration'],loc='best')
plt.subplot(3,1,3) plt.plot(t,T,'k.-',linewidth=3) plt.ylabel('T (K)') plt.xlabel('Time (min)') plt.legend(['Reactor Temperature'],loc='best')
plt.show() (:sourceend:)
(:title Model Predictive Control:)
(:title Nonlinear Model Predictive Control:)
(:title Nonlinear Model Predictive Control:) (:keywords Nonlinear, Model Predictive Control, real-time, tutorial:)
(:title Model Predictive Control:) (:keywords Nonlinear, PID, Model Predictive Control, real-time, tutorial:)
Solution
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/qH55ym1g-NM" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/PyrLMlht-PU" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/lBx10LvT8uA" frameborder="0" allowfullscreen></iframe> (:htmlend:)
See additional information on this application on the Process Control Class Web-page.
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/sfhHcSF2i90" frameborder="0" allowfullscreen></iframe> (:htmlend:)
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).
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.
Objective: Design a controller to maintain temperature of a chemical reactor. Develop 3 separate controllers (PID, Linear MPC, Nonlinear MPC) in Python or MATLAB/Simulink. Demonstrate controller performance with steps in the set point and disturbance changes. Estimated time: 3 hours.
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.
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/YHAA-uXhI0E?rel=0" frameborder="0" allowfullscreen></iframe> (:htmlend:)
Objective: Design a model predictive controller to maintain temperature of a chemical reactor. Develop a linear, first-order model of the reactor and implement the controller in Python or MATLAB/Simulink. Demonstrate controller performance with steps in the set point and disturbance changes. Estimated time: 3 hours.
Objective: Design a controller to maintain temperature of a chemical reactor. Develop 3 separate controllers (PID, Linear MPC, Nonlinear MPC) in Python or MATLAB/Simulink. Demonstrate controller performance with steps in the set point and disturbance changes. Estimated time: 3 hours.
(:htmlend:)
(:htmlend:)
Exercise
Objective: Design a model predictive controller to maintain temperature of a chemical reactor. Develop a linear, first-order model of the reactor and implement the controller in Python or MATLAB/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).
Solution
(:htmlend:)
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/qH55ym1g-NM" frameborder="0" allowfullscreen></iframe>
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.
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.
(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/YHAA-uXhI0E?rel=0" frameborder="0" allowfullscreen></iframe> (:htmlend:)
(:title Nonlinear Model Predictive Control:) (:keywords Nonlinear, Model Predictive Control, real-time, tutorial:) (:description Nonlinear Control (NLC) with predictive models is a dynamic optimization approach that seeks to follow a trajectory or drive certain values to maximum or minimum levels:)
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.