GEKKO MPC实现干扰抑制时被控变量无法回到稳态的问题问询

你当前代码实现的是开环最优控制逻辑，没有滚动反馈校正环节，这是阶跃干扰引入后被控变量无法回归稳态的核心原因，具体问题和修复方法如下：

核心问题点 #

• IMODE=6调用逻辑错误：单次调用solve()时，GEKKO会基于初始状态和预先定义的完整干扰序列一次性计算出全时间域的控制序列，不会在每个控制步引入实际被控量的测量值做校正，本质是开环运行，模型误差、未知干扰都会导致稳态偏差。
• CV参数配置错误：你将所有CV的FSTATUS设为1，代表启用外部测量值更新，但全程没有传入实际测量值，同时手动设置BIAS=1会引入固定预测偏差，无法抵消干扰影响。
• 干扰注入不符合实际运行逻辑：你预先把干扰d1的完整时间序列存入Param参数，相当于模型提前完全知晓干扰变化规律，实际场景下干扰是未知的，必须通过反馈校正抵消。

修复方案（闭环滚动MPC实现） #

需要改为滚动时域控制逻辑，每个控制步更新测量值、解算控制量、输出执行，示例实现代码如下：

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

# 配置参数
pred_steps = 30 # 预测时域步长
ctrl_steps = 100 # 总控制步数
dt = 3 # 控制周期3s，总时长300s和原设置一致

# 初始化GEKKO模型
m = GEKKO(remote=False)
m.time = np.linspace(0, pred_steps*dt, pred_steps)

# 系统参数（和原代码一致）
T1 = m.Param(value = 53.97272679974334)
T2 = m.Param(value = 48.06851424706475)
T3 = m.Param(value = 38.48651254747577)
T4 = m.Param(value = 31.018933652439845)
k1 = m.Param(value = 5.51)
k2 = m.Param(value = 6.58)
γ1bar = m.Param(value = 0.333)
γ2bar = m.Param(value = 0.307)
A1 = m.Param(value = 730)
A2 = m.Param(value = 730)
A3 = m.Param(value = 730)
A4 = m.Param(value = 730)
v1bar = m.Param(value = 60)
v2bar = m.Param(value = 60)

# 操纵变量
v1 = m.MV(value=0, lb=0, ub=100)
v1.STATUS = 1  
v2 = m.MV(value=0, lb=0, ub=100)
v2.STATUS = 1  
γ1 = m.MV(value=0, lb=0, ub=1)
γ1.STATUS = 1  
γ2 = m.MV(value=0, lb=0, ub=1)
γ2.STATUS = 1  

# 干扰（改为实时更新，不需要预定义全序列）
d1 = m.Param(value=0)
d2 = m.Param(value=0)

m.options.CV_TYPE = 2 # 平方误差
m.options.IMODE = 6 # 控制模式

# 被控变量
h1 = m.CV(value=0)
h1.STATUS = 1 
h1.SP = 1   
h1.TR_INIT = 0 # 滚动运行时关闭初始轨迹
h1.TAU = 1   
h1.FSTATUS = 1 # 接收测量值

h2 = m.CV(value=0)
h2.STATUS = 1  
h2.SP = 0  
h2.TR_INIT = 0
h2.TAU = 1     
h2.FSTATUS = 1

h3 = m.CV(value=0)
h3.STATUS = 1  
h3.SP = 0    
h3.TR_INIT = 0
h3.TAU = 1     
h3.FSTATUS = 1

h4 = m.CV(value=0)
h4.STATUS = 1  
h4.SP = 0    
h4.TR_INIT = 0
h4.TAU = 1     
h4.FSTATUS = 1

