If an MPC application is giving unexpected results, it can typically be diagnosed by showing the unbiased and biased model predictions and the reference trajectory. Here are a couple things that you should consider with you model.
- Change the equation to be in deviation variable form. The steady state point for temperature is
TC1=TC1_ss and Q1d=0. The updated equation is:
m.Equation(tau * TC1.dt() + (TC1-TC1_ss) == Kp * (Q1d-0))
- Use a steady-state simulation to initialize all of the variables with
IMODE=1.
# Steady-state initialization
m.options.IMODE = 1
m.solve()
# Global Options
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # Objective type
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 1 # 1=APOPT, 3=IPOPT
This helps the problem but does not eliminate it. Anytime there is plant / model mismatch, the MPC will be too sluggish if the model gain is too high or create oscillations if the model gain is too low. The oscillations are because the MPC has updated measurements and has to recalculate the move plan based on updated information.
import tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from gekko import GEKKO
import _thread
import queue
import threading
import json
# Connect to Arduino
a = tclab.TCLabModel()
R = threading.Lock()
# Turn LED on
print('LED On')
a.LED(100)
# Simulate a time delay
delay = 40
# Run time in minutes
run_time = 60.0
# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)
# Temperature (K)
T1 = np.ones(loops) * a.T1 # temperature (degC)
Tsp1 = np.ones(loops) * 40.0 # set point (degC)
# heater values
Q1s = np.ones(loops) * 0.0
#########################################################
# Initialize Model
#########################################################
# use remote=True for MacOS
m = GEKKO(name='tclab-mpc',remote=False)
# 30 second time horizon
m.time = np.linspace(0,240,121)
# Parameters
Q1_ss = m.Param(value=20)
TC1_ss = m.Param(value=a.T1)
Kp = m.Param(value=1.5)
tau = m.Param(value=160.0)
# Manipulated variable
Q1 = m.MV(value=Q1_ss.VALUE, name='q1')
Q1.STATUS = 1 # use to control temperature
Q1.FSTATUS = 0 # no feedback measurement
Q1.LOWER = 0.0
Q1.UPPER = 100.0
Q1.DMAX = 2
Q1.DCOST = 10.0
# Controlled variable
TC1 = m.CV(value=TC1_ss.VALUE, name='tc1')
TC1.STATUS = 1 # minimize error with setpoint range
TC1.FSTATUS = 1 # receive measurement
TC1.TR_INIT = 2 # reference trajectory
# TC1.TAU = 100 # time constant for response
# 添加延时
Q1d=m.Var()
m.delay(Q1, Q1d, 40)
# Equation
m.Equation(tau * TC1.dt() + (TC1-TC1_ss) == Kp * Q1d)
# Steady-state initialization
m.options.IMODE = 1
m.solve()
# Global Options
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # Objective type
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 1 # 1=APOPT, 3=IPOPT
##################################################################
# Create plot
plt.figure()
plt.ion()
plt.show()
filter_tc1 = []
q = queue.Queue()
flag = False
def setQ1():
global flag
global a
last_time = time.time()
while True:
sleep_max = 2.0
sleep = sleep_max - (time.time() - last_time)
print("q.qsize()", q.qsize())
if sleep >= 0.01:
time.sleep(sleep)
else:
time.sleep(0.01)
# Record time and change in time
last_time = time.time()
if q.qsize() >= delay:
q1 = q.get()
print(f'Q1: {q1} ')
try:
R.acquire()
print("write Q1：", a.Q1(float(q1)))
R.release()
except:
print("转换报错！")
pass
_thread.start_new_thread(setQ1, ())
# Rolling average filtering
def movefilter(predata, new, n):
if len(predata) < n:
predata.append(new)
else:
predata.pop(0)
predata.append(new)
return np.average(predata)
# Main Loop
start_time = time.time()
prev_time = start_time
last_Q1 = 0
try:
for i in range(1,loops):
# Sleep time
sleep_max = 2.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
R.acquire()
curr_T1 = a.T1
R.release()
last_T1 = curr_T1
avg_T1 = movefilter(filter_tc1, last_T1, 3)
T1[i] = curr_T1
# T1[i] = a.T1
# T2[i] = a.T2
###############################
### MPC CONTROLLER ###
###############################
TC1.MEAS = avg_T1
# input setpoint with deadband +/- DT
DT = 0.5
TC1.SPHI = Tsp1[i] + DT
TC1.SPLO = Tsp1[i] - DT
try:
# solve MPC
m.solve(disp=False)
except:
pass
# test for successful solution
if (m.options.APPSTATUS==1):
# retrieve the first Q value
Q1s[i] = Q1.NEWVAL
with open(m.path+'//results.json') as f:
results = json.load(f)
else:
# not successful, set heater to zero
Q1s[i] = 0
# Write output (0-100)
q.put(Q1s[i])
# Plot
plt.clf()
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm[0:i],T1[0:i],'ro',markersize=3,label=r'$T_1$')
plt.plot(tm[0:i],Tsp1[0:i],'b-',markersize=3,label=r'$T_1 Setpoint$')
plt.plot(tm[i]+m.time,results['tc1.bcv'],'k-.',
label=r'$T_1$ predicted',linewidth=3)
plt.plot(tm[i]+m.time,results['tc1.tr_hi'],'k--',
label=r'$T_1$ trajectory')
plt.plot(tm[i]+m.time,results['tc1.tr_lo'],'k--')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',linewidth=3,label=r'$Q_1$')
plt.plot(tm[i]+m.time,Q1.value,'k-.',
label=r'$Q_1$ plan',linewidth=3)
# 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')
a.close()
# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Shutting down')
a.close()
# 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
Figure 14 of this paper quantifies the MPC performance with model mismatch.
- Hedengren, J. D., Eaton, A. N., Overview of Estimation Methods for Industrial Dynamic Systems, Optimization and Engineering, Springer, Vol 18 (1), 2017, pp. 155-178, DOI: 10.1007/s11081-015-9295-9.
