from OMPython import ModelicaSystem
import numpy as np
import numpy.random as nr
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
LW1 = 2.5
LW2 = LW1/2
Cb1 = (0.3,0.3,1)
Cb2 = (0.7,0.7,1)
Cg1 = (0,0.6,0)
Cg2 = (0.5,0.8,0.5)
Cr1 = "Red"
Cr2 = (1,0.5,0.5)
LS1 = "solid"
LS2 = "dotted"
LS3 = "dashed"
figpath = "../figs/"
sr_org = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
sr_org.setSimulationOptions(stopTime=15,stepSize=0.05)
#
sr_org.setInputs(Tc=300)
sr_org.setInputs(Ti=350)
sr_org.setInputs(Vdi=100)
sr_org.setInputs(cAi=1)
sr_org.getLinearizationOptions()
A_org, B_org, C_org, D_org = sr_org.linearize()
A_org
B_org
C_org
D_org
sr_org.getLinearInputs()
sr_org.getLinearStates()
sr_org.getLinearOutputs()
With , and , the linearized model is:
np.linalg.eig(A_org)
Transfer function from to (from MATLAB due to some lack of support for control toolbox in Python):
At the initial state, the system is open loop unstable due to the eigenvalue at +2.83388381. Observe also that with as the control input/manipulatable variable and as the controlled output, the system has a zero at ; in other words, the zero dynamics is stable!
A P-controller for controlling output via control input will have form . This will change the expression for to:
EdR = 8750
k0 = np.exp(EdR/350)
VdidV = 1
cAi = 1
def cA(T):
return VdidV/(k0*np.exp(-EdR/T) + VdidV)*cAi
#
rho = 1e3
V = 100
cph = 0.239
Cp = rho*V*cph
DrHt = -5e4
UA = 5e4
Ti = 350
Tcast = 300
Tref = 350
Kp = 0
def Qdrm(T,Kp):
return -(VdidV*Cp*(Ti-T) + UA*(Tcast+Kp*(Tref-T)-T))
def Qdgen(T):
return (-DrHt)*k0*np.exp(-EdR/T)*cA(T)*V
TT = np.linspace(35,110)
plt.plot(TT, Qdgen(TT+273.15)/1e6,LW=LW1, color=Cr1, label=r"$\dot{Q}_\mathrm{gen}$")
plt.plot(TT, Qdrm(TT+273.15,Kp)/1e6,LW=LW1, color=Cb1, label=r"$\dot{Q}_\mathrm{rm}$, $K_\mathrm{p} = 0$")
plt.plot(TT, Qdrm(TT+273.15,1)/1e6,LW=LW1, color=Cb2, label=r"$\dot{Q}_\mathrm{rm}$, $K_\mathrm{p} = 1$")
plt.plot(TT, Qdrm(TT+273.15,3)/1e6,LW=LW1, color=Cb2, LS = LS2, label=r"$\dot{Q}_\mathrm{rm}$, $K_\mathrm{p} = 3$")
plt.plot(TT, Qdrm(TT+273.15,5)/1e6,LW=LW1, color=Cb1, LS = LS2, label=r"$\dot{Q}_\mathrm{rm}$, $K_\mathrm{p} = 5$")
#
plt.plot((350-273.15)*np.array([1,1]), [-1,6],LW=LW1,color=Cg1,LS=LS2)
plt.plot([30,110], Qdrm(350,Kp)/1e6*np.array([1,1]),LW=LW1,color=Cg1,LS=LS2)
plt.plot((50.8)*np.array([1,1]), [-1,6],LW=LW2,color=Cg2,LS=LS2)
plt.plot([30,110], Qdrm(50.8+273.15,Kp)/1e6*np.array([1,1]),LW=LW2,color=Cg2,LS=LS2)
plt.plot((97.2)*np.array([1,1]), [-1,6],LW=LW2,color=Cg2,LS=LS2)
plt.plot([30,110], Qdrm(97.2+273.15,Kp)/1e6*np.array([1,1]),LW=LW2,color=Cg2,LS=LS2)
#
plt.title("Generated vs. removed heat rates")
plt.xlabel(r"temperature $T$ [C]")
plt.ylabel(r"$\dot{Q}$ [MJ/min]")
plt.grid()
plt.xlim(35,110)
plt.ylim(np.min(Qdrm(TT+273.15,Kp)/1e6), np.max(Qdrm(TT+273.15,Kp)/1e6))
plt.legend()
figfile = "equil_controlled_SeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
plt.plot(TT, Qdgen(TT+273.15)/1e6,LW=LW1, color=Cr1, label=r"$\dot{Q}_\mathrm{gen}$")
plt.plot(TT, Qdrm(TT+273.15,Kp)/1e6,LW=LW1, color=Cb1, label=r"$\dot{Q}_\mathrm{rm}$, $K_\mathrm{p} = 0$")
plt.plot(TT, Qdrm(TT+273.15,1)/1e6,LW=LW1, color=Cb2, label=r"$\dot{Q}_\mathrm{rm}$, $K_\mathrm{p} = 1$")
plt.plot(TT, Qdrm(TT+273.15,3)/1e6,LW=LW1, color=Cb2, LS = LS2, label=r"$\dot{Q}_\mathrm{rm}$, $K_\mathrm{p} = 3$")
plt.plot(TT, Qdrm(TT+273.15,5)/1e6,LW=LW1, color=Cb1, LS = LS2, label=r"$\dot{Q}_\mathrm{rm}$, $K_\mathrm{p} = 5$")
#
plt.plot((350-273.15)*np.array([1,1]), [-1,6],LW=LW1,color=Cg1,LS=LS2)
plt.plot([30,110], Qdrm(350,Kp)/1e6*np.array([1,1]),LW=LW1,color=Cg1,LS=LS2)
plt.plot((50.8)*np.array([1,1]), [-1,6],LW=LW2,color=Cg2,LS=LS2)
plt.plot([30,110], Qdrm(50.8+273.15,Kp)/1e6*np.array([1,1]),LW=LW2,color=Cg2,LS=LS2)
plt.plot((97.2)*np.array([1,1]), [-1,6],LW=LW2,color=Cg2,LS=LS2)
plt.plot([30,110], Qdrm(97.2+273.15,Kp)/1e6*np.array([1,1]),LW=LW2,color=Cg2,LS=LS2)
#
plt.title("Generated vs. removed heat rates")
plt.xlabel(r"temperature $T$ [C]")
plt.ylabel(r"$\dot{Q}$ [MJ/min]")
plt.grid()
plt.xlim(76,78)
plt.ylim(2.2,2.8)
plt.legend()
figfile = "equil_detail_controlled_SeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
With , we have seen that the system is (open loop) unstable. With , the system appears to have a stable equilibrium point, and the larger the value for the gain is, the larger the stability margin becomes. (This is contrary to open-loop stable systems.) However, as we will see below: the choice is not sufficient go guarantee stability: we need to consider both the effect in temperature and in concentration!
A_ol = A_org
B_ol = np.zeros((A_ol.shape[0],1))
B_ol[:,0] = B_org[:,0]
C_ol = C_org
#
KKp = np.linspace(-1,1.2)
Eig = np.zeros((KKp.size,2),np.dtype(complex))
#
for (i,Kp) in enumerate(KKp):
A_cl = A_org - Kp*B_ol*C_ol
Eig[i,:] = np.linalg.eig(A_cl)[0]
#
plt.plot(np.real(Eig)[:,0],np.imag(Eig)[:,0],"x",color=Cr1,markersize=5)
plt.plot(np.real(Eig)[:,1],np.imag(Eig)[:,1],"x",color=Cr1,markersize=5, label=r"$K_\mathrm{p} < 1.2$")
#
KKp = np.linspace(1.2,8)
Eig = np.zeros((KKp.size,2),np.dtype(complex))
#
for (i,Kp) in enumerate(KKp):
A_cl = A_org - Kp*B_ol*C_ol
Eig[i,:] = np.linalg.eig(A_cl)[0]
#
plt.plot(np.real(Eig)[:,0],np.imag(Eig)[:,0],"+",color=Cb1,markersize=5)
plt.plot(np.real(Eig)[:,1],np.imag(Eig)[:,1],"+",color=Cb1,markersize=5, label=r"$K_\mathrm{p} > 1.2$")
#
Eig0 = np.linalg.eig(A_ol)[0]
plt.plot(np.real(Eig0)[0],np.imag(Eig0)[0],"x",color=Cg1,markersize=8)
plt.plot(np.real(Eig0)[1],np.imag(Eig0)[1],"x",color=Cg1,markersize=8, label=r"$K_\mathrm{p} = 0$")
#
Eig1 = np.linalg.eig(A_org - 5.7*B_ol*C_ol)[0]
plt.plot(np.real(Eig1)[0],np.imag(Eig1)[0],"o",color=Cg1,markersize=8)
plt.plot(np.real(Eig1)[1],np.imag(Eig1)[1],"o",color=Cg1,markersize=8, label=r"$K_\mathrm{p} = 5.7$")
#
plt.grid()
plt.xlabel(r"$\Re [ \lambda (A_\mathrm{cl}) ]$")
plt.ylabel(r"$\Im [ \lambda (A_\mathrm{cl}) ]$")
plt.title(r"Closed loop eigenvalues as a function of $K_\mathrm{p}\in\{-1,8\}$")
plt.legend()
figfile = "ClosedLoop_Eigenvals_SeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
The system appears to be stable with a P-controller when or so. With , both eigenvalues become real with similar values. Some modification may be needed if we add integral action, but probably not much. We could probably choose or so (some 10 times closed loop timeconstant of ca. ).
mo
model PIconSeborgCSTRorg
// Simulation of Seborg Reactor
// author: Bernt Lie
// University of Southeast Norway
// January 10, 2018
//
// Instantiate model of Seborg Reactor (sr)
ModSeborgCSTRorg srORG;
// Declaring parameters
parameter Real Kp = 0 "Proportional gain in P-controller";
parameter Real Tc_min = 273.15 "Minimum cooling temperature (freezing point), K";
// Declaring variables
// -- inputs
Real _Vdi "Volumetric flow rate through reactor, L/s";
Real _cAi "Influent molar concentration of A, mol/s";
Real _Ti "Influent temperature, K";
Real _Tc_nom "Nominal cooling temperature, K";
Real T_ref "Temperature reference, K";
// -- outputs
output Real _cA_org "Molar concentration of A, mol/L";
output Real _T_org "Reactor temperature, K";
output Real _Tc_org "Controlled cooling temperature, K";
// Equations
equation
// -- input values
_Vdi = 100;
_cAi = 1;
_Ti = 350;
_Tc_nom = 300;
T_ref = if time < 3 then 350 else 360;
// -- injecting input functions to model inputs
srORG.Vdi = _Vdi;
srORG.cAi = _cAi;
srORG.Ti = _Ti;
srORG.Tc = max(Tc_min,_Tc_nom + Kp*(T_ref - _T_org));
// -- outputs
_cA_org = srORG.cA;
_T_org = srORG.T;
_Tc_org = srORG.Tc;
end PIconSeborgCSTRorg;
c_sr_org = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.PIconSeborgCSTRorg")
pd.DataFrame(c_sr_org.getQuantities())
c_sr_org.setSimulationOptions(stopTime=15,stepSize=0.05)
c_sr_org.getParameters()
c_sr_org.setParameters(Tc_min = -1e3)
c_sr_org.setParameters(Kp=0)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,T_corg-273.15,LW=LW1,color=Cr1,label=r"$T; K_\mathrm{p} = 0$")
#
c_sr_org.setParameters(Kp=2)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,T_corg-273.15,LW=LW1,color=Cr2,label=r"$T; K_\mathrm{p} = 2$")
#
c_sr_org.setParameters(Kp=4)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,T_corg-273.15,LW=LW1,color=Cb1,label=r"$T; K_\mathrm{p} = 4$")
#
c_sr_org.setParameters(Kp=8)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,T_corg-273.15,LW=LW1,color=Cg1,label=r"$T; K_\mathrm{p} = 8$")
#
plt.plot(tm_corg,T_ref-273.15,LW=LW1,color=Cr1,LS=LS2,label=r"$T_{\mathrm{ref}}$")
#
plt.title("Reactor temperature")
plt.xlabel(r"time $t$ [min]")
plt.ylabel(r"$T$ [${}^\circ \mathrm{C}$]")
plt.legend()
plt.grid()
plt.xlim(0,15)
figfile = "T_PIconSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
c_sr_org.setParameters(Kp=0)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,cA_corg,LW=LW1,color=Cr1,label=r"$c_\mathrm{A}; K_\mathrm{p} = 0$")
#
c_sr_org.setParameters(Kp=2)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,cA_corg,LW=LW1,color=Cr2,label=r"$c_\mathrm{A}; K_\mathrm{p} = 2$")
#
c_sr_org.setParameters(Kp=4)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,cA_corg,LW=LW1,color=Cb1,label=r"$c_\mathrm{A}; K_\mathrm{p} = 4$")
#
c_sr_org.setParameters(Kp=8)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,cA_corg,LW=LW1,color=Cg1,label=r"$c_\mathrm{A}; K_\mathrm{p} = 8$")
#
plt.title("Reactor concentration ")
plt.xlabel(r"time $t$ [min]")
plt.ylabel(r"$c_\mathrm{A}$ [mol/L]")
plt.legend()
plt.grid()
plt.xlim(0,15)
figfile = "cA_PIconSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
c_sr_org.setParameters(Kp=0)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,Tc_corg-273.15,LW=LW1,color=Cr1,label=r"$T_\mathrm{c}; K_\mathrm{p} = 0$")
#
c_sr_org.setParameters(Kp=2)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,Tc_corg-273.15,LW=LW1,color=Cr2,label=r"$T_\mathrm{c}; K_\mathrm{p} = 2$")
#
c_sr_org.setParameters(Kp=4)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,Tc_corg-273.15,LW=LW1,color=Cb1,label=r"$T_\mathrm{c}; K_\mathrm{p} = 4$")
#
c_sr_org.setParameters(Kp=8)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,Tc_corg-273.15,LW=LW1,color=Cg1,label=r"$T_\mathrm{c}; K_\mathrm{p} = 8$")
#
plt.title("Control signal: cooling temperature ")
plt.xlabel(r"time $t$ [min]")
plt.ylabel(r"$T_\mathrm{c}$ [${}^\circ \mathrm{C}$]")
plt.legend()
plt.grid()
plt.xlim(0,15)
figfile = "Tc_PIconSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
The design value for based on a linear model does not give good performance; a larger control gain should be chosen to ensure decent performance for the nonlinear system. We also see that the controllers give negative cooling temperature, which is not realistic. Let us see what happens if we constrain the cooling temperature (the control input) to [C].
c_sr_org.setParameters(Tc_min = 4+273.15)
c_sr_org.setParameters(Kp=0)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,T_corg-273.15,LW=LW1,color=Cr1,label=r"$T; K_\mathrm{p} = 0$")
#
c_sr_org.setParameters(Kp=2)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,T_corg-273.15,LW=LW1,color=Cr2,label=r"$T; K_\mathrm{p} = 2$")
#
c_sr_org.setParameters(Kp=4)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,T_corg-273.15,LW=LW1,color=Cb1,label=r"$T; K_\mathrm{p} = 4$")
#
c_sr_org.setParameters(Kp=8)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,T_corg-273.15,LW=LW1,color=Cg1,label=r"$T; K_\mathrm{p} = 8$")
#
plt.plot(tm_corg,T_ref-273.15,LW=LW1,color=Cr1,LS=LS2,label=r"$T_{\mathrm{ref}}$")
#
plt.title("Reactor temperature")
plt.xlabel(r"time $t$ [min]")
plt.ylabel(r"$T$ [${}^\circ \mathrm{C}$]")
plt.legend()
plt.grid()
plt.xlim(0,15)
figfile = "T_constr_PIconSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
c_sr_org.setParameters(Kp=0)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,cA_corg,LW=LW1,color=Cr1,label=r"$c_\mathrm{A}; K_\mathrm{p} = 0$")
#
c_sr_org.setParameters(Kp=2)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,cA_corg,LW=LW1,color=Cr2,label=r"$c_\mathrm{A}; K_\mathrm{p} = 2$")
#
c_sr_org.setParameters(Kp=4)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,cA_corg,LW=LW1,color=Cb1,label=r"$c_\mathrm{A}; K_\mathrm{p} = 4$")
#
c_sr_org.setParameters(Kp=8)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,cA_corg,LW=LW1,color=Cg1,label=r"$c_\mathrm{A}; K_\mathrm{p} = 8$")
#
plt.title("Reactor concentration ")
plt.xlabel(r"time $t$ [min]")
plt.ylabel(r"$c_\mathrm{A}$ [mol/L]")
plt.legend()
plt.grid()
plt.xlim(0,15)
figfile = "cA_constr_PIconSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
c_sr_org.setParameters(Kp=0)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,Tc_corg-273.15,LW=LW1,color=Cr1,label=r"$T_\mathrm{c}; K_\mathrm{p} = 0$")
#
c_sr_org.setParameters(Kp=2)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,Tc_corg-273.15,LW=LW1,color=Cr2,label=r"$T_\mathrm{c}; K_\mathrm{p} = 2$")
#
c_sr_org.setParameters(Kp=4)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,Tc_corg-273.15,LW=LW1,color=Cb1,label=r"$T_\mathrm{c}; K_\mathrm{p} = 4$")
#
c_sr_org.setParameters(Kp=8)
c_sr_org.simulate()
tm_corg, T_corg, Tc_corg, cA_corg, T_ref = c_sr_org.getSolutions("time","_T_org","_Tc_org","_cA_org","T_ref")
plt.plot(tm_corg,Tc_corg-273.15,LW=LW1,color=Cg1,label=r"$T_\mathrm{c}; K_\mathrm{p} = 8$")
#
plt.title("Control signal: cooling temperature ")
plt.xlabel(r"time $t$ [min]")
plt.ylabel(r"$T_\mathrm{c}$ [${}^\circ \mathrm{C}$]")
plt.legend()
plt.grid()
plt.xlim(0,15)
figfile = "Tc_constr_PIconSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
The constraint in the control input dramatically reduces the performance of the controller.