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/"
srp = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
srm = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
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)
p = srp.getParameters()
p.keys()
p.values()
UA = p["UA"]
UA
nr.seed(0)
UUAA = UA*(1 + (nr.rand(20)-0.5)/2)
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)
sr_org = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
sr_org.setSimulationOptions(stopTime=20)
sr_org.setInputs(Tc=300)
sr_org.setInputs(Ti=350)
sr_org.setInputs(Vdi=100)
sr_org.setInputs(cAi=1)
sr_org.simulate()
tm = sr_org.getSolutions("time")
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
Sn, Sa = sensitivity(sr_org,["T","cA"],["UA","EdR"])
Sn
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)
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.