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

Nonlinear 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/"

Nonlinear Decoupling control designΒΆ

Control strategyΒΆ

The control strategy consists in specifying the closed loop response, in this particular case as:

dTdt=1Ο„cl(Trefβˆ’T)
while keeping the control input Tc indetermined. The controlled system is thus a DAE system, with two differential equations for temperature T and none for Tc -- the DAE solver will then compute Tc so that the two differential equations for T simultaneously are satisfied. Parameter Ο„cl is the specified closed loop time constant.

This strategy is sometimes referred to as Nonlinear Decoupling, in that the closed loop response often is specified as

dydt=Ξ›(yrefβˆ’y)
where Ξ› often is diagonal -- thus decoupled.

The Nonlinear Decoupling method only works if the system has stable zero dynamics. In our case, analysis of the linearized system indicated that locally, at the operating point, the system had a negative zero in the transfer function, hence stable zero dynamics.

The following is an example of Modelica code for specifying the controller:

mo
model NDconSeborgCSTRorg
    // Simulation of Seborg Reactor
    // author:  Bernt Lie
    //          University of Southeast Norway
    //          January 11, 2018
    //
    // Instantiate model of Seborg Reactor (sr)
    ModSeborgCSTRorg srORG;
    // Declaring parameters
    parameter Real T_cl = 2 "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 "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;
    T_ref = if time < 3 then 350 else 360;
    // -- response shaping
    der(_T_org) = (T_ref - _T_org)/T_cl;
    // -- injecting input functions to model inputs
    srORG.Vdi = _Vdi;
    srORG.cAi = _cAi;
    srORG.Ti = _Ti;
    srORG.Tc = _Tc;
    // -- outputs
    _cA_org = srORG.cA;
    _T_org = srORG.T;
    _Tc_org = srORG.Tc;
  end NDconSeborgCSTRorg;

Instantiating models and setting inputsΒΆ

In [2]:
c_sr_org = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.NDconSeborgCSTRorg")
2018-01-25 11:43:39,187 - OMPython - INFO - OMC Server is up and running at file:///c:/users/bernt_~1/appdata/local/temp/openmodelica.port.4df6e81e26b440eb89e1f5bfb143ee85
In [3]:
c_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)

Parameters & SimulationΒΆ

In [4]:
c_sr_org.getParameters()
Out[4]:
{'T_cl': 2.0,
 '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 [5]:
#c_sr_org.setParameters(Tc_min = -1e3)
c_sr_org.setParameters(T_cl=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=Cr1,label=r"$T; \tau_\mathrm{cl} = 2$")
#
c_sr_org.setParameters(T_cl=1)
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,LS=LS2,label=r"$T; \tau_\mathrm{cl} = 1$")
#
c_sr_org.setParameters(T_cl=0.5)
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; \tau_\mathrm{cl} = 0.5$")
#
c_sr_org.setParameters(T_cl=0.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,LS=LS2,label=r"$T; \tau_\mathrm{cl} = 0.2$")
#
plt.plot(tm_corg,T_ref-273.15,LW=LW1,color=Cb1,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_NDconSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
In [6]:
c_sr_org.setParameters(T_cl=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=Cr1,label=r"$c_\mathrm{A}; \tau_\mathrm{cl} = 2$")
#
c_sr_org.setParameters(T_cl=1)
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,LS=LS2,label=r"$c_\mathrm{A}; \tau_\mathrm{cl} = 1$")
#
c_sr_org.setParameters(T_cl=0.5)
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}; \tau_\mathrm{cl} = 0.5$")
#
c_sr_org.setParameters(T_cl=0.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,LS=LS2,label=r"$c_\mathrm{A}; \tau_\mathrm{cl} = 0.2$")
#
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_NDconSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
In [7]:
c_sr_org.setParameters(T_cl=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=Cr1,label=r"$T_\mathrm{c}; \tau_\mathrm{cl} = 2$")
#
c_sr_org.setParameters(T_cl=1)
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,LS=LS2,label=r"$T_\mathrm{c}; \tau_\mathrm{cl} = 1$")
#
c_sr_org.setParameters(T_cl=0.5)
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}; \tau_\mathrm{cl} = 0.5$")
#
c_sr_org.setParameters(T_cl=0.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,LS=LS2,label=r"$T_\mathrm{c}; \tau_\mathrm{cl} = 0.2$")
#
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 Nonlinear Decoupling controller gives "perfect" compensation of the nonlinearity. Even with relatively fast closed-loop performance, the control input is limited and never approaches freezing.