8.2. Simulation, Control, and Estimation using Pyomo#
8.2.1. Installations#
The following instructions show how to download and install pyomo and the ipopt solver. Execute the appropriate cell for your platform (if needed).
8.2.1.1. Google Colab#
!pip install -q pyomo
!wget -N -q "https://ampl.com/dl/open/ipopt/ipopt-linux64.zip"
!unzip -o -q ipopt-linux64
ipopt_executable = '/content/ipopt'
/bin/sh: wget: command not found
unzip: cannot find or open ipopt-linux64, ipopt-linux64.zip or ipopt-linux64.ZIP.
8.2.1.2. MacOS#
!pip install -q pyomo
!curl -s https://ampl.com/dl/open/ipopt/ipopt-osx.zip --output ipopt-osx.zip
!tar xf ipopt-osx.zip ipopt
ipopt_executable = "./ipopt"
!rm ipopt-osx.zip
8.2.1.3. Windows PC#
!conda install -c conda-forge pyomo pyomo.extras
!conda install -c conda-forge/label/cf201901 ipopt
Collecting package metadata: ...working... done
Solving environment: ...working... done
# All requested packages already installed.
Collecting package metadata: ...working... done
Solving environment: ...working... done
# All requested packages already installed.
8.2.2. Process Information#
8.2.2.1. Process Parameter Values#
P = 0.04 # power input when the system is turned
Ua = 0.068 # heat transfer coefficient from heater to environment
CpH = 6.50 # heat capacity of the heater (J/deg C)
CpS = 1.25 # heat capacity of the sensor (J/deg C)
Uc = 0.036 # heat transfer coefficient from heater to sensor
Tamb = 21.0 # ambient room temperature
8.2.2.2. Process Inputs#
The next cell defines some process inputs that will be used throughout the notebook to demonstrate aspects of process simulation, control, and estimation. These are gathered in one place to make it easier to modify the notebook to test the response under different conditions. These functions are implemented using the interp1d
from the scipy
library.
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import interpolate
import numpy as np
tclab_disturbance = interpolate.interp1d(
[ 0, 300, 400, 9999],
[ 0, 0, -.5, -.5])
tclab_input = interpolate.interp1d(
[ 0, 50, 51, 450, 451, 9999],
[ 0, 0, 80, 80, 25, 25])
tclab_setpoint = interpolate.interp1d(
[0, 50, 150, 450, 550, 9999],
[Tamb, Tamb, 60, 60, 35, 35])
t_sim = np.linspace(0, 1000, 201)
u_sim = tclab_input(t_sim)
d_sim = tclab_disturbance(t_sim)
setpoint_sim = tclab_setpoint(t_sim)
plt.figure(figsize=(10,8))
plt.subplot(3, 1, 1)
plt.plot(t_sim, setpoint_sim)
plt.title('setpoint')
plt.ylabel('deg C')
plt.grid(True)
plt.subplot(3, 1, 2)
plt.plot(t_sim, u_sim)
plt.title('heat power input')
plt.ylabel('percent of max')
plt.grid(True)
plt.subplot(3, 1, 3)
plt.plot(t_sim, d_sim)
plt.title('unmeasured disturbance')
plt.ylabel('watts')
plt.grid(True)
plt.tight_layout()
![../_images/85f5b848fb5f84edfc9c0d8b2dff540b8f2be32aaf2f60629d70b5807c961230.png](../_images/85f5b848fb5f84edfc9c0d8b2dff540b8f2be32aaf2f60629d70b5807c961230.png)
8.2.3. Pyomo Simulation#
Let’s see how well our initial guess at a control strategy will work for us.
subject to initial conditions
and prior specification of inputs \(u(t)\) and \(d(t)\).
from pyomo.environ import *
from pyomo.dae import *
m = ConcreteModel()
m.t = ContinuousSet(initialize = t_sim) # make sure the expt time grid are discretization points
m.Th = Var(m.t)
m.Ts = Var(m.t)
m.U = Var(m.t, bounds=(0, 100))
m.D = Var(m.t)
m.Thdot = DerivativeVar(m.Th, wrt = m.t)
m.Tsdot = DerivativeVar(m.Ts, wrt = m.t)
# differential equations
m.Th_ode = Constraint(m.t, rule = lambda m, t:
CpH*m.Thdot[t] == Ua*(Tamb - m.Th[t]) + Uc*(m.Ts[t] - m.Th[t]) + P*m.U[t] + m.D[t])
m.Ts_ode = Constraint(m.t, rule = lambda m, t:
CpS*m.Tsdot[t] == Uc*(m.Th[t] - m.Ts[t]))
# input specifications
m.Usim = Constraint(range(0, len(t_sim)), rule = lambda m, k: m.U[t_sim[k]] == u_sim[k])
m.Dsim = Constraint(range(0, len(t_sim)), rule = lambda m, k: m.D[t_sim[k]] == d_sim[k])
# initial conditions
m.Th[0].fix(Tamb)
m.Ts[0].fix(Tamb)
TransformationFactory('dae.finite_difference').apply_to(m, method='forward')
SolverFactory('ipopt').solve(m)
Th_sim = np.array([m.Th[t]() for t in t_sim])
Ts_sim = np.array([m.Ts[t]() for t in t_sim])
# visualization
plt.figure(figsize=(10,8))
plt.subplot(3, 1, 1)
plt.plot(t_sim, Th_sim)
plt.plot(t_sim, Ts_sim)
plt.plot(t_sim, setpoint_sim)
plt.title('temperatures')
plt.ylabel('deg C')
plt.legend(['T_heater', 'T_sensor', 'Heater Setpoint'])
plt.grid(True)
plt.subplot(3, 1, 2)
plt.plot(t_sim, np.array([m.U[t]() for t in t_sim]))
plt.title('heater power')
plt.ylabel('percent of max')
plt.grid(True)
plt.subplot(3, 1, 3)
plt.plot(t_sim, np.array([m.D[t]() for t in t_sim]))
plt.title('disturbance')
plt.ylabel('watts')
plt.grid(True)
plt.tight_layout()
WARNING: More finite elements were found in ContinuousSet 't' than the number
of finite elements specified in apply. The larger number of finite
elements will be used.
![../_images/9ff665ea97cf15626689b749169f448f7503a0b9ccd2fb9d7b92512b9a4fe921.png](../_images/9ff665ea97cf15626689b749169f448f7503a0b9ccd2fb9d7b92512b9a4fe921.png)
8.2.4. Optimal Control with Knowledge of Disturbances#
An optimal control policy minimizes the differences
subject to constraints
initial conditions
and prior knowledge of \(d(t)\).
#%%time #if you put this as the first line, jupyter will report how long the code took to execute.
from pyomo.environ import *
from pyomo.dae import *
m = ConcreteModel()
m.t = ContinuousSet(initialize = t_sim) # make sure the expt time grid are discretization points
m.Th = Var(m.t)
m.Ts = Var(m.t)
m.U = Var(m.t, bounds=(0, 100))
m.D = Var(m.t)
m.Thdot = DerivativeVar(m.Th, wrt = m.t)
m.Tsdot = DerivativeVar(m.Ts, wrt = m.t)
# differential equations
m.Th_ode = Constraint(m.t, rule = lambda m, t:
CpH*m.Thdot[t] == Ua*(Tamb - m.Th[t]) + Uc*(m.Ts[t] - m.Th[t]) + P*m.U[t] + m.D[t])
m.Ts_ode = Constraint(m.t, rule = lambda m, t:
CpS*m.Tsdot[t] == Uc*(m.Th[t] - m.Ts[t]))
# input specifications
m.Dsim = Constraint(range(0, len(t_sim)), rule = lambda m, k: m.D[t_sim[k]] == d_sim[k])
# with these two lines which provide an objective function to determine the input
m.ls_control = sum([(setpoint_sim[k] - m.Th[t_sim[k]])**2 for k in range(0, len(t_sim))])
m.obj = Objective(expr = m.ls_control, sense=minimize)
# initial conditions
m.Th[0].fix(Tamb)
m.Ts[0].fix(Tamb)
TransformationFactory('dae.finite_difference').apply_to(m, nfe=len(t_sim), method='forward')
SolverFactory('ipopt').solve(m).write()
# visualization
plt.figure(figsize=(10,8))
plt.subplot(3, 1, 1)
plt.plot(t_sim, np.array([m.Th[t]() for t in t_sim]))
plt.plot(t_sim, np.array([m.Ts[t]() for t in t_sim]))
plt.plot(t_sim, setpoint_sim)
plt.title('temperatures')
plt.ylabel('deg C')
plt.legend(['T_heater', 'T_sensor', 'Heater Setpoint'])
plt.grid(True)
plt.subplot(3, 1, 2)
plt.plot(t_sim, np.array([m.U[t]() for t in t_sim]))
plt.title('heater power')
plt.ylabel('percent of max')
plt.grid(True)
plt.subplot(3, 1, 3)
plt.plot(t_sim, np.array([m.D[t]() for t in t_sim]))
plt.title('disturbance')
plt.ylabel('watts')
plt.grid(True)
plt.tight_layout()
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Lower bound: -inf
Upper bound: inf
Number of objectives: 1
Number of constraints: 1007
Number of variables: 1210
Sense: unknown
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Message: Ipopt 3.11.1\x3a Optimal Solution Found
Termination condition: optimal
Id: 0
Error rc: 0
Time: 0.27823472023010254
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
![../_images/4d28b1d6f1de0d7d976f0b340e4d88e4502c7ec0af858da970cdc44f191d7aa0.png](../_images/4d28b1d6f1de0d7d976f0b340e4d88e4502c7ec0af858da970cdc44f191d7aa0.png)
8.2.4.1. Exercise#
The optimal control computed above requires rapid changes in power level. In process systems where control action requires movement of a valve stem position, there are often limits on how fast the manipulated variable can change. Modify the model to include differential inequalities that limit the time rate of change of control.
where \(\dot{u}_{max}\) is the maximum rate of change.
8.2.5. Estimation/Observation#
“… and now my watch begins …” ―The Night’s Watch oath, Game of Thrones
The trouble with open-loop optimal control is that we can’t anticipate or know the values of unmeasured disturbances, much less the future values of those disturbances. The best we can do is use available data and process models to estimate the process state and disturbances. The est
subject to
and initial conditions
%%time
from pyomo.environ import *
from pyomo.dae import *
m = ConcreteModel()
m.t = ContinuousSet(initialize = t_sim) # make sure the expt time grid are discretization points
m.Th = Var(m.t)
m.Ts = Var(m.t)
m.U = Var(m.t, bounds=(0, 100))
m.D = Var(m.t)
m.Thdot = DerivativeVar(m.Th, wrt = m.t)
m.Tsdot = DerivativeVar(m.Ts, wrt = m.t)
# differential equations
m.Th_ode = Constraint(m.t, rule = lambda m, t:
CpH*m.Thdot[t] == Ua*(Tamb - m.Th[t]) + Uc*(m.Ts[t] - m.Th[t]) + P*m.U[t] + m.D[t])
m.Ts_ode = Constraint(m.t, rule = lambda m, t:
CpS*m.Tsdot[t] == Uc*(m.Th[t] - m.Ts[t]))
m.Usim = Constraint(range(0, len(t_sim)), rule = lambda m, k: m.U[t_sim[k]] == u_sim[k])
#for k in range(0, len(t_sim)):
# m.U[t_sim[k]]
m.ls_observer = sum([(m.Ts[t_sim[k]] - Ts_sim[k])**2 for k in range(0, len(t_sim))])
m.obj = Objective(expr = m.ls_observer, sense=minimize)
m.Th[0].fix(Tamb)
m.Ts[0].fix(Tamb)
TransformationFactory('dae.finite_difference').apply_to(m, nfe=len(t_sim), method='forward')
SolverFactory('ipopt').solve(m).write()
# visualization
plt.figure(figsize=(10,8))
plt.subplot(3, 1, 1)
plt.plot(t_sim, np.array([m.Th[t]() for t in t_sim]))
plt.plot(t_sim, np.array([m.Ts[t]() for t in t_sim]))
plt.plot(t_sim, setpoint_sim)
plt.title('temperatures')
plt.ylabel('deg C')
plt.legend(['T_heater', 'T_sensor', 'Heater Setpoint'])
plt.grid(True)
plt.subplot(3, 1, 2)
plt.plot(t_sim, np.array([m.U[t]() for t in t_sim]))
plt.title('heater power')
plt.ylabel('percent of max')
plt.grid(True)
plt.subplot(3, 1, 3)
plt.plot(t_sim, np.array([m.D[t]() for t in t_sim]))
plt.title('disturbance')
plt.ylabel('watts')
plt.grid(True)
plt.tight_layout()
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Lower bound: -inf
Upper bound: inf
Number of objectives: 1
Number of constraints: 1007
Number of variables: 1210
Sense: unknown
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Message: Ipopt 3.11.1\x3a Optimal Solution Found
Termination condition: optimal
Id: 0
Error rc: 0
Time: 0.21533441543579102
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
Wall time: 718 ms
![../_images/a1ddbf5e18469fd1eb70b69a071d9075297b878fc9e68d698bd5e6b3ad39e108.png](../_images/a1ddbf5e18469fd1eb70b69a071d9075297b878fc9e68d698bd5e6b3ad39e108.png)
8.2.6. Coding the Observer as a Python Generator#
from pyomo.environ import *
from pyomo.dae import *
def tclab_observer(h=2):
t_hist = [-1]
u_hist = [0]
Ts_hist = [Tamb]
t_est = -1
Th_est = []
Ts_est = []
d_est = []
while True:
t_meas, u_meas, Ts_meas = yield t_est, Th_est, Ts_est, d_est
t_hist.append(t_meas)
u_hist.append(u_meas)
Ts_hist.append(Ts_meas)
t_hist = t_hist[-h:]
u_hist = u_hist[-h:]
Ts_hist = Ts_hist[-h:]
m = ConcreteModel()
m.t = ContinuousSet(initialize = t_hist)
m.Th = Var(m.t)
m.Ts = Var(m.t)
m.U = Var(m.t, bounds=(0, 100))
m.D = Var(m.t)
m.Thdot = DerivativeVar(m.Th, wrt = m.t)
m.Tsdot = DerivativeVar(m.Ts, wrt = m.t)
m.Th_ode = Constraint(m.t, rule = lambda m, t:
CpH*m.Thdot[t] == Ua*(Tamb - m.Th[t]) + Uc*(m.Ts[t] - m.Th[t]) + P*m.U[t] + m.D[t])
m.Ts_ode = Constraint(m.t, rule = lambda m, t:
CpS*m.Tsdot[t] == Uc*(m.Th[t] - m.Ts[t]))
m.Usim = Constraint(range(0, len(t_hist)), rule = lambda m, k: m.U[t_hist[k]] == u_hist[k])
m.ls_observer = sum([(m.Ts[t_hist[k]] - Ts_hist[k])**2 for k in range(0, len(t_hist))])
m.obj = Objective(expr = m.ls_observer, sense=minimize)
TransformationFactory('dae.finite_difference').apply_to(m, nfe=len(t_hist), method='forward')
SolverFactory('ipopt').solve(m)
t_est = t_hist[-1]
Th_est = m.Th[t_est]()
Ts_est = m.Ts[t_est]()
d_est = m.D[t_est]()
%%time
t_est = []
Th_est = []
Ts_est = []
d_est = []
observer = tclab_observer(5)
observer.send(None)
for k in range(0, len(t_sim)):
t, Th, Ts, d = observer.send([t_sim[k], u_sim[k], Ts_sim[k]])
t_est.append(t)
Th_est.append(Th)
Ts_est.append(Ts)
d_est.append(d)
plt.figure(figsize=(10,8))
plt.subplot(3,1,1)
plt.plot(t_est, Th_est)
plt.plot(t_sim, Th_sim)
plt.grid(True)
plt.subplot(3,1,2)
plt.plot(t_est, Ts_est)
plt.plot(t_sim, Ts_sim)
plt.grid(True)
plt.subplot(3,1,3)
plt.plot(t_est, d_est)
plt.plot(t_sim, d_sim)
plt.grid(True)
plt.tight_layout()
CPU times: user 2.53 s, sys: 1.37 s, total: 3.9 s
Wall time: 6.29 s
![../_images/9d1d12d6874aa4bc0f2617c027ff8d813487360a09d5af0bfc6fc7b427680be3.png](../_images/9d1d12d6874aa4bc0f2617c027ff8d813487360a09d5af0bfc6fc7b427680be3.png)
8.2.7. Coding the Controller as a Python Generator#
from pyomo.environ import *
from pyomo.dae import *
def tclab_control(h=100):
u = 0
while True:
t_est, Th_est, Ts_est, d_est = yield u
tf = t_est + h
m = ConcreteModel()
m.t = ContinuousSet(initialize=np.linspace(t, tf, 1 + round((tf-t)/2)))
m.Th = Var(m.t)
m.Ts = Var(m.t)
m.U = Var(m.t, bounds=(0, 100))
m.Thdot = DerivativeVar(m.Th, wrt = m.t)
m.Tsdot = DerivativeVar(m.Ts, wrt = m.t)
m.heater = Constraint(m.t, rule = lambda m, t:
CpH*m.Thdot[t] == Ua*(Tamb - m.Th[t]) + Uc*(m.Ts[t] - m.Th[t]) + P*m.U[t] + d)
m.sensor = Constraint(m.t, rule = lambda m, t:
CpS*m.Tsdot[t] == Uc*(m.Th[t] - m.Ts[t]))
m.Th[t].fix(Th_est)
m.Ts[t].fix(Ts_est)
m.obj = Objective(expr = sum([(tclab_setpoint(t) - m.Th[t])**2 for t in m.t]), sense=minimize)
TransformationFactory('dae.finite_difference').apply_to(m, nfe=len(m.t), method='backwards')
SolverFactory('ipopt').solve(m)
umpc = np.array([m.U[t]() for t in m.t])
u = umpc[1]
observer = tclab_observer(3)
observer.send(None)
controller = tclab_control()
controller.send(None)
t_mpc = []
u_mpc = []
Th_mpc = []
Ts_mpc = []
for k in range(0, len(t_sim)):
t, Th, Ts, d = observer.send([t_sim[k], u_sim[k], Ts_sim[k]])
u = controller.send([t, Th, Ts, d])
u_mpc.append(u)
Th_mpc.append(Th)
Ts_mpc.append(Ts)
plt.plot(t_sim, u_mpc)
plt.plot(t_sim, Th_mpc)
C:\Users\eliza\Anaconda3\lib\site-packages\ipykernel_launcher.py:13: DeprecationWarning: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.
del sys.path[0]
[<matplotlib.lines.Line2D at 0x1912230dba8>]
![../_images/34ba308305b823e03aef989120936093fb6108f72a1874238b7bcacad9608fdb.png](../_images/34ba308305b823e03aef989120936093fb6108f72a1874238b7bcacad9608fdb.png)
len(t_sim)
201
def tclab_control(t, h, ic, d):
tf = t + h
m = ConcreteModel()
m.t = ContinuousSet(initialize=np.linspace(t, tf, 1 + round((tf-t)/2)))
m.Th = Var(m.t)
m.Ts = Var(m.t)
m.U = Var(m.t, bounds=(0, 100))
m.Thdot = DerivativeVar(m.Th, wrt = m.t)
m.Tsdot = DerivativeVar(m.Ts, wrt = m.t)
m.heater = Constraint(m.t, rule = lambda m, t:
CpH*m.Thdot[t] == Ua*(Tamb - m.Th[t]) + Uc*(m.Ts[t] - m.Th[t]) + P*m.U[t] + d)
m.sensor = Constraint(m.t, rule = lambda m, t:
CpS*m.Tsdot[t] == Uc*(m.Th[t] - m.Ts[t]))
m.Th[t].fix(ic[0])
m.Ts[t].fix(ic[1])
m.obj = Objective(expr = sum([(tclab_setpoint(t) - m.Th[t])**2 for t in m.t]), sense=minimize)
TransformationFactory('dae.finite_difference').apply_to(m, nfe=len(m.t), method='backwards')
SolverFactory('ipopt').solve(m)
tsim = np.array([t for t in m.t])
Thsim = np.array([m.Th[t]() for t in m.t])
Tssim = np.array([m.Ts[t]() for t in m.t])
Usim = np.array([m.U[t]() for t in m.t])
return [tsim, Usim, Thsim, Tssim]
[tsim, Usim, Thsim, Tssim] = tclab_control(100, 800, [30, 30], .5)
# visualization
plt.subplot(2,1,1)
plt.plot(tsim, Thsim, tsim, Tssim)
plt.title('temperatures')
plt.ylabel('deg C')
plt.legend(['T_heater', 'T_sensor'])
plt.grid(True)
plt.subplot(2,1,2)
plt.plot(tsim, Usim)
plt.title('heater input')
plt.ylabel('percent of max')
plt.grid(True)
plt.tight_layout()
![../_images/2aba53862f217e9521f59bfa7ee75931473f188047853aad96a4e1b15692689b.png](../_images/2aba53862f217e9521f59bfa7ee75931473f188047853aad96a4e1b15692689b.png)
from tclab import setup, clock, Historian, Plotter
TCLab = setup(connected=False, speedup=20)
tf = 800 # run time
# create a controller instance
controller = tclab_control(100)
controller.send(None)
# create an model estimator
observer = tclab_observer(5)
observer.send(None)
# execute the event loop
tf = 800
with TCLab() as lab:
h = Historian([('T1', lambda: lab.T1), ('Q1', lab.Q1),
('Th', lambda: Th), ('Ts', lambda: Ts)])
p = Plotter(h, tf)
U1 = 0
for t in clock(tf, 5): # allow time for more calculations
T1 = lab.T1 # measure the sensor temperature
t, Th, Ts, d = observer.send([t, U1, T1]) # estimate the heater temperature # get setpoint
U1 = controller.send([t, Th, Ts, d]) # compute control action
lab.U1 = U1 # set manipulated variable
p.update(t) # log data
![../_images/840adad95ff8dd1697f3e2779aba558f88b1a07c5fc453e0bcdef2654324eac9.png](../_images/840adad95ff8dd1697f3e2779aba558f88b1a07c5fc453e0bcdef2654324eac9.png)
TCLab Model disconnected successfully.
![../_images/840adad95ff8dd1697f3e2779aba558f88b1a07c5fc453e0bcdef2654324eac9.png](../_images/840adad95ff8dd1697f3e2779aba558f88b1a07c5fc453e0bcdef2654324eac9.png)