如何解决为什么延迟参数会影响预测曲线?
我将 Gekko 的计算结果放入队列,延迟后,将队列的第一个值写入 TCLab Arduino。我用这种方法模拟了一个工厂大时延系统,然后我优化了 Gekko 参数以达到更好的控制效果。当我在模型中添加一个 delay=1 时,我得到了一个非常好的预测曲线:
最后的控制效果也不错:
但是当我设置delay=80时,没有修改其他参数,预测曲线不理想:
最后的控制效果也不好:
为什么延迟参数会影响预测曲线?我认为参考轨迹也应该随时间延迟而变化。我该如何解决这个问题?
代码如下:
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 = 80
# 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)
# 240 second time horizon
m.time = np.linspace(0,240,121)
# Parameters
Q1_ss = m.Param(value=0)
TC1_ss = m.Param(value=a.T1)
Kp = m.Param(value=1.0)
tau = m.Param(value=120.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 = 50
Q1.dcosT = 2.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 = 30 # time constant for response
# 添加延时
Q1d=m.Var()
m.delay(Q1,Q1d,delay)
# 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()
while q.qsize() >= delay:
q1 = q.get()
print(f'Q1: {q1} ')
try:
R.acquire()
print("write Q1:",a.Q1(float(q1)))
R.release()
except:
print("convertion error!")
_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)
print(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 = 1.0
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-',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,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,\
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
解决方法
控制器按预期工作。这个问题是一个众所周知的模型不匹配和延迟问题。经典的控制方法是使用史密斯预测器来补偿延迟。使用 MPC,延迟被内置到模型中,类似于 Smith Predictor 的功能。如果您将模型更改为:
Kp = m.Param(value=0.7)
tau = m.Param(value=180.0)
然后性能得到改善:
对参考轨迹的修改也提高了性能,因为控制器专注于长期目标而不是它无法控制的近期误差。
方法 1 - 在开始时打开参考轨迹
TC1.TR_INIT = 2 # reference trajectory
TC1.TR_OPEN = 20 # more focus on final destination
TC1.TAU = 60 # time constant for response
方法 2 - 稍后开始惩罚
使用 CV_WGT_START 仅在一定时间延迟后进行惩罚。
TC1.TR_INIT = 0 # reference trajectory
m.options.CV_WGT_START = 80
参数需要再调整一次以更好地匹配过程。
Kp = m.Param(value=0.8)
tau = m.Param(value=160.0)
方法 3:定义自己的参考轨迹
如果您确实需要自定义参考轨迹,那么您可以定义一个延迟为 SP
的新变量 SPd
。定义一个新的 CV error
,它是 TC1
和延迟参考轨迹之间的差值。
SP = m.Var()
m.Equation(30 * SP.dt() == -SP + 40)
SPd = m.Var()
m.delay(SP,SPd,delay)
error = m.CV(value=0)
m.Equation(error==SPd - TC1)
error.STATUS = 1 # minimize error with setpoint range
error.FSTATUS = 0 # receive measurement
error.TR_INIT = 0 # reference trajectory
error.SPHI = 0.5 # drive error close to 0
error.SPLO = -0.5 # drive error close to 0
使用这种方法,您需要找到一种方法来更新 error
测量以获得测量反馈。
也有很多时候 MPC 计算稍微超出了周期时间。这会造成模型不匹配,您可以通过将 MPC 周期时间增加到 3 秒来帮助解决此问题。
修改后的脚本
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 = 80
# 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)
# 240 second time horizon
m.time = np.linspace(0,240,121)
# Parameters
Q1_ss = m.Param(value=0)
TC1_ss = m.Param(value=a.T1)
Kp = m.Param(value=0.7)
tau = m.Param(value=180.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 = 50
Q1.DCOST = 2.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 = 0 # reference trajectory
m.options.CV_WGT_START = 80
# 添加延时
Q1d=m.Var()
m.delay(Q1,Q1d,delay)
# 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()
while q.qsize() >= delay:
q1 = q.get()
print(f'Q1: {q1} ')
try:
R.acquire()
print("write Q1:",a.Q1(float(q1)))
R.release()
except:
print("convertion error!")
_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)
print(time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep)
else:
print('Warning: Dynamic mismatch from excess MPC time')
print(' Consider increasing the cycle time to 3+ sec')
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 = 1.0
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-',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,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,\
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
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。