2025, Dec 02 21:01
Стабильная симуляция ARX‑моделей в Gekko при na>1: правильная инициализация
Разбираем, почему ARX‑модели в Gekko расходятся при na>1, и показываем, как стабилизировать симуляцию: инициализация IMODE=1 и sysid с параметром shift='init'.
При моделировании температурной динамики с ARX‑моделями в Gekko нередко удаётся получить приличные результаты при na=1, а стоит увеличить порядок (na>1) — симуляция «взрывается» или уходит в неправдоподобные значения. Парадокс в том, что первая сымитированная точка совпадает с данными, а дальше траектория расходится. Дело не в порядке как таковом, а в том, как модель инициализируется для динамического расчёта.
Минимальный пример, который воспроизводит проблему
Ниже приведён скрипт, который выполняет идентификацию на наборе данных ступенчатого теста TCLab, а затем сразу запускает имитацию ARX‑модели с IMODE=4, задавая в качестве начального условия лишь первое значение выхода. Для na=1 это работает, но при более высоких порядках выдаёт бессмысленные результаты.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from gekko import GEKKO
# Загрузка данных
src = 'http://apmonitor.com/dde/uploads/Main/tclab_step_test.txt'
df_raw = pd.read_csv(src)
t_idx = df_raw['Time']
u_mat = df_raw[['Q1', 'Q2']]
y_mat = df_raw[['T1', 'T2']]
# Идентификация системы
gm = GEKKO(remote=True)
na_ord, nb_ord = 2, 6
y_hat, theta, kval = gm.sysid(t_idx, u_mat, y_mat, na_ord, nb_ord, diaglevel=1, scale=False)
# Графики входов и выходов с предсказаниями
plt.figure(figsize=(8, 5))
plt.subplot(2, 1, 1)
plt.plot(t_idx, u_mat)
plt.legend([r'$Q_1$', r'$Q_2$'])
plt.ylabel('Inputs')
plt.subplot(2, 1, 2)
plt.plot(t_idx, y_mat)
plt.plot(t_idx, y_hat)
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‑модель для моделирования
y_comp, u_comp = gm.arx(theta)
gm.time = np.linspace(0, 600, 601)
gm.options.IMODE = 4
# Задать профили входов
u_comp[0].value = df_raw['Q1']
u_comp[1].value = df_raw['Q2']
# Задать начальные значения выходов
y_comp[0].value = df_raw['T1'][0]
y_comp[1].value = df_raw['T2'][0]
# Решение и построение графиков результатов моделирования
gm.solve(disp=False)
plt.figure(figsize=(10, 6))
plt.subplot(2, 2, 1)
plt.title('Step Test')
plt.plot(gm.time, u_comp[0].value, 'b-', label=r'$Q_1$')
plt.plot(gm.time, u_comp[1].value, 'r-', label=r'$Q_2$')
plt.ylabel('Heater (%)')
plt.legend()
plt.subplot(2, 2, 3)
plt.plot(gm.time, y_comp[0].value, 'b--', label=r'$T_1$')
plt.plot(gm.time, y_comp[1].value, 'r--', label=r'$T_2$')
plt.plot(t_idx, y_mat, 'k-', label='Measured')
plt.ylabel('Temperature (K)')
plt.xlabel('Time (sec)')
plt.legend()
plt.tight_layout()
plt.show()
Почему симуляция уходит вразнос
ARX‑моделям нужна инициализация для динамического запуска. Без явной инициализации внутренняя история, необходимая авторегрессионной части, не находится в согласованном установившемся состоянии, поэтому у динамики более высокого порядка нет корректной стартовой точки. Простого присвоения первого образца выхода в качестве начального значения недостаточно — это не заполняет историю, которую ожидает na>1. Попытка передать несколько начальных значений выходов короткими массивами тоже не спасает, потому что перед переходом к динамическому режиму модель нужно корректно вывести в установившееся состояние.
Решение и рабочий скрипт
Нужны два шага, чтобы симуляция стала устойчивой. Во‑первых, выполнить инициализацию в установившемся режиме с IMODE=1, используя первые значения входов. Во‑вторых, если начальные точки данных действительно соответствуют установившемуся состоянию, вызвать sysid с параметром shift='init'. В скрипте ниже оба шага применены.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from gekko import GEKKO
# Загрузка данных
src = 'http://apmonitor.com/dde/uploads/Main/tclab_step_test.txt'
df_raw = pd.read_csv(src)
t_idx = df_raw['Time']
u_mat = df_raw[['Q1', 'Q2']]
y_mat = df_raw[['T1', 'T2']]
# Идентификация системы
gm = GEKKO(remote=True)
na_ord, nb_ord = 2, 6
y_hat, theta, kval = gm.sysid(t_idx, u_mat, y_mat, na_ord, nb_ord, diaglevel=1, scale=False, shift='init')
# Графики входов и выходов с предсказаниями
plt.figure(figsize=(8, 5))
plt.subplot(2, 1, 1)
plt.plot(t_idx, u_mat)
plt.legend([r'$Q_1$', r'$Q_2$'])
plt.ylabel('Inputs')
plt.subplot(2, 1, 2)
plt.plot(t_idx, y_mat)
plt.plot(t_idx, y_hat)
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‑модель для моделирования
y_comp, u_comp = gm.arx(theta)
# Инициализация в установившемся режиме
gm.options.IMODE = 1
u_comp[0].value = df_raw['Q1'].values[0]
u_comp[1].value = df_raw['Q2'].values[0]
gm.solve(disp=False)
# Динамическое моделирование
gm.time = np.linspace(0, 600, 601)
gm.options.IMODE = 4
u_comp[0].value = df_raw['Q1']
u_comp[1].value = df_raw['Q2']
y_comp[0].value = df_raw['T1'][0]
y_comp[1].value = df_raw['T2'][0]
gm.solve(disp=False)
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.title('Step Test')
plt.plot(gm.time, u_comp[0].value, 'b-', label=r'$Q_1$')
plt.plot(gm.time, u_comp[1].value, 'r-', label=r'$Q_2$')
plt.ylabel('Heater (%)')
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(gm.time, y_comp[0].value, 'b--', label=r'$T_1$')
plt.plot(gm.time, y_comp[1].value, 'r--', label=r'$T_2$')
plt.plot(t_idx, y_mat, 'k-', label='Measured')
plt.ylabel('Temperature (K)')
plt.xlabel('Time (sec)')
plt.legend()
plt.tight_layout()
plt.savefig('results.png', dpi=300)
plt.show()
Почему это важно
Выбор порядка в ARX‑модели напрямую определяет, сколько прошлой информации она использует. Если старт не из согласованного установившегося состояния, члены истории более высокого порядка сразу искажают траекторию. Корректная инициализация выравнивает динамическую историю с идентифицированной моделью, поэтому ступенчатые тесты и другие сценарии ведут себя адекватно при разных значениях na.
Итоги
Если симуляция Gekko ARX «взрывается» при na>1, перед динамическим запуском проинициализируйте модель в установившемся режиме и, когда первые отсчёты действительно стационарны, вызывайте sysid с shift='init'. Так стартовые условия согласуются с идентифицированной динамикой: можно безопасно повышать na без потери устойчивости и получать ступенчатые тесты, отражающие ваши реальные начальные условия.