Gekko ARX模型输出与SysId回归结果不符问题排查

使用Python Gekko库的sysid工具构建动态过程模型后，验证环节出现异常：

• 经sysid得到系数字典p与回归输出yp，用p创建ARX模型并输入原始数据时，预测输出yc与实际输出y偏差极大
• 跳过阶跃测试代码直接执行验证，结果匹配正常；但运行阶跃测试后再执行验证，结果就会异常
• 重启Jupyter后首次验证结果正常，但运行包含阶跃测试的完整脚本后再次验证仍出问题，推测是模型初始化状态或残差残留导致。

系统辨识代码 #

# 加载数据并解析为列
url = 'http://apmonitor.com/do/uploads/Main/cstr_step_tests.txt'
data = pd.read_csv(url)
print(data.head())

# 生成时间序列模型
t = data['Time']
u = data['Tc']
y = data['T']
m = GEKKO(remote=False)

# 系统辨识
na = 2 # 输出系数阶数
nb = 2 # 输入系数阶数
yp,p,K = m.sysid(t,u,y,na,nb,shift='init',scale=True,objf=100,diaglevel=1)

# 绘制拟合结果
plt.figure()
plt.subplot(2,1,1)
plt.plot(t,u)
plt.legend([r'$T_c$'])
plt.ylabel('MV')
plt.subplot(2,1,2)
plt.plot(t,y)
plt.plot(t,yp)
plt.legend([r'$T_{meas}$',r'$T_{pred}$'])
plt.ylabel('CV')
plt.xlabel('Time')
plt.savefig('sysid.png')

阶跃测试代码 #

# 基于ARX模型构建阶跃测试模型
yc,uc = m.arx(p)

# 重命名操作变量与被控变量
Tc = uc[0]
T = yc[0]

# 稳态初始化
m.options.IMODE = 1
Tc.value = 300
m.solve(disp=False)

# 动态仿真（阶跃测试验证）
m.time = np.linspace(0,2,21)
m.options.IMODE = 4
Tc.value = np.ones(21)*300
Tc.value[5:] = 305
m.solve(disp=False)

plt.figure()
plt.subplot(2,1,1)
plt.title('阶跃测试')
plt.plot(m.time,Tc.value,'b-',label='冷却夹套温度')
plt.ylabel(r'$T_c (K)$')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m.time,T.value,'r-',label='反应器温度')
plt.ylabel('T (K)')
plt.xlabel('时间 (min)')
plt.legend()

plt.show()

验证代码（原异常版本） #

# 基于ARX模型构建验证模型
yc,uc = m.arx(p)

# 重命名操作变量与被控变量
Tc = uc[0]
T = yc[0]

# 稳态初始化
m.options.IMODE = 1
Tc.value = 300
m.solve(disp=False)

# 动态仿真
m.time = np.linspace(0,50,501)
Tc.value = data['Tc'].values
m.options.IMODE = 4
m.solve(disp=False)

plt.figure()
plt.subplot(2,1,1)
plt.title('模型验证')
plt.plot(m.time,Tc.value,'b-',label='冷却夹套温度')
plt.ylabel(r'$T_c (K)$')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m.time,T.value,'r-',label='反应器温度')
plt.ylabel('T (K)')
plt.xlabel('时间 (min)')
plt.legend()

plt.show()

问题原因与解决方案 #

原因 #

Gekko是有状态的建模环境，同一个Gekko实例m在运行阶跃测试后，内部的时间序列、变量状态、模型结构会被修改并残留。后续再用这个实例创建ARX模型并执行验证时，残留的状态会干扰新的模型计算，导致预测结果异常。

解决方案 #

最稳妥的方式是为不同的建模/仿真任务创建独立的Gekko实例，避免状态交叉污染。修改后的验证代码如下：

# 新建独立的Gekko实例，彻底隔离状态
m_val = GEKKO(remote=False)

# 基于ARX模型构建验证模型
yc,uc = m_val.arx(p)

# 重命名操作变量与被控变量
Tc = uc[0]
T = yc[0]

# 稳态初始化
m_val.options.IMODE = 1
Tc.value = 300
m_val.solve(disp=False)

# 动态仿真
m_val.time = np.linspace(0,50,501)
Tc.value = data['Tc'].values
m_val.options.IMODE = 4
m_val.solve(disp=False)

plt.figure()
plt.subplot(2,1,1)
plt.title('模型验证')
plt.plot(m_val.time,Tc.value,'b-',label='冷却夹套温度')
plt.ylabel(r'$T_c (K)$')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m_val.time,T.value,'r-',label='反应器温度')
plt.ylabel('T (K)')
plt.xlabel('时间 (min)')
plt.legend()

plt.show()

同理，阶跃测试代码也可以使用独立的Gekko实例，进一步避免状态干扰：

# 新建Gekko实例用于阶跃测试
m_step = GEKKO(remote=False)

# 基于ARX模型构建阶跃测试模型
yc,uc = m_step.arx(p)

# 重命名操作变量与被控变量
Tc = uc[0]
T = yc[0]

# 稳态初始化
m_step.options.IMODE = 1
Tc.value = 300
m_step.solve(disp=False)

# 动态仿真（阶跃测试验证）
m_step.time = np.linspace(0,2,21)
m_step.options.IMODE = 4
Tc.value = np.ones(21)*300
Tc.value[5:] = 305
m_step.solve(disp=False)

plt.figure()
plt.subplot(2,1,1)
plt.title('阶跃测试')
plt.plot(m_step.time,Tc.value,'b-',label='冷却夹套温度')
plt.ylabel(r'$T_c (K)$')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m_step.time,T.value,'r-',label='反应器温度')
plt.ylabel('T (K)')
plt.xlabel('时间 (min)')
plt.legend()

plt.show()

内容的提问来源于stack exchange，提问作者Shawn B