How to estimate theta value in FOPDT equation using gekko?

Question

I'm trying to use GEKKO to fit to a certain dataset, using the FOPDT Optimization Method to estimate k, tau and theta.

I saw the example using odeint on https://apmonitor.com/pdc/index.php/Main/FirstOrderOptimization and tried to do the same thing with GEKKO, but I can't use the value of theta in the equation.

I saw this question where should the delay call be placed inside a gekko code? and the docs https://apmonitor.com/wiki/index.php/Apps/TimeDelay, but in this case I wanted to estimate the value of theta and not use the initial guess value. I tried to use gekko's delay, but I get an error that it only works if the delay is an int value (and not a gekko FV). I also tried to use time directly in the equation, but I can't figure out how to place x(t-theta) in there, since I can't do that syntax with gekko variables.

import pandas as pd
import numpy as np
from gekko import GEKKO
import plotly.express as px 

data = pd.read_csv('data.csv',sep=',',header=0,index_col=0)

xm1 = data['x']
ym1 = data['y']
xm = xm1.to_numpy()
ym = ym1.to_numpy()

xm_r = len(xm)
tm = np.linspace(0,xm_r-1,xm_r)

m = GEKKO()
m.options.IMODE=5
m.time = tm

k = m.FV()
k.STATUS=1
tau = m.FV()
tau.STATUS=1
theta = m.FV()
theta.STATUS=1

x = m.Param(value=xm)
y = m.CV()
y.FSTATUS = 1
yObj = m.Param(value=ym)

xtheta = m.Var()
m.delay(x,xtheta,theta)

m.Equation(y.dt()==(-y + k * xtheta)/tau)
m.Minimize((y-yObj)**2)
m.options.EV_TYPE=2

m.solve(disp=True)

Answer 1

Here are some strategies for implementing variable time-delay in a model such as when an optimizer adjusts the time delay in a First Order Plus Dead Time (FOPDT) model.

• Create a cubic spline (continuous approximation) of the relationship between time t and the input u. This allows a fractional time delay that is not restricted to an integer multiple of the sample interval.
• Create time as a variable with derivative equal to 1.
• Define tc with an equation tc==time-theta to get the time shifted value. This will lookup the spline uc value that corresponds to this tc value.

You can also fit the FOPDT model to data with Excel or other tools.

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

# load data
url = 'http://apmonitor.com/do/uploads/Main/tclab_siso_data.txt'
data = pd.read_csv(url)
t = data['time'].values
u = data['voltage'].values
y = data['temperature'].values

m = GEKKO(remote=False)
m.time = t; time = m.Var(0); m.Equation(time.dt()==1)

K = m.FV(lb=0,ub=1);      K.STATUS=1
tau = m.FV(lb=1,ub=300);  tau.STATUS=1
theta = m.FV(lb=2,ub=30); theta.STATUS=1

# create cubic spline with t versus u
uc = m.Var(u); tc = m.Var(t); m.Equation(tc==time-theta)
m.cspline(tc,uc,t,u,bound_x=False)

ym = m.Param(y)
yp = m.Var(y); m.Equation(tau*yp.dt()+(yp-y[0])==K*(uc-u[0]))

m.Minimize((yp-ym)**2)

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

print('K: ', K.value[0])
print('tau: ',  tau.value[0])
print('theta: ', theta.value[0])

plt.figure()
plt.subplot(2,1,1)
plt.plot(t,u)
plt.legend([r'$V_1$ (mV)'])
plt.ylabel('MV Voltage (mV)')
plt.subplot(2,1,2)
plt.plot(t,y)
plt.plot(t,yp)
plt.legend([r'$T_{1meas}$',r'$T_{1pred}$'])
plt.ylabel('CV Temp (degF)')
plt.xlabel('Time')
plt.savefig('sysid.png')
plt.show()

K:  0.25489655932
tau:  229.06377617
theta:  2.0

Another way to approach this is to estimate a higher-order ARX model and then determine the statistical significance of the beta terms. Here is an example of using the Gekko sysid function.

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_siso_data.txt'
data = pd.read_csv(url)
t = data['time']
u = data['voltage']
y = data['temperature']

# generate time-series model
m = GEKKO()

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

print('alpha: ', p['a'])
print('beta: ',  p['b'])
print('gamma: ', p['c'])

plt.figure()
plt.subplot(2,1,1)
plt.plot(t,u)
plt.legend([r'$V_1$ (mV)'])
plt.ylabel('MV Voltage (mV)')
plt.subplot(2,1,2)
plt.plot(t,y)
plt.plot(t,yp)
plt.legend([r'$T_{1meas}$',r'$T_{1pred}$'])
plt.ylabel('CV Temp (degF)')
plt.xlabel('Time')
plt.savefig('sysid.png')
plt.show()

With results:

alpha:  [[0.525143  ]
 [0.19284469]
 [0.08177381]
 [0.06152181]
 [0.12918898]]
beta:  [[[-8.51804876e-05]
  [ 5.88425202e-04]
  [ 1.99205676e-03]
  [-2.81456773e-03]
  [ 2.38110003e-03]]]
gamma:  [0.75189199]

The first two beta terms are nearly zero but they can also be left in the model for a higher-order representation of the system.