基于Gekko的MPC多目标控制：RNA-MEAS偏差校正及动力学求解求助

一、将RNA作为MEAS接入MPC实现偏差校正 #

Gekko的MEAS功能用于引入测量值校正模型预测偏差，针对你的RNA虚拟分析仪，按以下步骤操作：

• 定义测量参数：将RNA的输出作为m.Param()传入Gekko，实时运行时每次求解前更新该参数的value即可。
• 切换变量类型：把被控变量（CV1/CV2/CV3）从m.Var()改为m.CV()，设置.MEAS属性为RNA输出值，同时开启.STATUS=1让MPC启用测量校正。
• 配置校正参数：通过.TAU调整校正速度，值越小校正越快，通常设为过程时间常数的1/10~1/5。

示例代码片段：

# 模拟RNA预测逻辑（实际替换为你的训练好的模型）
def rna_predict(B_val, h2_val, T_val, E_val, DT_val, H1_val):
    cv1_pred = 1000*(-0.149776*B_val+0.300058*B_val**2-0.237433*B_val**3+0.0001952*(h2_val)+0.000149*(T_val)+0.91439)
    cv2_pred = math.exp(2.05433906*B_val+0.08761687*(h2_val)+0.19555572*E_val-0.02127295*H1_val+0.02838329*DT_val-5.59437532)
    cv3_pred = (((8.6576*E_val + 112.66-(30+4))/(((16.15-0.019*(8.6576*E_val + 112.66+30)/2))*E_val))-0.06*B_val)*100
    return cv1_pred, cv2_pred, cv3_pred

# 生成模拟RNA测量数据
rna_cv1, rna_cv2, rna_cv3 = [], [], []
for i in range(len(t)):
    cv1p, cv2p, cv3p = rna_predict(0.12, 1, T1[i], 20, 8.6576*20+112.66-30, 230)
    rna_cv1.append(cv1p)
    rna_cv2.append(cv2p)
    rna_cv3.append(cv3p)

# 定义测量参数
meas_cv1 = m.Param(value=rna_cv1)
meas_cv2 = m.Param(value=rna_cv2)
meas_cv3 = m.Param(value=rna_cv3)

# 定义CV并启用MEAS
CV1 = m.CV(value=rna_cv1[0])
CV1.STATUS = 1
CV1.MEAS = meas_cv1
CV1.TAU = 2  # 校正时间常数
CV1.SP = m.Param(value=sp_d)

CV2 = m.CV(value=rna_cv2[0])
CV2.STATUS = 1
CV2.MEAS = meas_cv2
CV2.TAU = 2
CV2.SP = m.Param(value=sp_m)

CV3 = m.CV(value=rna_cv3[0])
CV3.STATUS = 1
CV3.MEAS = meas_cv3
CV3.TAU = 2
CV3.SP = m.Param(value=sp_q)

二、添加动力学后无可行解的问题分析与解决 #

你的代码存在几个核心问题导致无解：

1. 代数环问题：h2 -> k_mi -> h2和B -> k_d -> B的循环依赖会让求解器无法收敛。解决方法是将这类循环转换为动态方程：

   # 将h2的代数方程改为动态形式
   m.Equation(m.dt(h2) == (40*k_mi*(1-m.exp(-(ti_m-tau_m_d)/tau_m)) - h2)/tau_m)
   # 将k_d的循环改为动态更新
   m.Equation(m.dt(k_d) == (0.0909/B - k_d)/tau_d)
   

2. MV约束过严/冗余：MV数量（4个）多于CV数量（3个），且部分MV的DMAX设置过小，限制了求解空间。可以：

   • 移除冗余MV（比如取消ti_d的MV属性）
   • 放宽DMAX参数（如ti_m的DMAX从2改为5）
   • 检查MV上下限是否符合实际过程范围

3. 初始值不可行：部分变量初始值未满足方程约束（如CV2未设初始值），导致求解器初始点无效。需给所有Var/CV设置符合方程的初始值。

三、修改后的完整代码示例 #

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt  
import math 

m = GEKKO(remote=False)
t= np.linspace(0,500,100)
m.time=t

m.options.CV_TYPE = 2 
m.options.IMODE = 6  # MPC模式
m.options.SOLVER = 2
m.options.COLDSTART = 1  # 冷启动帮助找到初始可行解

# 参数定义
T1 = np.ones(len(t))*260
T1[90:] = 265

tau_m_d= m.Param(value=1)
tau_m= m.Param(value=5)
tau_d_d= m.Param(value=5)
tau_d= m.Param(value=15)
T = m.Param(value=T1)
Ti=m.Param(value=30)
H1=m.Param(value=230)

# 操纵变量
B = m.MV(value=0.12, lb=0, ub=0.8) 
B.STATUS = 1  
B.DCOST = 0.01 
B.DMAX = 0.2

ti_m = m.MV(value=1, lb=0, ub=300)
ti_m.STATUS = 1  
ti_m.DCOST = 0.1 
ti_m.DMAX = 5  # 放宽移动约束

E = m.MV(value=20, lb=18, ub=22)
E.STATUS = 1  
E.DCOST = 0.1 
E.DMAX = 0.2

To = m.Intermediate(8.6576*E + 112.66)

# 设置点
sp_d = np.ones(len(t))*940
sp_d[20:] =945
sp_d[60:]=930

