Run ARX Simulation with initial condition with GEKKO

I want to do a temperature simulation with ARX models and found Gekko as framework for it. I was able to run basically the example of apmonitor and give my simulation an initial condition when using na=1. These are my plots with na=1 and using the first measured value as inital condition in the step test:

Looks fine so far.

However when i increase na to any value higher than 1 lets say 2 my step test now looks like:

The simulation values rise to nonsensical values, but they do in fact start at the correct value for the first value.

My understanding is that with increasing na the simulation takes more past steps into account. How can i start the simulation with multiple past steps as initial state so what i am doing works?

I dont want to use a steady state inital condition because it does not reflect the real starting point of the simulation. I just want to use na real past measurements as starting point.

Here is my code:

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()

I already tried to set to

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

so 2 past steps are available for the simulation, but then I get an exception that the size is not correct.

When I initialize the values with arrays like [20,20,0,0,0,0,0,...] or [20,20,None,None,None,...] the simulation runs without exception but still with the same completely wrong values and also only one 20 is after the simulation in the yc value instead of the set two.

ARX models need to be initialized for dynamic simulation. Add the following steady-state initialization before a dynamic simulation.

# steady-state initialization
m.options.IMODE = 1
# Set input profiles
uc[0].value = data['Q1'].values[0]
uc[1].value = data['Q2'].values[0]
m.solve(disp=False)

Also, use the m.sysid() option shift='init' instead of the default calc if there are good steady-state initial conditions as the first data points in the training set. Here is the complete code:

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, shift='init')

# 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)

# steady-state initialization
m.options.IMODE = 1
# Set input profiles
uc[0].value = data['Q1'].values[0]
uc[1].value = data['Q2'].values[0]
m.solve(disp=False)

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, 1, 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, 1, 2)
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.savefig('results.png',dpi=300)
plt.show()