GEKKO中ARX模型na>1时多初始状态设置问题求助

问题场景 #

使用GEKKO框架基于ARX模型进行温度仿真：

• na=1时，可正常运行示例并设置初始条件，仿真结果符合预期；
• na增大至2（或更大）时，仿真结果出现无意义数值异常，仅首个值匹配初始值。

需求：用真实历史测量值作为na个初始状态（而非稳态初始条件），但设置yc数组初始值时要么出现尺寸错误，要么结果仍异常。

原代码 #

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

# Load data
url = 'http://apmonitor.com/dde/uploads/Main/tclab_step_test.txt'
data = pd.read_csv(url)
t = data['Time']
u = data[['Q1', 'Q2']]
y = data[['T1', 'T2']]

# System identification
m = GEKKO(remote=True)
na, nb = 2, 6
yp, p, K = m.sysid(t, u, y, na, nb, diaglevel=1, scale=False)

# Plot input and output data with predictions
plt.figure(figsize=(8, 5))
plt.subplot(2, 1, 1)
plt.plot(t, u)
plt.legend([r'$Q_1$', r'$Q_2$'])
plt.ylabel('Inputs')
plt.subplot(2, 1, 2)
plt.plot(t, y)
plt.plot(t, yp)
plt.legend([r'$T_{1,meas}$', r'$T_{2,meas}$',
           r'$T_{1,pred}$', r'$T_{2,pred}$'])
plt.ylabel('Outputs')
plt.xlabel('Time')

# ARX model for simulation
yc, uc = m.arx(p)
m.time = np.linspace(0, 600, 601)
m.options.IMODE = 4

# Set input profiles
uc[0].value = data['Q1']
uc[1].value = data['Q2']

# Set initial output values
yc[0].value = data['T1'][0]
yc[1].value = data['T2'][0]

# Solve and plot simulation results
m.solve(disp=False)

plt.figure(figsize=(10, 6))
plt.subplot(2, 2, 1)
plt.title('Step Test')
plt.plot(m.time, uc[0].value, 'b-', label=r'$Q_1$')
plt.plot(m.time, uc[1].value, 'r-', label=r'$Q_2$')
plt.ylabel('Heater (%)')
plt.legend()

plt.subplot(2, 2, 3)
plt.plot(m.time, yc[0].value, 'b--', label=r'$T_1$')
plt.plot(m.time, yc[1].value, 'r--', label=r'$T_2$')
plt.plot(t, y, 'k-', label='Measured')
plt.ylabel('Temperature (K)')
plt.xlabel('Time (sec)')
plt.legend()

plt.tight_layout()
plt.show()

尝试过的错误方法 #

• 直接设置短数组：

  yc[0].value = [20, 20]
  yc[1].value = [20, 20]
  

  出现尺寸错误，因为yc的value长度必须与m.time的长度一致。

• 设置长数组但填充非空值：

  yc[0].value = [20,20,0,0,...]  # 或 [20,20,None,...]
  

  无报错但结果异常，yc仅保留首个20值。

正确解决方案 #

ARX模型当na>1时，需要提供na个历史测量值作为初始条件，剩余时间步用None让模型自动计算。具体操作：

1. 获取前na个历史输出测量值；
2. 构造与m.time长度一致的数组，前na位为历史值，后续为None；
3. 赋值给yc的value属性。

修正后的初始值设置代码 #

# 设置初始输出值：前na个为历史测量值，后续由模型计算
n_steps = len(m.time)
# T1的初始历史值：取前na个测量值
yc0_init = list(data['T1'][:na]) + [None]*(n_steps - na)
yc[0].value = yc0_init
# T2的初始历史值：取前na个测量值
yc1_init = list(data['T2'][:na]) + [None]*(n_steps - na)
yc[1].value = yc1_init

完整修正后的仿真代码片段 #

替换原代码中Set initial output values部分为上述代码即可。

原理说明 #

• ARX模型的递归特性：当na=2时，模型需要t=0和t=-1时刻的输出值来计算t=1时刻的预测值；
• GEKKO变量的value数组长度必须与m.time匹配，None表示该位置的值由模型求解得到，而非手动赋值；
• 用真实历史测量值填充前na位，确保模型从正确的初始历史状态开始仿真。

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