sp_m = np.ones(len(t))*4
sp_m[40:] = 40
sp_m[80:]= 10

sp_q = np.ones(len(t))*94 
sp_q[15:] =94.5
sp_q[70:]=95

# 模拟RNA测量值
def rna_predict(B_val, h2_val, T_val, E_val, DT_val, H1_val):
    cv1_pred = 1000*(-0.149776*B_val+0.300058*B_val**2-0.237433*B_val**3+0.0001952*(h2_val)+0.000149*(T_val)+0.91439)
    cv2_pred = math.exp(2.05433906*B_val+0.08761687*(h2_val)+0.19555572*E_val-0.02127295*H1_val+0.02838329*DT_val-5.59437532)
    cv3_pred = (((8.6576*E_val + 112.66-(30+4))/(((16.15-0.019*(8.6576*E_val + 112.66+30)/2))*E_val))-0.06*B_val)*100
    return cv1_pred, cv2_pred, cv3_pred

rna_cv1, rna_cv2, rna_cv3 = [], [], []
for i in range(len(t)):
    cv1p, cv2p, cv3p = rna_predict(0.12, 1, T1[i], 20, 8.6576*20+112.66-30, 230)
    rna_cv1.append(cv1p)
    rna_cv2.append(cv2p)
    rna_cv3.append(cv3p)

# 测量参数
meas_cv1 = m.Param(value=rna_cv1)
meas_cv2 = m.Param(value=rna_cv2)
meas_cv3 = m.Param(value=rna_cv3)

# 被控变量与中间变量
h2 = m.Var(value=1, lb=0, ub=300)
CV1 = m.CV(value=rna_cv1[0]) 
CV1.STATUS = 1
CV1.MEAS = meas_cv1
CV1.TAU = 2
CV1.SP = m.Param(value=sp_d)

CV2 = m.CV(value=rna_cv2[0]) 
CV2.STATUS = 1
CV2.MEAS = meas_cv2
CV2.TAU = 2
CV2.SP = m.Param(value=sp_m)

DT = m.Var(value=8.6576*20+112.66-30) 
k_mi = m.Var(value=0.4)
k_d = m.Var(value=0.0909/0.12)
CV3 = m.CV(value=rna_cv3[0])
CV3.STATUS = 1
CV3.MEAS = meas_cv3
CV3.TAU = 2
CV3.SP = m.Param(value=sp_q)

# 过程方程（消除代数环）
m.Equation(m.dt(h2) == (40*k_mi*(1-m.exp(-(ti_m-tau_m_d)/tau_m)) - h2)/tau_m)
m.Equation(m.dt(k_d) == (0.0909/B - k_d)/tau_d)
m.Equation(CV1==(1000*(-0.149776*B+0.300058*B**2-0.237433*B**3+0.0001952*(h2)+0.000149*(T)+0.91439)))
m.Equation(CV2==m.exp(2.05433906*B+0.08761687*(h2)+0.19555572*E-0.02127295*H1+0.02838329*DT-5.59437532)) 
m.Equation(CV3==(((To-(Ti+4))/(((16.15-0.019*(To+Ti)/2))*E))-0.06*B)*100)
m.Equation(k_mi==CV2*0.1)
m.Equation(DT==To-Ti)

m.solve(disp=True)

# 绘图
plt.figure(figsize=(15,20))
plt.subplot(6,1,1)
plt.plot(m.time,B.value,'b-',label='MV B')
plt.ylabel('B')
plt.legend(loc='best')

plt.subplot(6,1,2)
plt.plot(m.time,CV1.value,'k-',label='CV1 Response')
plt.plot(m.time,sp_d,'r--',label='Setpoint')
plt.plot(m.time,rna_cv1,'g:',label='RNA Measurement')
plt.legend(loc='best')
plt.ylabel('CV1')

plt.subplot(6,1,3)
plt.plot(m.time,E.value,'b-',label='MV E')
plt.legend(loc='best')
plt.ylabel('E')

plt.subplot(6,1,4)
plt.plot(m.time,CV3.value,'b-',label='CV3 Response')
plt.plot(m.time,sp_q,'r--',label='CV3 Setpoint')
plt.plot(m.time,rna_cv3,'g:',label='RNA Measurement')
plt.ylim([93, 96])
plt.ylabel('CV3')
plt.legend(loc='best')

plt.subplot(6,1,5)
plt.plot(m.time,h2.value,'b-',label='h2')
plt.legend(loc='best')
plt.ylabel('H2')

plt.subplot(6,1,6)
plt.plot(m.time,CV2.value,'k-',label='CV2 Response')
plt.plot(m.time,sp_m,'r--',label='Setpoint')
plt.plot(m.time,rna_cv2,'g:',label='RNA Measurement')
plt.ylabel('CV2')
plt.legend(loc='best')

plt.xlabel('Time (min)')
plt.show()

关键注意事项 #

• 代数环处理：必须消除循环依赖的代数方程，改用动态方程或重构变量关系。
• MEAS限制：只有CV类型变量支持.MEAS属性，Var类型无法使用该功能。
• 冷启动设置：启用COLDSTART=1可帮助求解器快速找到初始可行点，解决启动阶段无解问题。
• 参数调优：根据实际过程调整TAU、DCOST、DMAX等参数，平衡控制性能与求解稳定性。

内容的提问来源于stack exchange，提问作者Adilton Lopes da Silva