TCLab G - Nonlinear MPC
Main.TCLabG History
Show minor edits - Show changes to output
Changed line 220 from:
plt.plot(tm[0:i],Tsp1[0:i],'k-',LineWidth=2,label=r'$T_1 SP$')
to:
plt.plot(tm[0:i],Tsp1[0:i],'k-',lw=2,label=r'$T_1 SP$')
Changed line 226 from:
plt.plot(tm[0:i],Tsp2[0:i],'g-',LineWidth=2,label=r'$T_2 SP$')
to:
plt.plot(tm[0:i],Tsp2[0:i],'g-',lw=2,label=r'$T_2 SP$')
Changed lines 231-232 from:
plt.plot(tm[0:i],Q1s[0:i],'r-',LineWidth=3,label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',LineWidth=3,label=r'$Q_2$')
plt.plot(tm[0:i],Q2s[0:i],'b:',
to:
plt.plot(tm[0:i],Q1s[0:i],'r-',lw=3,label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',lw=3,label=r'$Q_2$')
plt.plot(tm[0:i],Q2s[0:i],'b:',lw=3,label=r'$Q_2$')
Changed line 563 from:
plt.plot(tm[0:i],T1sp[0:i],'k-',LineWidth=2,label=r'$T_1 SP$')
to:
plt.plot(tm[0:i],T1sp[0:i],'k-',lw=2,label=r'$T_1 SP$')
Changed line 569 from:
plt.plot(tm[0:i],T2sp[0:i],'g-',LineWidth=2,label=r'$T_2 SP$')
to:
plt.plot(tm[0:i],T2sp[0:i],'g-',lw=2,label=r'$T_2 SP$')
Changed lines 574-575 from:
plt.plot(tm[0:i],Q1s[0:i],'r-',LineWidth=3,label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',LineWidth=3,label=r'$Q_2$')
plt.plot(tm[0:i],Q2s[0:i],'b:',
to:
plt.plot(tm[0:i],Q1s[0:i],'r-',lw=3,label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',lw=3,label=r'$Q_2$')
plt.plot(tm[0:i],Q2s[0:i],'b:',lw=3,label=r'$Q_2$')
Added lines 619-621:
[[https://github.com/APMonitor/dynopt/blob/master/DynamicOptimization.ipynb|Virtual TCLab on Google Colab]]
-> [[https://colab.research.google.com/github/APMonitor/dynopt/blob/master/TCLabA.ipynb|Lab A]] | [[https://colab.research.google.com/github/APMonitor/dynopt/blob/master/TCLabB.ipynb|Lab B]] | [[https://colab.research.google.com/github/APMonitor/dynopt/blob/master/TCLabC.ipynb|Lab C]] | [[https://colab.research.google.com/github/APMonitor/dynopt/blob/master/TCLabD.ipynb|Lab D]] | [[https://colab.research.google.com/github/APMonitor/dynopt/blob/master/TCLabE.ipynb|Lab E]] | [[https://colab.research.google.com/github/APMonitor/dynopt/blob/master/TCLabF.ipynb|Lab F]] | [[https://colab.research.google.com/github/APMonitor/dynopt/blob/master/TCLabG.ipynb|Lab G]] | [[https://colab.research.google.com/github/APMonitor/dynopt/blob/master/TCLabH.ipynb|Lab H]]
Changed lines 56-60 from:
to:
Tsp1[3:] = 40.0
Tsp2[40:] = 30.0
Tsp1[80:] = 32.0
Tsp2[120:] = 35.0
Tsp1[160:] = 45.0
Tsp2[40:] = 30.0
Tsp1[80:] = 32.0
Tsp2[120:] = 35.0
Tsp1[160:] = 45.0
Changed line 220 from:
plt.plot(tm[0:i],Tsp1[0:i],'k-',LineWidth=2,label=r'$T_1 Setpoint$')
to:
plt.plot(tm[0:i],Tsp1[0:i],'k-',LineWidth=2,label=r'$T_1 SP$')
Changed line 226 from:
plt.plot(tm[0:i],Tsp2[0:i],'g-',LineWidth=2,label=r'$T_2 Setpoint$')
to:
plt.plot(tm[0:i],Tsp2[0:i],'g-',LineWidth=2,label=r'$T_2 SP$')
Changed line 563 from:
plt.plot(tm[0:i],T1sp[0:i],'k-',LineWidth=2,label=r'$T_1 Setpoint$')
to:
plt.plot(tm[0:i],T1sp[0:i],'k-',LineWidth=2,label=r'$T_1 SP$')
Changed line 569 from:
plt.plot(tm[0:i],T2sp[0:i],'g-',LineWidth=2,label=r'$T_2 Setpoint$')
to:
plt.plot(tm[0:i],T2sp[0:i],'g-',LineWidth=2,label=r'$T_2 SP$')
Changed line 272 from:
* [[Attach:tclab_ss_data1.txt|Steady state data set 1 (low power heaters)]]
to:
* [[Attach:tclab_ss_data3.txt|Steady state data set 1 (low power heaters)]]
Changed line 292 from:
url = 'https://apmonitor.com/do/uploads/Main/tclab_ss_data1.txt'
to:
url = 'https://apmonitor.com/do/uploads/Main/tclab_ss_data3.txt'
Added lines 269-274:
Training data is generated from the script in [[Main/TCLabB|TCLab B]] or use one of the data files below.
* [[Attach:tclab_ss_data1.txt|Steady state data set 1 (low power heaters)]]
* [[Attach:tclab_ss_data2.txt|Steady state data set 2 (high power heaters)]]
Changed line 292 from:
url = 'https://apmonitor.com/do/uploads/Main/tclab_ss_data2.txt'
to:
url = 'https://apmonitor.com/do/uploads/Main/tclab_ss_data1.txt'
Changed line 67 from:
# Initialize Model as Estimator
to:
# Initialize Model
Added lines 252-588:
# Make sure serial connection still closes when there's an error
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Error: Shutting down')
a.close()
raise
(:sourceend:)
(:divend:)
----
%width=550px%Attach:tclab_neural_network_mpc.png
(:toggle hide gekko_labGNN button show="Lab G: Python TCLab MIMO MPC with Neural Network Model":)
(:div id=gekko_labGNN:)
(:source lang=python:)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from gekko import GEKKO
import tclab
import time
# -------------------------------------
# connect to Arduino
# -------------------------------------
a = tclab.TCLab()
# -------------------------------------
# import data
# -------------------------------------
url = 'https://apmonitor.com/do/uploads/Main/tclab_ss_data2.txt'
data = pd.read_csv(url)
# -------------------------------------
# scale data
# -------------------------------------
s = MinMaxScaler(feature_range=(0,1))
sc_train = s.fit_transform(data)
# partition into inputs and outputs
xs = sc_train[:,0:2] # 2 heaters
ys = sc_train[:,2:4] # 2 temperatures
# -------------------------------------
# build neural network
# -------------------------------------
nin = 2 # inputs
n1 = 2 # hidden layer 1 (linear)
n2 = 2 # hidden layer 2 (nonlinear)
n3 = 2 # hidden layer 3 (linear)
nout = 2 # outputs
# Initialize gekko models
train = GEKKO()
mpc = GEKKO(remote=False)
model = [train,mpc]
for m in model:
# use APOPT solver
m.options.SOLVER = 1
# input(s)
if m==train:
# parameter for training
m.inpt = [m.Param() for i in range(nin)]
else:
# variable for MPC
m.inpt = [m.Var() for i in range(nin)]
# layer 1 (linear)
m.w1 = m.Array(m.FV, (nout,nin,n1))
m.l1 = [[m.Intermediate(sum([m.w1[k,j,i]*m.inpt[j] \
for j in range(nin)])) for i in range(n1)] \
for k in range(nout)]
# layer 2 (tanh)
m.w2 = m.Array(m.FV, (nout,n1,n2))
m.l2 = [[m.Intermediate(sum([m.tanh(m.w2[k,j,i]*m.l1[k][j]) \
for j in range(n1)])) for i in range(n2)] \
for k in range(nout)]
# layer 3 (linear)
m.w3 = m.Array(m.FV, (nout,n2,n3))
m.l3 = [[m.Intermediate(sum([m.w3[k,j,i]*m.l2[k][j] \
for j in range(n2)])) for i in range(n3)] \
for k in range(nout)]
# outputs
m.outpt = [m.CV() for i in range(nout)]
m.Equations([m.outpt[k]==sum([m.l3[k][i] for i in range(n3)]) \
for k in range(nout)])
# flatten matrices
m.w1 = m.w1.flatten()
m.w2 = m.w2.flatten()
m.w3 = m.w3.flatten()
# -------------------------------------
# fit parameter weights
# -------------------------------------
m = train
for i in range(nin):
m.inpt[i].value=xs[:,i]
for i in range(nout):
m.outpt[i].value = ys[:,i]
m.outpt[i].FSTATUS = 1
for i in range(len(m.w1)):
m.w1[i].FSTATUS=1
m.w1[i].STATUS=1
m.w1[i].MEAS=1.0
for i in range(len(m.w2)):
m.w2[i].STATUS=1
m.w2[i].FSTATUS=1
m.w2[i].MEAS=0.5
for i in range(len(m.w3)):
m.w3[i].FSTATUS=1
m.w3[i].STATUS=1
m.w3[i].MEAS=1.0
m.options.IMODE = 2
m.options.EV_TYPE = 2
# solve for weights to minimize loss (objective)
m.solve(disp=True)
# -------------------------------------
# Create Model Predictive Controller
# -------------------------------------
m = mpc
# 60 second time horizon, steps of 3 sec
m.time = np.linspace(0,60,21)
# load neural network parameters
for i in range(len(m.w1)):
m.w1[i].MEAS=train.w1[i].NEWVAL
m.w1[i].FSTATUS = 1
for i in range(len(m.w2)):
m.w2[i].MEAS=train.w2[i].NEWVAL
m.w2[i].FSTATUS = 1
for i in range(len(m.w3)):
m.w3[i].MEAS=train.w3[i].NEWVAL
m.w3[i].FSTATUS = 1
# MVs and CVs
Q1 = m.MV(value=0)
Q2 = m.MV(value=0)
TC1 = m.CV(value=a.T1)
TC2 = m.CV(value=a.T2)
# scaled inputs to neural network
m.Equation(m.inpt[0] == Q1 * s.scale_[0] + s.min_[0])
m.Equation(m.inpt[1] == Q2 * s.scale_[1] + s.min_[1])
# define Temperature output
Q0 = 0 # initial heater
T0 = 23 # ambient temperature
# scaled steady state ouput
T1_ss = m.Var(value=T0)
T2_ss = m.Var(value=T0)
m.Equation(T1_ss == (m.outpt[0]-s.min_[2])/s.scale_[2])
m.Equation(T2_ss == (m.outpt[1]-s.min_[3])/s.scale_[3])
# time constants
tauA = m.Param(value=80)
tauB = m.Param(value=20)
TH1 = m.Var(a.T1)
TH2 = m.Var(a.T2)
# additional model equation for dynamics
m.Equation(tauA*TH1.dt()==-TH1+T1_ss)
m.Equation(tauA*TH2.dt()==-TH2+T2_ss)
m.Equation(tauB*TC1.dt()==-TC1+TH1)
m.Equation(tauB*TC2.dt()==-TC2+TH2)
# Manipulated variable tuning
Q1.STATUS = 1 # use to control temperature
Q1.FSTATUS = 0 # no feedback measurement
Q1.LOWER = 0.0
Q1.UPPER = 100.0
Q1.DMAX = 40.0
Q1.COST = 0.0
Q1.DCOST = 0.0
Q2.STATUS = 1 # use to control temperature
Q2.FSTATUS = 0 # no feedback measurement
Q2.LOWER = 0.0
Q2.UPPER = 100.0
Q2.DMAX = 40.0
Q2.COST = 0.0
Q2.DCOST = 0.0
# Controlled variable tuning
TC1.STATUS = 1 # minimize error with setpoint range
TC1.FSTATUS = 1 # receive measurement
TC1.TR_INIT = 1 # reference trajectory
TC1.TAU = 10 # time constant for response
TC2.STATUS = 1 # minimize error with setpoint range
TC2.FSTATUS = 1 # receive measurement
TC2.TR_INIT = 1 # reference trajectory
TC2.TAU = 10 # time constant for response
# Global Options
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # Objective type
m.options.NODES = 3 # Collocation nodes
m.options.SOLVER = 3 # 1=APOPT, 3=IPOPT
# -------------------------------------
# Initialize model and data storage
# -------------------------------------
# Get Version
print(a.version)
# Turn LED on
print('LED On')
a.LED(100)
# Run time in minutes
run_time = 10.0
# Number of cycles with 3 second intervals
loops = int(20.0*run_time)
tm = np.zeros(loops)
# Temperature (K)
T1 = np.ones(loops) * a.T1 # temperature (degC)
T1sp = np.ones(loops) * 35.0 # set point (degC)
T2 = np.ones(loops) * a.T2 # temperature (degC)
T2sp = np.ones(loops) * 23.0 # set point (degC)
# Set point changes
T1sp[3:] = 40.0
T2sp[40:] = 30.0
T1sp[80:] = 32.0
T2sp[120:] = 35.0
T1sp[160:] = 45.0
# heater values
Q1s = np.ones(loops) * 0.0
Q2s = np.ones(loops) * 0.0
# Create plot
plt.figure()
plt.ion()
plt.show()
# Main Loop
start_time = time.time()
prev_time = start_time
try:
for i in range(1,loops):
# Sleep time
sleep_max = 3.0
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep)
else:
time.sleep(0.01)
# Record time and change in time
t = time.time()
dt = t - prev_time
prev_time = t
tm[i] = t - start_time
# Read temperatures in Kelvin
T1[i] = a.T1
T2[i] = a.T2
###############################
### MPC CONTROLLER ###
###############################
TC1.MEAS = T1[i]
TC2.MEAS = T2[i]
# input setpoint with deadband +/- DT
DT = 0.1
TC1.SPHI = T1sp[i] + DT
TC1.SPLO = T1sp[i] - DT
TC2.SPHI = T2sp[i] + DT
TC2.SPLO = T2sp[i] - DT
# solve MPC
m.solve(disp=False)
# test for successful solution
if (m.options.APPSTATUS==1):
# retrieve the first Q value
Q1s[i] = Q1.NEWVAL
Q2s[i] = Q2.NEWVAL
else:
# not successful, set heater to zero
Q1s[i] = 0
Q2s[i] = 0
# Write output (0-100)
a.Q1(Q1s[i])
a.Q2(Q2s[i])
# Plot
plt.clf()
ax=plt.subplot(3,1,1)
ax.grid()
plt.plot(tm[0:i],T1[0:i],'ro',MarkerSize=3,label=r'$T_1$')
plt.plot(tm[0:i],T1sp[0:i],'k-',LineWidth=2,label=r'$T_1 Setpoint$')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')
ax=plt.subplot(3,1,2)
ax.grid()
plt.plot(tm[0:i],T2[0:i],'ro',MarkerSize=3,label=r'$T_2$')
plt.plot(tm[0:i],T2sp[0:i],'g-',LineWidth=2,label=r'$T_2 Setpoint$')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')
ax=plt.subplot(3,1,3)
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',LineWidth=3,label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',LineWidth=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.draw()
plt.pause(0.05)
# Turn off heaters
a.Q1(0)
a.Q2(0)
print('Shutting down')
# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Shutting down')
a.close()
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Error: Shutting down')
a.close()
raise
(:sourceend:)
(:divend:)
----
%width=550px%Attach:tclab_neural_network_mpc.png
(:toggle hide gekko_labGNN button show="Lab G: Python TCLab MIMO MPC with Neural Network Model":)
(:div id=gekko_labGNN:)
(:source lang=python:)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from gekko import GEKKO
import tclab
import time
# -------------------------------------
# connect to Arduino
# -------------------------------------
a = tclab.TCLab()
# -------------------------------------
# import data
# -------------------------------------
url = 'https://apmonitor.com/do/uploads/Main/tclab_ss_data2.txt'
data = pd.read_csv(url)
# -------------------------------------
# scale data
# -------------------------------------
s = MinMaxScaler(feature_range=(0,1))
sc_train = s.fit_transform(data)
# partition into inputs and outputs
xs = sc_train[:,0:2] # 2 heaters
ys = sc_train[:,2:4] # 2 temperatures
# -------------------------------------
# build neural network
# -------------------------------------
nin = 2 # inputs
n1 = 2 # hidden layer 1 (linear)
n2 = 2 # hidden layer 2 (nonlinear)
n3 = 2 # hidden layer 3 (linear)
nout = 2 # outputs
# Initialize gekko models
train = GEKKO()
mpc = GEKKO(remote=False)
model = [train,mpc]
for m in model:
# use APOPT solver
m.options.SOLVER = 1
# input(s)
if m==train:
# parameter for training
m.inpt = [m.Param() for i in range(nin)]
else:
# variable for MPC
m.inpt = [m.Var() for i in range(nin)]
# layer 1 (linear)
m.w1 = m.Array(m.FV, (nout,nin,n1))
m.l1 = [[m.Intermediate(sum([m.w1[k,j,i]*m.inpt[j] \
for j in range(nin)])) for i in range(n1)] \
for k in range(nout)]
# layer 2 (tanh)
m.w2 = m.Array(m.FV, (nout,n1,n2))
m.l2 = [[m.Intermediate(sum([m.tanh(m.w2[k,j,i]*m.l1[k][j]) \
for j in range(n1)])) for i in range(n2)] \
for k in range(nout)]
# layer 3 (linear)
m.w3 = m.Array(m.FV, (nout,n2,n3))
m.l3 = [[m.Intermediate(sum([m.w3[k,j,i]*m.l2[k][j] \
for j in range(n2)])) for i in range(n3)] \
for k in range(nout)]
# outputs
m.outpt = [m.CV() for i in range(nout)]
m.Equations([m.outpt[k]==sum([m.l3[k][i] for i in range(n3)]) \
for k in range(nout)])
# flatten matrices
m.w1 = m.w1.flatten()
m.w2 = m.w2.flatten()
m.w3 = m.w3.flatten()
# -------------------------------------
# fit parameter weights
# -------------------------------------
m = train
for i in range(nin):
m.inpt[i].value=xs[:,i]
for i in range(nout):
m.outpt[i].value = ys[:,i]
m.outpt[i].FSTATUS = 1
for i in range(len(m.w1)):
m.w1[i].FSTATUS=1
m.w1[i].STATUS=1
m.w1[i].MEAS=1.0
for i in range(len(m.w2)):
m.w2[i].STATUS=1
m.w2[i].FSTATUS=1
m.w2[i].MEAS=0.5
for i in range(len(m.w3)):
m.w3[i].FSTATUS=1
m.w3[i].STATUS=1
m.w3[i].MEAS=1.0
m.options.IMODE = 2
m.options.EV_TYPE = 2
# solve for weights to minimize loss (objective)
m.solve(disp=True)
# -------------------------------------
# Create Model Predictive Controller
# -------------------------------------
m = mpc
# 60 second time horizon, steps of 3 sec
m.time = np.linspace(0,60,21)
# load neural network parameters
for i in range(len(m.w1)):
m.w1[i].MEAS=train.w1[i].NEWVAL
m.w1[i].FSTATUS = 1
for i in range(len(m.w2)):
m.w2[i].MEAS=train.w2[i].NEWVAL
m.w2[i].FSTATUS = 1
for i in range(len(m.w3)):
m.w3[i].MEAS=train.w3[i].NEWVAL
m.w3[i].FSTATUS = 1
# MVs and CVs
Q1 = m.MV(value=0)
Q2 = m.MV(value=0)
TC1 = m.CV(value=a.T1)
TC2 = m.CV(value=a.T2)
# scaled inputs to neural network
m.Equation(m.inpt[0] == Q1 * s.scale_[0] + s.min_[0])
m.Equation(m.inpt[1] == Q2 * s.scale_[1] + s.min_[1])
# define Temperature output
Q0 = 0 # initial heater
T0 = 23 # ambient temperature
# scaled steady state ouput
T1_ss = m.Var(value=T0)
T2_ss = m.Var(value=T0)
m.Equation(T1_ss == (m.outpt[0]-s.min_[2])/s.scale_[2])
m.Equation(T2_ss == (m.outpt[1]-s.min_[3])/s.scale_[3])
# time constants
tauA = m.Param(value=80)
tauB = m.Param(value=20)
TH1 = m.Var(a.T1)
TH2 = m.Var(a.T2)
# additional model equation for dynamics
m.Equation(tauA*TH1.dt()==-TH1+T1_ss)
m.Equation(tauA*TH2.dt()==-TH2+T2_ss)
m.Equation(tauB*TC1.dt()==-TC1+TH1)
m.Equation(tauB*TC2.dt()==-TC2+TH2)
# Manipulated variable tuning
Q1.STATUS = 1 # use to control temperature
Q1.FSTATUS = 0 # no feedback measurement
Q1.LOWER = 0.0
Q1.UPPER = 100.0
Q1.DMAX = 40.0
Q1.COST = 0.0
Q1.DCOST = 0.0
Q2.STATUS = 1 # use to control temperature
Q2.FSTATUS = 0 # no feedback measurement
Q2.LOWER = 0.0
Q2.UPPER = 100.0
Q2.DMAX = 40.0
Q2.COST = 0.0
Q2.DCOST = 0.0
# Controlled variable tuning
TC1.STATUS = 1 # minimize error with setpoint range
TC1.FSTATUS = 1 # receive measurement
TC1.TR_INIT = 1 # reference trajectory
TC1.TAU = 10 # time constant for response
TC2.STATUS = 1 # minimize error with setpoint range
TC2.FSTATUS = 1 # receive measurement
TC2.TR_INIT = 1 # reference trajectory
TC2.TAU = 10 # time constant for response
# Global Options
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # Objective type
m.options.NODES = 3 # Collocation nodes
m.options.SOLVER = 3 # 1=APOPT, 3=IPOPT
# -------------------------------------
# Initialize model and data storage
# -------------------------------------
# Get Version
print(a.version)
# Turn LED on
print('LED On')
a.LED(100)
# Run time in minutes
run_time = 10.0
# Number of cycles with 3 second intervals
loops = int(20.0*run_time)
tm = np.zeros(loops)
# Temperature (K)
T1 = np.ones(loops) * a.T1 # temperature (degC)
T1sp = np.ones(loops) * 35.0 # set point (degC)
T2 = np.ones(loops) * a.T2 # temperature (degC)
T2sp = np.ones(loops) * 23.0 # set point (degC)
# Set point changes
T1sp[3:] = 40.0
T2sp[40:] = 30.0
T1sp[80:] = 32.0
T2sp[120:] = 35.0
T1sp[160:] = 45.0
# heater values
Q1s = np.ones(loops) * 0.0
Q2s = np.ones(loops) * 0.0
# Create plot
plt.figure()
plt.ion()
plt.show()
# Main Loop
start_time = time.time()
prev_time = start_time
try:
for i in range(1,loops):
# Sleep time
sleep_max = 3.0
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep)
else:
time.sleep(0.01)
# Record time and change in time
t = time.time()
dt = t - prev_time
prev_time = t
tm[i] = t - start_time
# Read temperatures in Kelvin
T1[i] = a.T1
T2[i] = a.T2
###############################
### MPC CONTROLLER ###
###############################
TC1.MEAS = T1[i]
TC2.MEAS = T2[i]
# input setpoint with deadband +/- DT
DT = 0.1
TC1.SPHI = T1sp[i] + DT
TC1.SPLO = T1sp[i] - DT
TC2.SPHI = T2sp[i] + DT
TC2.SPLO = T2sp[i] - DT
# solve MPC
m.solve(disp=False)
# test for successful solution
if (m.options.APPSTATUS==1):
# retrieve the first Q value
Q1s[i] = Q1.NEWVAL
Q2s[i] = Q2.NEWVAL
else:
# not successful, set heater to zero
Q1s[i] = 0
Q2s[i] = 0
# Write output (0-100)
a.Q1(Q1s[i])
a.Q2(Q2s[i])
# Plot
plt.clf()
ax=plt.subplot(3,1,1)
ax.grid()
plt.plot(tm[0:i],T1[0:i],'ro',MarkerSize=3,label=r'$T_1$')
plt.plot(tm[0:i],T1sp[0:i],'k-',LineWidth=2,label=r'$T_1 Setpoint$')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')
ax=plt.subplot(3,1,2)
ax.grid()
plt.plot(tm[0:i],T2[0:i],'ro',MarkerSize=3,label=r'$T_2$')
plt.plot(tm[0:i],T2sp[0:i],'g-',LineWidth=2,label=r'$T_2 Setpoint$')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')
ax=plt.subplot(3,1,3)
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',LineWidth=3,label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',LineWidth=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.draw()
plt.pause(0.05)
# Turn off heaters
a.Q1(0)
a.Q2(0)
print('Shutting down')
# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Shutting down')
a.close()
Added line 26:
import tclab
Deleted lines 29-31:
# get gekko package with:
# pip install gekko
Changed lines 31-34 from:
# pip install tclab
from tclab import TCLab
to:
Changed lines 33-51 from:
a = TCLab()
#Final time
tf = 10 # min
# number of data points (every 3 seconds)
n = tf * 20 + 1
# Percent Heater(0-100%)
Q1s = np.zeros(n)
Q2s = np.zeros(n)
# Temperatures (degC)
T1m = a.T1 * np.ones(n)
T2m = a.T2 * np.ones(n)
# Temperature setpoints
T1sp = T1m[0] * np.ones(n)
T2sp = T2m[0] * np.ones(n)
# Heater set point steps about every 150 sec
#
tf = 10 # min
n = tf * 20 + 1
# Percent Heater
T2m = a.T2 * np.ones
T1sp = T1m[0] * np.ones
# Heater set point steps about every 150 sec
to:
a = tclab.TCLab()
# Get Version
print(a.version)
# Turn LED on
print('LED On')
a.LED(100)
# Run time in minutes
run_time = 10.0
# Number of cycles with 3 second intervals
loops = int(20.0*run_time)
tm = np.zeros(loops)
# Temperature (K)
T1 = np.ones(loops) * a.T1 # temperature (degC)
Tsp1 = np.ones(loops) * 35.0 # set point (degC)
T2 = np.ones(loops) * a.T2 # temperature (degC)
Tsp2 = np.ones(loops) * 23.0 # set point (degC)
# Set point changes
# Get Version
print(a.version)
# Turn LED on
print('LED On')
a.LED(100)
# Run time in minutes
run_time = 10.0
# Number of cycles with 3 second intervals
loops = int(20.0*run_time)
tm = np.zeros(loops)
# Temperature (K)
T1 = np.ones(loops) * a.T1 # temperature (degC)
Tsp1 = np.ones(loops) * 35.0 # set point (degC)
T2 = np.ones(loops) * a.T2 # temperature (degC)
Tsp2 = np.ones(loops) * 23.0 # set point (degC)
# Set point changes
Added lines 62-65:
# heater values
Q1s = np.ones(loops) * 0.0
Q2s = np.ones(loops) * 0.0
Q1s = np.ones(loops) * 0.0
Q2s = np.ones(loops) * 0.0
Changed line 69 from:
# use remote=True for MacOS
to:
# remote=True for MacOS
Changed lines 72-90 from:
# Control horizon, non-uniform time steps
m.time =[0,3,6,10,14,18,22,27,32,38,45,55,65, \
75,90,110,130,150]
# Parameters from Estimation
K1 = m.FV(value=0.607)
K2 = m.FV(value=0.293)
K3 = m.FV(value=0.24)
tau12 = m.FV(value=192)
tau3 = m.FV(value=15)
# don't update parameters with optimizer
K1.STATUS = 0
K2.STATUS = 0
K3.STATUS = 0
tau12.STATUS = 0
tau3.STATUS = 0
# Measured inputs
m.time =
# Parameters from Estimation
K1 =
tau3 = m.FV(value=15)
# don't update parameters with optimizer
K1.STATUS = 0
K2.STATUS = 0
K3.STATUS = 0
tau12.STATUS = 0
tau3.STATUS = 0
# Measured inputs
to:
# 60 second time horizon, steps of 3 sec
m.time = np.linspace(0,60,21)
# Parameters
U = m.FV(value=10,name='u')
tau = m.FV(value=5,name='tau')
alpha1 = m.FV(value=0.01) # W / % heater
alpha2 = m.FV(value=0.0075) # W / % heater
Kp = m.Param(value=0.5)
tau = m.Param(value=50.0)
zeta = m.Param(value=1.5)
# Manipulated variables
m.time = np.linspace(0,60,21)
# Parameters
U = m.FV(value=10,name='u')
tau = m.FV(value=5,name='tau')
alpha1 = m.FV(value=0.01) # W / % heater
alpha2 = m.FV(value=0.0075) # W / % heater
Kp = m.Param(value=0.5)
tau = m.Param(value=50.0)
zeta = m.Param(value=1.5)
# Manipulated variables
Changed lines 87-91 from:
Q1.STATUS = 1 # manipulated
Q1.FSTATUS = 0 # measured
Q1.DMAX = 20.0
Q1.DCOST = 0.1
Q1.UPPER = 100.0
Q1.FSTATUS = 0 # measured
Q1.
Q1.UPPER = 100.0
to:
Q1.STATUS = 1 # use to control temperature
Q1.FSTATUS = 0 # no feedback measurement
Q1.FSTATUS = 0 # no feedback measurement
Changed lines 90-94 from:
to:
Q1.UPPER = 100.0
Q1.DMAX = 40.0
Q1.COST = 0.0
Q1.DCOST = 0.0
Q1.DMAX = 40.0
Q1.COST = 0.0
Q1.DCOST = 0.0
Changed lines 96-100 from:
Q2.STATUS = 1 # manipulated
Q2.FSTATUS = 0 # measured
Q2.DMAX = 30.0
Q2.DCOST = 0.1
Q2.UPPER = 100.0
Q2.FSTATUS = 0 # measured
Q2.
Q2.UPPER = 100.0
to:
Q2.STATUS = 1 # use to control temperature
Q2.FSTATUS = 0 # no feedback measurement
Q2.FSTATUS = 0 # no feedback measurement
Changed lines 99-106 from:
# State variables
TH1
TH2
# Measurements for model alignment
TC1 = m.CV(value=
TC1.STATUS = 1 # minimize
to:
Q2.UPPER = 100.0
Q2.DMAX = 40.0
Q2.COST = 0.0
Q2.DCOST = 0.0
# Controlled variable
TC1 = m.CV(value=T1[0])
TC1.STATUS = 1 # minimize error with setpoint range
Q2.DMAX = 40.0
Q2.COST = 0.0
Q2.DCOST = 0.0
# Controlled variable
TC1 = m.CV(value=T1[0])
TC1.STATUS = 1 # minimize error with setpoint range
Changed lines 108-112 from:
TC1.TAU = 10 # response speed (time constant)
TC1.TR_INIT = 2 # reference trajectory
TC2 = m.CV(value=T2m[0])
TC2.STATUS = 1 # minimizeerror
TC1.TR_INIT = 2
TC2 =
TC2.STATUS = 1 # minimize
to:
TC1.TR_INIT = 1 # reference trajectory
TC1.TAU = 10 # time constant for response
# Controlled variable
TC2 = m.CV(value=T2[0])
TC2.STATUS = 1 # minimize error with setpoint range
TC1.TAU = 10 # time constant for response
# Controlled variable
TC2 = m.CV(value=T2[0])
TC2.STATUS = 1 # minimize error with setpoint range
Changed lines 115-119 from:
TC2.TAU = 10 # response speed (time constant)
TC2.TR_INIT = 2 # reference trajectory
Ta = m.Param(value=23.0) # degC
TC2.TR_INIT = 2
Ta =
to:
TC2.TR_INIT = 1 # reference trajectory
TC2.TAU = 10 # time constant for response
# State variables
TH1 = m.SV(value=T1[0])
TH2 = m.SV(value=T2[0])
Ta = m.Param(value=23.0+273.15) # K
mass = m.Param(value=4.0/1000.0) # kg
Cp = m.Param(value=0.5*1000.0) # J/kg-K
A = m.Param(value=10.0/100.0**2) # Area not between heaters in m^2
As = m.Param(value=2.0/100.0**2) # Area between heaters in m^2
eps = m.Param(value=0.9) # Emissivity
sigma = m.Const(5.67e-8) # Stefan-Boltzmann
# Heater temperatures
T1i = m.Intermediate(TH1+273.15)
T2i = m.Intermediate(TH2+273.15)
TC2.TAU = 10 # time constant for response
# State variables
TH1 = m.SV(value=T1[0])
TH2 = m.SV(value=T2[0])
Ta = m.Param(value=23.0+273.15) # K
mass = m.Param(value=4.0/1000.0) # kg
Cp = m.Param(value=0.5*1000.0) # J/kg-K
A = m.Param(value=10.0/100.0**2) # Area not between heaters in m^2
As = m.Param(value=2.0/100.0**2) # Area between heaters in m^2
eps = m.Param(value=0.9) # Emissivity
sigma = m.Const(5.67e-8) # Stefan-Boltzmann
# Heater temperatures
T1i = m.Intermediate(TH1+273.15)
T2i = m.Intermediate(TH2+273.15)
Changed lines 135-142 from:
# Empirical correlations
m.
m.Equation(
m.Equation(tau3 * TC2.dt()
to:
Q_C12 = m.Intermediate(U*As*(T2i-T1i)) # Convective
Q_R12 = m.Intermediate(eps*sigma*As*(T2i**4-T1i**4)) # Radiative
# Semi-fundamental correlations (energy balances)
m.Equation(TH1.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T1i) \
+ eps * sigma * A * (Ta**4 - T1i**4) \
+ Q_C12 + Q_R12 \
+ alpha1*Q1))
m.Equation(TH2.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T2i) \
+ eps * sigma * A * (Ta**4 - T2i**4) \
- Q_C12 - Q_R12 \
+ alpha2*Q2))
# Empirical correlations (lag equations to emulate conduction)
m.Equation(tau * TC1.dt() == -TC1 + TH1)
m.Equation(tau * TC2.dt() == -TC2 + TH2)
Q_R12 = m.Intermediate(eps*sigma*As*(T2i**4-T1i**4)) # Radiative
# Semi-fundamental correlations (energy balances)
m.Equation(TH1.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T1i) \
+ eps * sigma * A * (Ta**4 - T1i**4) \
+ Q_C12 + Q_R12 \
+ alpha1*Q1))
m.Equation(TH2.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T2i) \
+ eps * sigma * A * (Ta**4 - T2i**4) \
- Q_C12 - Q_R12 \
+ alpha2*Q2))
# Empirical correlations (lag equations to emulate conduction)
m.Equation(tau * TC1.dt() == -TC1 + TH1)
m.Equation(tau * TC2.dt() == -TC2 + TH2)
Changed lines 154-155 from:
m.options.IMODE = 6 # MHE
m.options.EV_TYPE = 1 # Objective type
m.options.
to:
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # Objective type
m.options.CV_TYPE = 1 # Objective type
Changed lines 157-158 from:
m.options.SOLVER = 3 # IPOPT
m.options.COLDSTART= 1 # COLDSTART on first cycle
m.options.COLDSTART
to:
m.options.SOLVER = 3 # 1=APOPT, 3=IPOPT
Added line 159:
Changed line 161 from:
plt.figure(figsize=(10,7))
to:
plt.figure()
Deleted lines 167-168:
Changed line 169 from:
for i in range(1,n):
to:
for i in range(1,loops):
Changed line 174 from:
time.sleep(sleep-0.01)
to:
time.sleep(sleep)
Changed lines 183-197 from:
to:
# Read temperatures in Kelvin
T1[i] = a.T1
T2[i] = a.T2
###############################
### MPC CONTROLLER ###
###############################
TC1.MEAS = T1[i]
TC2.MEAS = T2[i]
# input setpoint with deadband +/- DT
DT = 0.1
TC1.SPHI = Tsp1[i] + DT
TC1.SPLO = Tsp1[i] - DT
TC2.SPHI = Tsp2[i] + DT
TC2.SPLO = Tsp2[i] - DT
# solve MPC
m.solve(disp=False)
# test for successful solution
if (m.options.APPSTATUS==1):
# retrieve the first Q value
Q1s[i] = Q1.NEWVAL
Q2s[i] = Q2.NEWVAL
else:
# not successful, set heater to zero
Q1s[i] = 0
Q2s[i] = 0
Changed lines 211-224 from:
# Predict Parameters and Temperatures with MHE
# use remote=False for local solve
m.solve()
if m.options.APPSTATUS == 1:
# Retrieve new values
Q1s[i] = Q1.NEWVAL
Q2s[i] = Q2.NEWVAL
else:
# Solution failed
Q1s[i] = 0.0
Q2s[i] = 0.0
# Write new heater values (0-100)
# use remote=False for local solve
m.solve()
if m.options.APPSTATUS == 1:
# Retrieve new values
Q1s[i] = Q1.NEWVAL
Q2s[i] = Q2.NEWVAL
else:
# Solution failed
Q1s[i] = 0.0
Q2s[i] = 0.0
# Write new heater values
to:
# Write output (0-100)
Changed line 217 from:
ax=plt.subplot(2,1,1)
to:
ax=plt.subplot(3,1,1)
Changed lines 219-222 from:
plt.plot(tm[0:i],T1m[0:i],'ro',label=r'$T_1$ measured')
plt.plot(tm[0:i],T1sp[0:i],'k-',label=r'$T_1$ setpoint')
plt.plot(tm[0:i],T2m[0:i],'bx',label=r'$T_2$ measured')
plt.plot(tm[0:i],T2sp[0:i],'k--',label=r'$T_2$ setpoint')
plt.plot(tm[0:i],
plt.plot(tm[0:i],T2m[0:i],'bx',label=r'$T_2$ measured')
plt.plot(tm[0:i],T2sp[0:i],'k--',label=r'$T_2$ setpoint
to:
plt.plot(tm[0:i],T1[0:i],'ro',MarkerSize=3,label=r'$T_1$')
plt.plot(tm[0:i],Tsp1[0:i],'k-',LineWidth=2,label=r'$T_1 Setpoint$')
plt.plot(tm[0:i],Tsp1[0:i],'k-',LineWidth=2,label=r'$T_1 Setpoint$')
Changed lines 222-223 from:
plt.legend(loc=2)
ax=plt.subplot(2,1,2)
ax=plt.subplot(
to:
plt.legend(loc='best')
ax=plt.subplot(3,1,2)
ax=plt.subplot(3,1,2)
Changed lines 225-228 from:
plt.plot(tm[0:i],Q1s[0:i],'r-',label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.plot(tm[0:i],
plt.ylabel('
plt.xlabel
to:
plt.plot(tm[0:i],T2[0:i],'ro',MarkerSize=3,label=r'$T_2$')
plt.plot(tm[0:i],Tsp2[0:i],'g-',LineWidth=2,label=r'$T_2 Setpoint$')
plt.ylabel('Temperature (degC)')
plt.plot(tm[0:i],Tsp2[0:i],'g-',LineWidth=2,label=r'$T_2 Setpoint$')
plt.ylabel('Temperature (degC)')
Added lines 229-235:
ax=plt.subplot(3,1,3)
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',LineWidth=3,label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',LineWidth=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',LineWidth=3,label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',LineWidth=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
Changed line 239 from:
# Turn off heaters and close connection
to:
# Turn off heaters
Changed lines 242-245 from:
plt.savefig('tclab_mpc.png')
to:
print('Shutting down')
Changed line 246 from:
# Turn off heaters and close connection
to:
# Disconnect from Arduino
Added line 249:
print('Shutting down')
Deleted lines 250-251:
plt.savefig('tclab_mpc.png')
Added line 257:
print('Error: Shutting down')
Deleted lines 258-259:
plt.savefig('tclab_mpc.png')
Changed line 1 from:
(:title TCLab G - Nonlinear Model Predictive Control:)
to:
(:title TCLab G - Nonlinear MPC:)
Added lines 1-308:
(:title TCLab G - Nonlinear Model Predictive Control:)
(:keywords Arduino, Empirical, Hybrid, MPC, Regression, temperature, control, process control, course:)
(:description Nonlinear Model Predictive Control for Temperature Control with Arduino Temperature Sensors and Heaters in MATLAB and Python:)
The TCLab is a hands-on application of machine learning and advanced temperature control with two heaters and two temperature sensors. The labs reinforce principles of model development, estimation, and advanced control methods. This is the seventh exercise and it involves nonlinear model predictive control with an energy balance model that is augmented with empirical elements. The predictions were previously aligned to the measured values through an estimator. This model predictive controller uses those parameters and a nonlinear MIMO (Multiple Input, Multiple Output) model of the TCLab input to output response to control temperatures to a set point.
!!!! Lab Problem Statement
* [[Attach:Lab_G_Nonlinear_MPC.pdf|Lab G - Nonlinear Model Predictive Control]]
!!!! Data and Solutions
* [[https://youtu.be/eoZRcbilKTU|Solution with MATLAB and Python]]
(:html:)
<iframe width="560" height="315" src="https://www.youtube.com/embed/eoZRcbilKTU" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
(:htmlend:)
----
%width=550px%Attach:tclab_hybrid_mpc.png
(:toggle hide gekko_labG button show="Lab G: Python TCLab MIMO MPC with Hybrid Model":)
(:div id=gekko_labG:)
(:source lang=python:)
import numpy as np
import time
import matplotlib.pyplot as plt
import random
# get gekko package with:
# pip install gekko
from gekko import GEKKO
# get tclab package with:
# pip install tclab
from tclab import TCLab
# Connect to Arduino
a = TCLab()
# Final time
tf = 10 # min
# number of data points (every 3 seconds)
n = tf * 20 + 1
# Percent Heater (0-100%)
Q1s = np.zeros(n)
Q2s = np.zeros(n)
# Temperatures (degC)
T1m = a.T1 * np.ones(n)
T2m = a.T2 * np.ones(n)
# Temperature setpoints
T1sp = T1m[0] * np.ones(n)
T2sp = T2m[0] * np.ones(n)
# Heater set point steps about every 150 sec
T1sp[3:] = 40.0
T2sp[40:] = 30.0
T1sp[80:] = 32.0
T2sp[120:] = 35.0
T1sp[160:] = 45.0
#########################################################
# Initialize Model as Estimator
#########################################################
# use remote=True for MacOS
m = GEKKO(name='tclab-mpc',remote=False)
# Control horizon, non-uniform time steps
m.time = [0,3,6,10,14,18,22,27,32,38,45,55,65, \
75,90,110,130,150]
# Parameters from Estimation
K1 = m.FV(value=0.607)
K2 = m.FV(value=0.293)
K3 = m.FV(value=0.24)
tau12 = m.FV(value=192)
tau3 = m.FV(value=15)
# don't update parameters with optimizer
K1.STATUS = 0
K2.STATUS = 0
K3.STATUS = 0
tau12.STATUS = 0
tau3.STATUS = 0
# Measured inputs
Q1 = m.MV(value=0)
Q1.STATUS = 1 # manipulated
Q1.FSTATUS = 0 # measured
Q1.DMAX = 20.0
Q1.DCOST = 0.1
Q1.UPPER = 100.0
Q1.LOWER = 0.0
Q2 = m.MV(value=0)
Q2.STATUS = 1 # manipulated
Q2.FSTATUS = 0 # measured
Q2.DMAX = 30.0
Q2.DCOST = 0.1
Q2.UPPER = 100.0
Q2.LOWER = 0.0
# State variables
TH1 = m.SV(value=T1m[0])
TH2 = m.SV(value=T2m[0])
# Measurements for model alignment
TC1 = m.CV(value=T1m[0])
TC1.STATUS = 1 # minimize error
TC1.FSTATUS = 1 # receive measurement
TC1.TAU = 10 # response speed (time constant)
TC1.TR_INIT = 2 # reference trajectory
TC2 = m.CV(value=T2m[0])
TC2.STATUS = 1 # minimize error
TC2.FSTATUS = 1 # receive measurement
TC2.TAU = 10 # response speed (time constant)
TC2.TR_INIT = 2 # reference trajectory
Ta = m.Param(value=23.0) # degC
# Heat transfer between two heaters
DT = m.Intermediate(TH2-TH1)
# Empirical correlations
m.Equation(tau12 * TH1.dt() + (TH1-Ta) == K1*Q1 + K3*DT)
m.Equation(tau12 * TH2.dt() + (TH2-Ta) == K2*Q2 - K3*DT)
m.Equation(tau3 * TC1.dt() + TC1 == TH1)
m.Equation(tau3 * TC2.dt() + TC2 == TH2)
# Global Options
m.options.IMODE = 6 # MHE
m.options.EV_TYPE = 1 # Objective type
m.options.NODES = 3 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.options.COLDSTART = 1 # COLDSTART on first cycle
##################################################################
# Create plot
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
# Main Loop
start_time = time.time()
prev_time = start_time
tm = np.zeros(n)
try:
for i in range(1,n):
# Sleep time
sleep_max = 3.0
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep-0.01)
else:
time.sleep(0.01)
# Record time and change in time
t = time.time()
dt = t - prev_time
prev_time = t
tm[i] = t - start_time
# Read temperatures in Celsius
T1m[i] = a.T1
T2m[i] = a.T2
# Insert measurements
TC1.MEAS = T1m[i]
TC2.MEAS = T2m[i]
# Adjust setpoints
TC1.SPHI = T1sp[i] + 0.1
TC1.SPLO = T1sp[i] - 0.1
TC2.SPHI = T2sp[i] + 0.1
TC2.SPLO = T2sp[i] - 0.1
# Predict Parameters and Temperatures with MHE
# use remote=False for local solve
m.solve()
if m.options.APPSTATUS == 1:
# Retrieve new values
Q1s[i] = Q1.NEWVAL
Q2s[i] = Q2.NEWVAL
else:
# Solution failed
Q1s[i] = 0.0
Q2s[i] = 0.0
# Write new heater values (0-100)
a.Q1(Q1s[i])
a.Q2(Q2s[i])
# Plot
plt.clf()
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm[0:i],T1m[0:i],'ro',label=r'$T_1$ measured')
plt.plot(tm[0:i],T1sp[0:i],'k-',label=r'$T_1$ setpoint')
plt.plot(tm[0:i],T2m[0:i],'bx',label=r'$T_2$ measured')
plt.plot(tm[0:i],T2sp[0:i],'k--',label=r'$T_2$ setpoint')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.draw()
plt.pause(0.05)
# Turn off heaters and close connection
a.Q1(0)
a.Q2(0)
a.close()
# Save figure
plt.savefig('tclab_mpc.png')
# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
# Turn off heaters and close connection
a.Q1(0)
a.Q2(0)
a.close()
print('Shutting down')
plt.savefig('tclab_mpc.png')
# Make sure serial connection still closes when there's an error
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
a.close()
print('Error: Shutting down')
plt.savefig('tclab_mpc.png')
raise
(:sourceend:)
(:divend:)
----
See also:
[[http://apmonitor.com/do/index.php/Main/AdvancedTemperatureControl|Advanced Control Lab Overview]]
-> [[Main/TCLabA|Lab A - Single Heater Modeling]]
-> [[Main/TCLabB|Lab B - Dual Heater Modeling]]
-> [[Main/TCLabC|Lab C - Parameter Estimation]]
-> [[Main/TCLabD|Lab D - Empirical Model Estimation]]
-> [[Main/TCLabE|Lab E - Hybrid Model Estimation]]
-> [[Main/TCLabF|Lab F - Linear Model Predictive Control]]
-> [[Main/TCLabG|Lab G - Nonlinear Model Predictive Control]]
-> [[Main/TCLabH|Lab H - Moving Horizon Estimation with MPC]]
[[https://gekko.readthedocs.io/en/latest/|GEKKO Documentation]]
[[https://tclab.readthedocs.io/en/latest/README.html|TCLab Documentation]]
[[https://github.com/APMonitor/arduino|TCLab Files on GitHub]]
[[https://apmonitor.com/heat.htm|Basic (PID) Control Lab]]
(:html:)
<style>
.button {
border-radius: 4px;
background-color: #0000ff;
border: none;
color: #FFFFFF;
text-align: center;
font-size: 28px;
padding: 20px;
width: 300px;
transition: all 0.5s;
cursor: pointer;
margin: 5px;
}
.button span {
cursor: pointer;
display: inline-block;
position: relative;
transition: 0.5s;
}
.button span:after {
content: '\00bb';
position: absolute;
opacity: 0;
top: 0;
right: -20px;
transition: 0.5s;
}
.button:hover span {
padding-right: 25px;
}
.button:hover span:after {
opacity: 1;
right: 0;
}
</style>
(:htmlend:)
(:keywords Arduino, Empirical, Hybrid, MPC, Regression, temperature, control, process control, course:)
(:description Nonlinear Model Predictive Control for Temperature Control with Arduino Temperature Sensors and Heaters in MATLAB and Python:)
The TCLab is a hands-on application of machine learning and advanced temperature control with two heaters and two temperature sensors. The labs reinforce principles of model development, estimation, and advanced control methods. This is the seventh exercise and it involves nonlinear model predictive control with an energy balance model that is augmented with empirical elements. The predictions were previously aligned to the measured values through an estimator. This model predictive controller uses those parameters and a nonlinear MIMO (Multiple Input, Multiple Output) model of the TCLab input to output response to control temperatures to a set point.
!!!! Lab Problem Statement
* [[Attach:Lab_G_Nonlinear_MPC.pdf|Lab G - Nonlinear Model Predictive Control]]
!!!! Data and Solutions
* [[https://youtu.be/eoZRcbilKTU|Solution with MATLAB and Python]]
(:html:)
<iframe width="560" height="315" src="https://www.youtube.com/embed/eoZRcbilKTU" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
(:htmlend:)
----
%width=550px%Attach:tclab_hybrid_mpc.png
(:toggle hide gekko_labG button show="Lab G: Python TCLab MIMO MPC with Hybrid Model":)
(:div id=gekko_labG:)
(:source lang=python:)
import numpy as np
import time
import matplotlib.pyplot as plt
import random
# get gekko package with:
# pip install gekko
from gekko import GEKKO
# get tclab package with:
# pip install tclab
from tclab import TCLab
# Connect to Arduino
a = TCLab()
# Final time
tf = 10 # min
# number of data points (every 3 seconds)
n = tf * 20 + 1
# Percent Heater (0-100%)
Q1s = np.zeros(n)
Q2s = np.zeros(n)
# Temperatures (degC)
T1m = a.T1 * np.ones(n)
T2m = a.T2 * np.ones(n)
# Temperature setpoints
T1sp = T1m[0] * np.ones(n)
T2sp = T2m[0] * np.ones(n)
# Heater set point steps about every 150 sec
T1sp[3:] = 40.0
T2sp[40:] = 30.0
T1sp[80:] = 32.0
T2sp[120:] = 35.0
T1sp[160:] = 45.0
#########################################################
# Initialize Model as Estimator
#########################################################
# use remote=True for MacOS
m = GEKKO(name='tclab-mpc',remote=False)
# Control horizon, non-uniform time steps
m.time = [0,3,6,10,14,18,22,27,32,38,45,55,65, \
75,90,110,130,150]
# Parameters from Estimation
K1 = m.FV(value=0.607)
K2 = m.FV(value=0.293)
K3 = m.FV(value=0.24)
tau12 = m.FV(value=192)
tau3 = m.FV(value=15)
# don't update parameters with optimizer
K1.STATUS = 0
K2.STATUS = 0
K3.STATUS = 0
tau12.STATUS = 0
tau3.STATUS = 0
# Measured inputs
Q1 = m.MV(value=0)
Q1.STATUS = 1 # manipulated
Q1.FSTATUS = 0 # measured
Q1.DMAX = 20.0
Q1.DCOST = 0.1
Q1.UPPER = 100.0
Q1.LOWER = 0.0
Q2 = m.MV(value=0)
Q2.STATUS = 1 # manipulated
Q2.FSTATUS = 0 # measured
Q2.DMAX = 30.0
Q2.DCOST = 0.1
Q2.UPPER = 100.0
Q2.LOWER = 0.0
# State variables
TH1 = m.SV(value=T1m[0])
TH2 = m.SV(value=T2m[0])
# Measurements for model alignment
TC1 = m.CV(value=T1m[0])
TC1.STATUS = 1 # minimize error
TC1.FSTATUS = 1 # receive measurement
TC1.TAU = 10 # response speed (time constant)
TC1.TR_INIT = 2 # reference trajectory
TC2 = m.CV(value=T2m[0])
TC2.STATUS = 1 # minimize error
TC2.FSTATUS = 1 # receive measurement
TC2.TAU = 10 # response speed (time constant)
TC2.TR_INIT = 2 # reference trajectory
Ta = m.Param(value=23.0) # degC
# Heat transfer between two heaters
DT = m.Intermediate(TH2-TH1)
# Empirical correlations
m.Equation(tau12 * TH1.dt() + (TH1-Ta) == K1*Q1 + K3*DT)
m.Equation(tau12 * TH2.dt() + (TH2-Ta) == K2*Q2 - K3*DT)
m.Equation(tau3 * TC1.dt() + TC1 == TH1)
m.Equation(tau3 * TC2.dt() + TC2 == TH2)
# Global Options
m.options.IMODE = 6 # MHE
m.options.EV_TYPE = 1 # Objective type
m.options.NODES = 3 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.options.COLDSTART = 1 # COLDSTART on first cycle
##################################################################
# Create plot
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
# Main Loop
start_time = time.time()
prev_time = start_time
tm = np.zeros(n)
try:
for i in range(1,n):
# Sleep time
sleep_max = 3.0
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep-0.01)
else:
time.sleep(0.01)
# Record time and change in time
t = time.time()
dt = t - prev_time
prev_time = t
tm[i] = t - start_time
# Read temperatures in Celsius
T1m[i] = a.T1
T2m[i] = a.T2
# Insert measurements
TC1.MEAS = T1m[i]
TC2.MEAS = T2m[i]
# Adjust setpoints
TC1.SPHI = T1sp[i] + 0.1
TC1.SPLO = T1sp[i] - 0.1
TC2.SPHI = T2sp[i] + 0.1
TC2.SPLO = T2sp[i] - 0.1
# Predict Parameters and Temperatures with MHE
# use remote=False for local solve
m.solve()
if m.options.APPSTATUS == 1:
# Retrieve new values
Q1s[i] = Q1.NEWVAL
Q2s[i] = Q2.NEWVAL
else:
# Solution failed
Q1s[i] = 0.0
Q2s[i] = 0.0
# Write new heater values (0-100)
a.Q1(Q1s[i])
a.Q2(Q2s[i])
# Plot
plt.clf()
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm[0:i],T1m[0:i],'ro',label=r'$T_1$ measured')
plt.plot(tm[0:i],T1sp[0:i],'k-',label=r'$T_1$ setpoint')
plt.plot(tm[0:i],T2m[0:i],'bx',label=r'$T_2$ measured')
plt.plot(tm[0:i],T2sp[0:i],'k--',label=r'$T_2$ setpoint')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.draw()
plt.pause(0.05)
# Turn off heaters and close connection
a.Q1(0)
a.Q2(0)
a.close()
# Save figure
plt.savefig('tclab_mpc.png')
# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
# Turn off heaters and close connection
a.Q1(0)
a.Q2(0)
a.close()
print('Shutting down')
plt.savefig('tclab_mpc.png')
# Make sure serial connection still closes when there's an error
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
a.close()
print('Error: Shutting down')
plt.savefig('tclab_mpc.png')
raise
(:sourceend:)
(:divend:)
----
See also:
[[http://apmonitor.com/do/index.php/Main/AdvancedTemperatureControl|Advanced Control Lab Overview]]
-> [[Main/TCLabA|Lab A - Single Heater Modeling]]
-> [[Main/TCLabB|Lab B - Dual Heater Modeling]]
-> [[Main/TCLabC|Lab C - Parameter Estimation]]
-> [[Main/TCLabD|Lab D - Empirical Model Estimation]]
-> [[Main/TCLabE|Lab E - Hybrid Model Estimation]]
-> [[Main/TCLabF|Lab F - Linear Model Predictive Control]]
-> [[Main/TCLabG|Lab G - Nonlinear Model Predictive Control]]
-> [[Main/TCLabH|Lab H - Moving Horizon Estimation with MPC]]
[[https://gekko.readthedocs.io/en/latest/|GEKKO Documentation]]
[[https://tclab.readthedocs.io/en/latest/README.html|TCLab Documentation]]
[[https://github.com/APMonitor/arduino|TCLab Files on GitHub]]
[[https://apmonitor.com/heat.htm|Basic (PID) Control Lab]]
(:html:)
<style>
.button {
border-radius: 4px;
background-color: #0000ff;
border: none;
color: #FFFFFF;
text-align: center;
font-size: 28px;
padding: 20px;
width: 300px;
transition: all 0.5s;
cursor: pointer;
margin: 5px;
}
.button span {
cursor: pointer;
display: inline-block;
position: relative;
transition: 0.5s;
}
.button span:after {
content: '\00bb';
position: absolute;
opacity: 0;
top: 0;
right: -20px;
transition: 0.5s;
}
.button:hover span {
padding-right: 25px;
}
.button:hover span:after {
opacity: 1;
right: 0;
}
</style>
(:htmlend:)