Use of Modelica + Python in Process Systems Engineering Education

Sensitivity analysis 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/"

Reactor temperature sensitivity to UA from Monte Carlo simulations

Instantiating models and setting inputs

In [2]:
srp = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
srm = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
2018-01-25 11:37:32,867 - OMPython - INFO - OMC Server is up and running at file:///c:/users/bernt_~1/appdata/local/temp/openmodelica.port.0c15a45e519643eca65546bb8b7fdeae
2018-01-25 11:37:40,729 - OMPython - INFO - OMC Server is up and running at file:///c:/users/bernt_~1/appdata/local/temp/openmodelica.port.daedfb72664f4ddcb0871737b3e0be74
In [3]:
srp.setSimulationOptions(stopTime=10,stepSize=0.05)
srm.setSimulationOptions(stopTime=10,stepSize=0.05)
#
srp.setInputs(Tc=300+5)
srm.setInputs(Tc=300-10)
srp.setInputs(Ti=350,Vdi=100,cAi=1)
srm.setInputs(Ti=350,Vdi=100,cAi=1)

Finding parameter names, and generating random variations of UA with uniform distribution and 5% variation

In [4]:
p = srp.getParameters()
p.keys()
Out[4]:
['a', 'cph', 'DrHt', 'cA0', 'T0', 'EdR', 'k0', 'rho', 'V', 'UA']
In [5]:
p.values()
Out[5]:
[1.0, 0.239, -50000.0, 0.5, 350.0, 8750.0, None, 1000.0, 100.0, 50000.0]
In [6]:
UA = p["UA"]
UA
Out[6]:
50000.0
In [7]:
nr.seed(0)
UUAA = UA*(1 + (nr.rand(20)-0.5)/2)

Simulating for nominal values UA¯ of UA, as well as for random variations in UA

In [8]:
srp.simulate()
tm0p,T0p = srp.getSolutions("time","T")
srm.simulate()
tm0m,T0m = srm.getSolutions("time","T")
#
for ua in UUAA:
    srp.setParameters(UA=ua);
    srp.simulate()
    tm0p1,T0p1 = srp.getSolutions("time","T")
    srm.setParameters(UA=ua);
    srm.simulate()
    tm0m1,T0m1 = srm.getSolutions("time","T")
    plt.plot(tm0p1,T0p1-273.15,LW=LW2,color=Cr1,LS=LS2,label="_nolabel_")
    plt.plot(tm0m1,T0m1-273.15,LW=LW2,color=Cb1,LS=LS2,label="_nolabel_")
#
plt.plot(tm0p,T0p-273.15,LW=LW1,color=Cr1,label=r"$\overline{UA} : T_c^{\ast}+5$")
plt.plot(tm0m,T0m-273.15,LW=LW1,color=Cb1,label=r"$\overline{UA} : T_c^{\ast}-10$")
plt.title("Reactor temperature Monte Carlo study")
plt.xlabel(r"time $t$ [s]")
plt.ylabel(r"$T$ [${}^\circ \mathrm{C}$]")
plt.grid()
plt.xlim(0,10)
plt.legend()
figfile = "sensitivitySeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)

Sensitivity of quantities wrt. parameters

In [9]:
sr_org = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
2018-01-25 11:38:12,302 - OMPython - INFO - OMC Server is up and running at file:///c:/users/bernt_~1/appdata/local/temp/openmodelica.port.8161dd3a3c9f43b09eb6f225aab4fc57
In [10]:
sr_org.setSimulationOptions(stopTime=20)
In [11]:
sr_org.setInputs(Tc=300)
sr_org.setInputs(Ti=350)
sr_org.setInputs(Vdi=100)
sr_org.setInputs(cAi=1)
In [12]:
sr_org.simulate()
tm = sr_org.getSolutions("time")

Function for numerical sensitivity computation

