Use of Modelica + Python in Process Systems Engineering EducationΒΆ

Basic control of "Seborg reactor"ΒΆ

Bernt LieΒΆ

University College of Southeast NorwayΒΆ

Basic import and definitionsΒΆ

In [1]:
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/"

Basic control designΒΆ

Instantiating models and setting inputsΒΆ

In [2]:
sr_org = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
2018-01-25 11:41:26,236 - OMPython - INFO - OMC Server is up and running at file:///c:/users/bernt_~1/appdata/local/temp/openmodelica.port.f56fce75276643738295fc91dba643c9
In [3]:
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)

Linearizing modelΒΆ

In [4]:
sr_org.getLinearizationOptions()
Out[4]:
{'numberOfIntervals': 500.0,
 'startTime': 0.0,
 'stepSize': 0.002,
 'stopTime': 1.0,
 'tolerance': 1e-08}
In [5]:
A_org, B_org, C_org, D_org = sr_org.linearize()
2018-01-25 11:41:50,404 - OMPython - INFO - OMC Server is up and running at file:///c:/users/bernt_~1/appdata/local/temp/openmodelica.port.dcf837addd2a442fa4c749c9de580e17
In [6]:
A_org
Out[6]:
array([[  4.37955768e+00,   2.09205021e+02],
       [ -3.57142857e-02,  -2.00000000e+00]])
In [7]:
B_org
Out[7]:
array([[  2.09205021e+00,   1.00000000e+00,   1.70530257e-15,
          0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   5.00000000e-03,
          1.00000000e+00]])
In [8]:
C_org
Out[8]:
array([[ 1.,  0.]])
In [9]:
D_org
Out[9]:
array([[ 0.,  0.,  0.,  0.]])
In [10]:
sr_org.getLinearInputs()
Out[10]:
['Tc', 'Ti', 'Vdi', 'cAi']
In [11]:
sr_org.getLinearStates()
Out[11]:
['T', 'cA']
In [12]:
sr_org.getLinearOutputs()
Out[12]:
['y_T']

Theoretical resultsΒΆ

With x=(T,nA), u=(Tc,i,Ti,VΛ™,cA,i) and y=cA, the linearized model is:

dxΞ΄dt=AxΞ΄+BuΞ΄
yΞ΄=CxΞ΄+DuΞ΄
where:
A=(4.37962.0921βˆ’3.5714βˆ’2),B=(2.0921100000.5100)C=(00.01),D=(0000)
This means that OMPython/the Python API has correctly linearized the model.

In [13]:
np.linalg.eig(A_org)
Out[13]:
(array([ 2.83388381, -0.45432613]), array([[ 0.99997271, -0.99973316],
        [-0.00738812,  0.0230998 ]]))

Transfer function from Tc to T (from MATLAB due to some lack of support for control toolbox in Python):

Observe...ΒΆ

At the initial state, the system is open loop unstable due to the eigenvalue at +2.83388381. Observe also that with Tc as the control input/manipulatable variable and T as the controlled output, the system has a zero at 2.092β‹…s+4.184=0β‡’s=βˆ’2; in other words, the zero dynamics is stable!

Tuning P-controllerΒΆ

A P-controller for controlling output y=T via control input u=Tc will have form u=Tcβˆ—+Kp(Trefβˆ’T). This will change the expression for QΛ™rm to:

βˆ’QΛ™rm=VΛ™iVCp(Tiβˆ’T)+UA(Tcβˆ—+Kp(Trefβˆ’T)βˆ’T)

In [14]:
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
In [15]:
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)
In [16]:
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)
Observe...ΒΆ

With Kp=0, we have seen that the system is (open loop) unstable. With Kp=1, 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 Kp=1 is not sufficient go guarantee stability: we need to consider both the effect in temperature and in concentration!

Root locus diagram for system with P-controllerΒΆ

In [17]:
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)
Observe...ΒΆ

The system appears to be stable with a P-controller when Kp>1.2 or so. With Kpβ‰ˆ5.7, 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 Ti∈{2,5} or so (some 10 times closed loop timeconstant of ca. 1/5).

Simulation of system with P-controllerΒΆ

Modelica model for PI control systemΒΆ

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;

Instantiating controlled systemΒΆ

