[Python] Why is output of Gekko ARX model not matching up with SysId regression output?

*Edited on 9/23/24 to add some additional information

• I have run this again after restarting my Jupyter notebook, and it returned the results I expected. However, if I run the full script referenced in the link below, which includes the step test, and then try running the validation procedure again, I get the erroneous results described below. Could his have something to do with how he model is initialized, and residual settings?

I am trying to validate a model of a dynamic process created using the Python Gekko library sysid routine, however, the results I'm seeing don't make sense, or at least are not what I'm expecting.

To demonstrate results I'm getting, I will be using an example from APMonitor tutorials: https://apmonitor.com/do/index.php/Main/NonlinearControl

I run the example and use the output of the sysid routine to generate a coefficient dictionary, p, and a regression output, y. I then create the arx model using p, and define the input (u) to be the same as the input used to generate p and y. After running the solver, I would expect the predicted output yc to be the same a y, but instead I get a something drastically different.

Is there something I am doing wrong?

Here is the code used to generate the original output, and coefficient dictionary, p.

# load data and parse into columns
url = 'http://apmonitor.com/do/uploads/Main/cstr_step_tests.txt'
data = pd.read_csv(url)
print(data.head())

# generate time-series model
t = data['Time']
u = data['Tc']
y = data['T']
m = GEKKO(remote=False)

# system identification
na = 2 # output coefficients
nb = 2 # input coefficients
yp,p,K = m.sysid(t,u,y,na,nb,shift='init',scale=True,objf=100,diaglevel=1)

# plot results of fitting
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')




I then run the code for the step test:

# step test model
yc,uc = m.arx(p)

# rename MV and CV
Tc = uc[0]
T = yc[0]

# steady state initialization
m.options.IMODE = 1
Tc.value = 300
m.solve(disp=False)

# dynamic simulation (step test validation)
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('Step Test')
plt.plot(m.time,Tc.value,'b-',label='Cooling Jacket')
plt.ylabel(r'$T_c (K)$')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m.time,T.value,'r-',label='Reactor')
plt.ylabel('T (K)')
plt.xlabel('Time (min)')
plt.legend()

plt.show()


And finally, I run the subsequent code used to attempt a validation of the coefficients with the original input:

# step test model
yc,uc = m.arx(p)

# rename MV and CV
Tc = uc[0]
T = yc[0]

# steady state initialization
m.options.IMODE = 1
Tc.value = 300
m.solve(disp=False)

# dynamic simulation
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('Step Test')
plt.plot(m.time,Tc.value,'b-',label='Cooling Jacket')
plt.ylabel(r'$T_c (K)$')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m.time,T.value,'r-',label='Reactor')
plt.ylabel('T (K)')
plt.xlabel('Time (min)')
plt.legend()

plt.show()




As you can see, the output of the "validation" doesn't match up with the regression output from the system identification step.

If I skip the step test code, then I get an output that matches the original prediction. Why does running the step test cause the results to become erroneous?

Continue reading...