Gekko ARX model crash possible to load last iteration

Question

I am pretty new to model predictive controls modeling with Gekko and in general.

I have created an ARX MPC in Gekko, which is working great. I, however, noticed that in the first 50-80 iterations, the results are well.. disappointing. However, after the first iterations, I get good results (I guess the ARX algorithm is at play here or possible BIAS?). Now my problem is that the model might crash after some time, and I have to redo the 50-80 iteration to get good results again, is there a way to "save" the last calculated model and use that when rebooting the calculations?

Answer 1

The issue that you are likely encountering is that the "prior" values have not yet been initialized. Try solving once with a steady-state initialization as shown in the example MPC application with the TCLab that is the final source block on for TCLab F.

m.options.IMODE=1
m.solve()

You can then switch to control or simulation mode:

# set up MPC
m.options.IMODE   = 6 # MPC
m.time=np.linspace(0,120,61)

Background information on using ARX models

Identification of the ARX model and prediction or control with the ARX model are two separate applications.

Identify ARX Model

The m.sysid() function to identify an ARX model does not save an archive but does return the model as output arguments:

yp,p,K = m.sysid(t,u,y,na,nb,pred='meas')

The model is returned as p.

# see https://apmonitor.com/wiki/index.php/Apps/ARXTimeSeries
from gekko import GEKKO
import numpy as np
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']
y = data['T1']

m = GEKKO(remote=False)

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

plt.figure()
plt.subplot(2,1,1)
plt.plot(t,u,label=r'$Heater_1$')
plt.legend()
plt.ylabel('Heater')
plt.subplot(2,1,2)
plt.plot(t,y)
plt.plot(t,yp)
plt.legend([r'$T_{meas}$',r'$T_{pred}$'])
plt.ylabel('Temperature (°C)')
plt.xlabel('Time (sec)')
plt.show()

Predict with ARX Model

Below is an example of prediction with the ARX model.

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

na = 2 # Number of A coefficients
nb = 1 # Number of B coefficients
ny = 2 # Number of outputs
nu = 2 # Number of inputs

# A (na x ny)
A = np.array([[0.36788,0.36788],\
              [0.223,-0.136]]) 
# B (ny x (nb x nu))
B1 = np.array([0.63212,0.18964]).T
B2 = np.array([0.31606,1.26420]).T
B = np.array([[B1],[B2]])

C = np.array([0,0])

# create parameter dictionary
# parameter dictionary p['a'], p['b'], p['c']
# a (coefficients for a polynomial, na x ny)
# b (coefficients for b polynomial, ny x (nb x nu))
# c (coefficients for output bias, ny)
p = {'a':A,'b':B,'c':C}

# Create GEKKO model
m = GEKKO(remote=False)

# Build GEKKO ARX model
y,u = m.arx(p)

# load inputs
tf = 20 # final time
u1 = np.zeros(tf+1)
u2 = u1.copy()
u1[5:] = 3.0
u2[10:] = 5.0
u[0].value = u1
u[1].value = u2

# customize names
mv1 = u[0]
mv2 = u[1]
cv1 = y[0]
cv2 = y[1]

# options
m.time = np.linspace(0,tf,tf+1)
m.options.imode = 4
m.options.nodes = 2

# simulate
m.solve()

m.open_folder()

plt.figure(1)
plt.subplot(2,1,1)
plt.plot(m.time,mv1.value,'r-',label=r'$MV_1$')
plt.plot(m.time,mv2.value,'b--',label=r'$MV_2$')
plt.ylabel('MV')
plt.legend(loc='best')
plt.subplot(2,1,2)
plt.plot(m.time,cv1.value,'r:',label=r'$CV_1$')
plt.plot(m.time,cv2.value,'b.-',label=r'$CV_2$')
plt.ylabel('CV')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.show()

The model is saved in the m.path folder that can be viewed with m.open_folder(). Set m = GEKKO(remote=False) to calculate locally and observe all of the files that are used to generate the model and the solution.