# 状态方程（和原代码一致）
m.Equation(h1.dt() == -(1/T1)*h1 + (A3/(A1*T3))*h3 + (γ1bar*k1*v1)/A1 + (γ1*k1*v1bar)/A1)
m.Equation(h2.dt() == -(1/T2)*h2 + (A4/(A2*T4))*h4 + (γ2bar*k2*v2)/A2 + (γ2*k2*v2bar)/A2)
m.Equation(h3.dt() == -(1/T3)*h3 + ((1-γ2bar)*k2*v2)/A3 - k2*v2bar*γ2/A3 - (k1*d1)/A3) 
m.Equation(h4.dt() == -(1/T4)*h4 + ((1-γ1bar)*k1*v1)/A4 - k1*v1bar*γ1/A4 - (k2*d2)/A4)

# 存储仿真结果
h1_store = []
h2_store = []
h3_store = []
h4_store = []
v1_store = []
v2_store = []
gamma1_store = []
gamma2_store = []

# 真实被控对象仿真（模拟实际系统，和模型完全一致用于验证）
def real_system(h1_prev, h2_prev, h3_prev, h4_prev, v1_val, v2_val, gamma1_val, gamma2_val, d1_val, d2_val, dt):
    dh1 = -(1/T1.value[0])*h1_prev + (A3.value[0]/(A1.value[0]*T3.value[0]))*h3_prev + (γ1bar.value[0]*k1.value[0]*v1_val)/A1.value[0] + (gamma1_val*k1.value[0]*v1bar.value[0])/A1.value[0]
    dh2 = -(1/T2.value[0])*h2_prev + (A4.value[0]/(A2.value[0]*T4.value[0]))*h4_prev + (γ2bar.value[0]*k2.value[0]*v2_val)/A2.value[0] + (gamma2_val*k2.value[0]*v2bar.value[0])/A2.value[0]
    dh3 = -(1/T3.value[0])*h3_prev + ((1-γ2bar.value[0])*k2.value[0]*v2_val)/A3.value[0] - k2.value[0]*v2bar.value[0]*gamma2_val/A3.value[0] - (k1.value[0]*d1_val)/A3.value[0]
    dh4 = -(1/T4.value[0])*h4_prev + ((1-γ1bar.value[0])*k1.value[0]*v1_val)/A4.value[0] - k1.value[0]*v1bar.value[0]*gamma1_val/A4.value[0] - (k2.value[0]*d2_val)/A4.value[0]
    h1_new = h1_prev + dh1*dt
    h2_new = h2_prev + dh2*dt
    h3_new = h3_prev + dh3*dt
    h4_new = h4_prev + dh4*dt
    return h1_new, h2_new, h3_new, h4_new

# 初始状态
h1_real = 0
h2_real = 0
h3_real = 0
h4_real = 0

# 滚动控制循环
for i in range(ctrl_steps):
    # 注入阶跃干扰：第10步（30s）后d1=1
    if i >= 10:
        d1.value = 1
    else:
        d1.value = 0
    
    # 更新CV测量值
    h1.MEAS = h1_real
    h2.MEAS = h2_real
    h3.MEAS = h3_real
    h4.MEAS = h4_real
    
    # 解算控制量
    m.solve(disp=False)
    
    # 取第一个时刻的控制量输出
    v1_out = v1.value[0]
    v2_out = v2.value[0]
    gamma1_out = γ1.value[0]
    gamma2_out = γ2.value[0]
    
    # 存储结果
    h1_store.append(h1_real)
    h2_store.append(h2_real)
    h3_store.append(h3_real)
    h4_store.append(h4_real)
    v1_store.append(v1_out)
    v2_store.append(v2_out)
    gamma1_store.append(gamma1_out)
    gamma2_store.append(gamma2_out)
    
    # 模拟实际系统运行一个周期
    h1_real, h2_real, h3_real, h4_real = real_system(h1_real, h2_real, h3_real, h4_real, v1_out, v2_out, gamma1_out, gamma2_out, d1.value[0], 0, dt)

验证说明 #

修改后MPC会每个控制步引入实际被控量的测量值校正预测偏差，阶跃干扰引入后被控变量可以逐步回归设定值，无稳态偏差。

内容的提问来源于stack exchange，提问作者Cameron Ernst Bolt