GEKKO ARX Model Forecasting: data arrays must have the same length

user_11313341 0
Viewed 79

I am trying to use sysid to regress an ARX model periodically, then evaluate the predictive ability of that model by simulating with the future inputs and comparing the output with the experimental data. When I try to solve using m.solve() I get the following error: Exception: Data arrays must have the same length, and match time discretization in dynamic problems

The following is an MRE:

X = [[ 0.9,  0.], [ 0.9, 0.],[ 0.9,0.],[ 0.9,0.],[ 0.9, 0.],[ 0.5,0.],[0.5,0.],[0.5,0.],[0.5,0.], [ 0.5, 0.]] # 2 values for inputs at each time step
Y = [20.3, 20.3, 20.2, 20.2, 20.1, 20.1, 20.1, 20.,  19.9, 19.8,] # 1 output at each time step

t = np.linspace(0, 9*300, 10)  # 10 points 5 minutes apart each

na = 1  # output coefficients
nb = 2  # input coefficients
res, p, K = m.sysid(t, X, Y, na, nb, pred='meas')

m.time = t - t[0]
y_, u_ = m.arx(p)
u_[0].value = X[0]
u_[1].value = X[1]

m.options.imode = 4
m.options.nodes = 2

# simulate
m.solve()

I don't want to control, rather apply the experimental values to future timesteps and see how the ARX model extrapolates.

Thanks for your help

gekko time-series python
1 Answers

The problem was with this section and how the data is loaded into value:

u_[0].value = X[:,0]
u_[1].value = X[:,1]
y_[0].value = Y

Try printing X[0] and X[1] to see that they are only the first and second elements of the original list. Converting to a numpy array helps with the slicing.

import numpy as np
from gekko import GEKKO

m = GEKKO()

X = np.array([[ 0.9,  0.], [ 0.9, 0.],\
     [ 0.9,0.],[ 0.9,0.],\
     [ 0.9, 0.],[ 0.5,0.],\
     [0.5,0.],[0.5,0.],\
     [0.5,0.], [ 0.5, 0.]]) # 2 values for inputs at each time step
Y = np.array([20.3, 20.3, 20.2, 20.2,\
     20.1, 20.1, 20.1, 20., \
     19.9, 19.8]) # 1 output at each time step

t = np.linspace(0, 9*300, 10)  # 10 points 5 minutes apart each

na = 1  # output coefficients
nb = 2  # input coefficients
res, p, K = m.sysid(t, X, Y, na, nb, pred='meas')

m.time = t - t[0]
y_, u_ = m.arx(p)
u_[0].value = X[:,0]
u_[1].value = X[:,1]
y_[0].value = Y

print(X[0])
print(X[1])

m.options.imode = 4
m.options.nodes = 2

# simulate
m.solve()

Here is another example with Pandas DataFrames:

system identification

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

# load data and parse into columns
url = 'http://apmonitor.com/do/uploads/Main/tclab_dyn_data2.txt'
data = pd.read_csv(url)
t = data['Time']
u = data[['H1','H2']]
y = data['T1']

# generate time-series model
m = GEKKO(remote=False) # remote=True for MacOS

# system identification
na = 2 # output coefficients
nb = 2 # input coefficients
yp,p,K = m.sysid(t,u,y,na,nb,diaglevel=1)

plt.figure()
plt.subplot(2,1,1)
plt.plot(t,u)
plt.legend([r'$u_0$',r'$u_1$'])
plt.ylabel('MVs')
plt.subplot(2,1,2)
plt.plot(t,y)
plt.plot(t,yp)
plt.legend([r'$y_0$',r'$z_0$'])
plt.ylabel('CVs')
plt.xlabel('Time')
plt.savefig('sysid.png')
plt.show()

We're also working on a package with a Seeq add-on for system identification that runs in Python and Jupyter notebooks.

user_23669410
Related