Search code examples
pythonloopsregressionautoregressive-modelsgekko

Adaptive modelling using GEKKO sysid

I have 100 points of data that I'm trying to describe using sysid in GEKKO. At some point ( t = 50 in this case) the data changes significantly and the prediction is no longer accurate. I'm trying to include an if statement that evaluates actual vs predicted and generates a new model (new yp ) if the prediction is x times bigger than the model. Here's my sample code. The loop continues evaluating yp at each point in time but it should be evaluating yp_new now. 

from gekko import GEKKO
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
# load data and parse into columns
t = np.linspace(0,1,101)
u = np.linspace(0,1,101)
y = np.zeros(len(t))
y[:50] = np.sin(u[:50])
y[50:] = np.exp(u[50:]/500)
# generate time-series model
m = GEKKO(remote=False)
# system identification
na = 2 # output coefficients
nb = 2 # input coefficients
yp,p,K = m.sysid(t,u,y,na,nb,diaglevel=1)
print(yp)
for i in range(len(t)):
    difference = np.abs((yp[i]-y[i])/max(0.01,y[i]))
    if difference>=0.2:   #If the difference is >20%
        yp_new,p_new,K = m.sysid(t,u,y,na,nb,diaglevel=0)
        print('Recalculating at i  = ' + str(i))
print(yp_new)
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.plot(t,y)
plt.plot(t,yp_new)
plt.show()
Solution

  • You would need to update yp to be the new yp value yp = yp_new or else just return yp when doing the next system identification. However, the data that you are using to redo the sysid calculation is the same as the data that you used previously so there is no change in the model prediction. Are you trying to update the time-series model only with the most recent data such as yp_new,p_new,K = m.sysid(t[i:],u[i:],y[i:],na,nb)?

    System Identification

    The model currently updates on cycles that do not agree with the original time-series model prediction.

    Recalculating at i  = 10
Recalculating at i  = 11
Recalculating at i  = 12
Recalculating at i  = 13
Recalculating at i  = 14
Recalculating at i  = 15
Recalculating at i  = 16
Recalculating at i  = 17
Recalculating at i  = 18
Recalculating at i  = 19
Recalculating at i  = 20
Recalculating at i  = 21
Recalculating at i  = 22
Recalculating at i  = 23
Recalculating at i  = 24
Recalculating at i  = 25
Recalculating at i  = 40
Recalculating at i  = 41
Recalculating at i  = 42
Recalculating at i  = 43
Recalculating at i  = 44
Recalculating at i  = 45
Recalculating at i  = 46
Recalculating at i  = 47
Recalculating at i  = 48
Recalculating at i  = 49
Recalculating at i  = 50
Recalculating at i  = 51
Recalculating at i  = 52

    If you include just the most recent data then it only recalculates at cycles 10, 42, 46, and 48.

    from gekko import GEKKO
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
# load data and parse into columns
t = np.linspace(0,1,101)
u = np.linspace(0,1,101)
y = np.zeros(len(t))
y[:50] = np.sin(u[:50])
y[50:] = np.exp(u[50:]/500)
# generate time-series model
m = GEKKO(remote=False)
# system identification
na = 2 # output coefficients
nb = 2 # input coefficients
yp,p,K = m.sysid(t,u,y,na,nb)
yp_init = yp.copy()
print(yp)
j = 0
for i in range(len(t)):
    difference = np.abs((yp[i-j]-y[i])/max(0.01,y[i]))
    if difference>=0.2:   #If the difference is >20%
        j = i # get cycle where the last update occurred
        yp,p,K = m.sysid(t[i:],u[i:],y[i:],na,nb)
        print('Recalculating at i  = ' + str(i))
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_init)
plt.plot(t,y)
plt.plot(t[j:],yp)
plt.show()

    Updated sysid