In [18]:
c_sr_org = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.PIconSeborgCSTRorg")
2018-01-25 11:42:03,115 - OMPython - INFO - OMC Server is up and running at file:///c:/users/bernt_~1/appdata/local/temp/openmodelica.port.6b8124ddaec242e98f85cfcf71211afb
In [19]:
pd.DataFrame(c_sr_org.getQuantities())
Out[19]:
Changeable Description Name Value Variability alias aliasvariable
0 false Initializing temperature in reactor, K srORG.T None continuous noAlias None
1 false Initializing concentration of A in reactor, mol/L srORG.cA None continuous noAlias None
2 false der(Initializing temperature in reactor, K) der(srORG.T) None continuous noAlias None
3 false der(Initializing concentration of A in reactor... der(srORG.cA) None continuous noAlias None
4 false None $cse1 None continuous noAlias None
5 false Temperature reference, K T_ref 350.0 continuous noAlias None
6 false Reactor temperature, K _T_org None continuous noAlias None
7 false Nominal cooling temperature, K _Tc_nom None continuous noAlias None
8 false Controlled cooling temperature, K _Tc_org None continuous noAlias None
9 false Influent temperature, K _Ti None continuous noAlias None
10 false Volumetric flow rate through reactor, L/s _Vdi None continuous noAlias None
11 false Molar concentration of A, mol/L _cA_org None continuous noAlias None
12 false Influent molar concentration of A, mol/s _cAi None continuous noAlias None
13 false Heat flow rate, J/min srORG.Qd None continuous noAlias None
14 false Reaction 'constant', ... srORG.k None continuous noAlias None
15 false Rate of reaction, mol/(L.s) srORG.r None continuous noAlias None
16 true Proportional gain in P-controller Kp 0.0 parameter noAlias None
17 true Minimum cooling temperature (freezing point), K Tc_min 273.15 parameter noAlias None
18 true Molar enthalpy of reaction, J/mol srORG.DrHt -50000.0 parameter noAlias None
19 true Activation temperature, K srORG.EdR 8750.0 parameter noAlias None
20 true Initial temperature, K srORG.T0 350.0 parameter noAlias None
21 true Heat transfer parameter, J/(min.K) srORG.UA 50000.0 parameter noAlias None
22 true Reactor volume, L srORG.V 100.0 parameter noAlias None
23 true Stoichiometric constant, - srORG.a 1.0 parameter noAlias None
24 true Initial concentration of A, mol/L srORG.cA0 0.5 parameter noAlias None
25 true Specific heat capacity of mixture, J.g-1.K-1 srORG.cph 0.239 parameter noAlias None
26 false Pre-exponential factor, 1/min srORG.k0 None parameter noAlias None
27 true Liquid density, g/L srORG.rho 1000.0 parameter noAlias None
28 false Cooling temperature', K srORG.Tc None continuous alias _Tc_org
29 false Influent temperature, K srORG.Ti None continuous alias _Ti
30 false Volumetric flow rate through reactor, L/min srORG.Vdi None continuous alias _Vdi
31 false Influent molar concentration of A, mol/L srORG.cAi None continuous alias _cAi
32 false Reactor temperature, K srORG.y_T None continuous alias srORG.T
In [20]:
c_sr_org.setSimulationOptions(stopTime=15,stepSize=0.05)

Parameters & SimulationΒΆ

In [21]:
c_sr_org.getParameters()
Out[21]:
{'Kp': 0.0,
 'Tc_min': 273.15,
 'srORG.DrHt': -50000.0,
 'srORG.EdR': 8750.0,
 'srORG.T0': 350.0,
 'srORG.UA': 50000.0,
 'srORG.V': 100.0,
 'srORG.a': 1.0,
 'srORG.cA0': 0.5,
 'srORG.cph': 0.239,
 'srORG.k0': None,
 'srORG.rho': 1000.0}

Plotting resultsΒΆ

In [22]:
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)
In [23]:
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)
In [24]:
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)
Observe...ΒΆ

The design value for Kp 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 Tc>4 [C].

In [25]:
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)
In [26]:
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)
In [27]:
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)
Observe...ΒΆ

The constraint in the control input dramatically reduces the performance of the controller.