In [13]:
def sensitivity(obj,Lv,Lp,Le=[1e-2]):
    """
    Method for computing numeric sensitivity of OpenModelica object
    
    Arguments:
    ----------
    1st arg: obj # OMPython API model object -- should be removed in method?
    2nd arg: Lv  # List of strings of Modelica Variable names
    3rd arg: Lp  # List of strings of Modelica Parameter names
    4th arg: Le  # List of float Excitations of parameters; defaults to single 1e-2
    
    Returns:
    --------
    1st return: LSname # List of Sensitivity names
    2nd return: Sarray # 2D array of sensitivities, one row per sensitivity name
    """
    # Production quality code should check type and form of input arguments
    #
    nLp = len(Lp) # number of parameter names
    nLe = len(Le) # number of excitations in parameters
    # Adjusting size of Le to that of Lp
    if nLe < nLp:
        for i in range(nLe,nLp):
            Le.append(Le[nLe-1]) # expands Le with the last element of Le 
    elif nLe > nLp:
        Le = Le[0:nLp] # truncates Le to same length as Lp
    # Nominal parameters p0
    par0 = obj.getParameters(*Lp)
    # eXcitation parameters parX
    parX = []
    for i,p in enumerate(par0):
        parX.append(p*(1.+Le[i]))
    # Zip parameter names and parameter values into list of tuples
    # --- preparation for setting excitated parameters via keyword assignment
    Lpar0 = zip(Lp,par0) # List of parameter name/value pairs for resetting to nominal case
    LparX = zip(Lp,parX) # List of eXcited parameter name/value pairs
    # Simulate nominal system
    obj.simulate()
    # Get nominal SOLutions of variables of interest (Lv), converted to 2D array
    sol0 = np.asarray(obj.getSolutions(*Lv))
    # Get list of eXcited SOLutions (2D arrays), one for each parameter (Lp)
    solX = []
    for i,d in enumerate(LparX):
        # change to excited parameter
        obj.setParameters(**dict([d]))
        # simulate perturbed system
        obj.simulate()
        # get eXcited SOLutions (Lv) as 2D array, and append to list
        solX.append(np.asarray(obj.getSolutions(*Lv)))
        # reset parameters to nominal values
        obj.setParameters(**dict(Lpar0))
    # Compute sensitivities and add to list, one 2D array per parameter (Lp)
    LSname = [] # Preparing list of names of sensitivities
    LSarray = [] # Preparing list of 2d sensitivity arrays, 1 per parameter
    for i,sol in enumerate(solX):
        LSarray.append((sol-sol0)/(par0[i]*Le[i]))
        for var in Lv:
            # NOTE: In analytic sensitivities, "var" and "LP[i]" are reversed
            # ... I find the order below more natural
            LSname.append("$Sensitivities." + var + "." + Lp[i])
    # Converting list of 2D arrays to a 2D array by stacking elements in list
    Sarray = np.vstack(LSarray)
    return LSname, Sarray
In [14]:
Sn, Sa = sensitivity(sr_org,["T","cA"],["UA","EdR"])
In [15]:
Sn
Out[15]:
['$Sensitivities.T.UA',
 '$Sensitivities.cA.UA',
 '$Sensitivities.T.EdR',
 '$Sensitivities.cA.EdR']
In [16]:
plt.subplot(2,1,1)
plt.plot(tm,Sa[0],linewidth=LW1,color=Cb1)
plt.ylabel(r'${\partial c_\mathrm{A}}/{\partial \mathrm{UA}}$')
plt.title(r'CSTR: Sensitivities in $\mathrm{UA}$')
plt.grid()
plt.xlim(0,20)
plt.subplot(2,1,2)
plt.plot(tm,Sa[1]*1e3,linewidth=LW1,color=Cr1)
plt.ylabel(r'$10^3 \times{\partial T}/{\partial \mathrm{UA}}$')
plt.xlabel(r'$t $ [min]')
plt.grid()
plt.xlim(0,20)
figfile = "TcA_UA_NumSensitivityReactororg.pdf"
plt.savefig(figpath+figfile)
In [17]:
plt.subplot(2,1,1)
plt.plot(tm,Sa[2],linewidth=LW1,color=Cb1)
plt.ylabel(r'${\partial c_\mathrm{A}}/{\partial (E/R)}$')
plt.title(r'CSTR: Sensitivities in $E/R$')
plt.grid()
plt.xlim(0,20)
plt.subplot(2,1,2)
plt.plot(tm,Sa[3]*1e3,linewidth=LW1,color=Cr1)
plt.ylabel(r'$10^3 \times{\partial T}/{\partial (E/R)}$')
plt.xlabel(r'$t $ [min]')
plt.grid()
plt.xlim(0,20)
figfile = "TcA_EdR_NumSensitivitySeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)

Here, numerical sensitivity has been computed. OpenModelica has support for computing analytic sensitivities, but that has not been properly interfaced to the Python API at the moment.