GEKKO: Parameter estimation on a delayed system with known delay

Here is a simple delayed system, or say a DDE

I believe the parameter k can be estimated with direct collocation method when τ is given. And the problem is how to do this with GEKKO.

Here is what I have done:

Firstly, the system is simulated with jitcdde package

from jitcdde import t, y, jitcdde
import numpy as np

# the constants in the equation
k = -0.2
tau = 0.5

# the equation
f = [    
    k*y(0, t-tau)
    ]

# initialising the integrator
DDE = jitcdde(f)

# enter initial conditions
DDE.constant_past([1.0])

# short pre-integration to take care of discontinuities
DDE.step_on_discontinuities()

# create timescale
stoptime = 10.5
numpoints = 100
times = np.arange(DDE.t, stoptime, 1/10)

# integrating
data = []
for time in times:
    data.append( DDE.integrate(time) ) # The warning is due to small sample step, no worry

Then, do the estimation

# The sample is taken every 0.1s and the delay is 0.5s
# so ydelayed should be taken 5 steps before ydata
tdata = times[5:]
ydata =  np.array(data)[5:]
ydelayed = np.array(data)[0:-5]

m = GEKKO(remote=False)
m.time = tdata
x = m.CV(value=ydata); x.FSTATUS = 1 # fit to measurement
xdelayed = m.CV(value=ydelayed); x.FSTATUS = 1  # fit to measurement
k = m.FV(); k.STATUS = 1               # adjustable parameter
m.Equation(x.dt()== k * xdelayed)            # differential equation

m.options.IMODE = 5   # dynamic estimation
m.options.NODES = 5   # collocation nodes
m.solve(disp=False)   # display solver output
k = k.value[0]

which show k=-0.03253 but the true value is -0.2. How to fix the identification?

Use the m.delay() function in Gekko as shown in this Time-Delay example. It creates a discrete state-space model to simulate the delay. Because it is a mixed differential equation and discrete state-space model, the number of NODES must be 2. The result is -0.2034 instead of -0.2, likely due to the discretization error with NODES=2.

from gekko import GEKKO
# The sample is taken every 0.1s and the delay is 0.5s
# so ydelayed should be taken 5 steps before ydata
tdata = times[5:]
ydata =  np.array(data)[5:]
ydelayed = np.array(data)[0:-5]

m = GEKKO(remote=False)
m.time = tdata
xdata = m.Param(ydata)
x = m.Var(value=ydata)
xdelayed = m.Var(value=ydelayed)
m.delay(x,xdelayed,5)

m.Minimize((x-xdata)**2)

k = m.FV(); k.STATUS = 1               # adjustable parameter
m.Equation(x.dt()== k * xdelayed)            # differential equation

m.options.IMODE = 5   # dynamic estimation
m.options.NODES = 2   # collocation nodes

for i in range(5):
  xdata.value = ydata
  m.solve(disp=False)   # display solver output
k = k.value[0]
print(k)

There is also an initialization issue because the delay states are unknown from [0:4] because there is no past information on the states. The problem is solved repeatedly with m.options.TIME_SHIFT=1 (default) to simulate the receding horizon.

import matplotlib.pyplot as plt
plt.figure(figsize=(6,3))
plt.plot(m.time,x.value,'ro',label='Meas')
plt.plot(m.time,xdata,'b--',label='Pred')
plt.legend(); plt.grid()
plt.show()

It is also possible to estimate the time-delay as shown with exogenous inputs as shown with this self-contained example of another first-order linear system.

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

# Import CSV data file
# Column 1 = time (t)
# Column 2 = input (u)
# Column 3 = output (yp)
url = 'http://apmonitor.com/pdc/uploads/Main/data_fopdt.txt'
data = pd.read_csv(url)
t = data['time'].values - data['time'].values[0]
u = data['u'].values
y = data['y'].values

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

K = m.FV(2,lb=0,ub=10);      K.STATUS=1
tau = m.FV(3,lb=1,ub=200);  tau.STATUS=1
theta_ub = 30 # upper bound to dead-time
theta = m.FV(0,lb=0,ub=theta_ub); theta.STATUS=1

# add extrapolation points
td = np.concatenate((np.linspace(-theta_ub,min(t)-1e-5,5),t))
ud = np.concatenate((u[0]*np.ones(5),u))
# create cubic spline with t versus u
uc = m.Var(u); tc = m.Var(t); m.Equation(tc==time-theta)
m.cspline(tc,uc,td,ud,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('Kp: ', K.value[0])
print('taup: ',  tau.value[0])
print('thetap: ', theta.value[0])

# plot results
plt.figure()
plt.subplot(2,1,1)
plt.plot(t,y,'k.-',lw=2,label='Process Data')
plt.plot(t,yp.value,'r--',lw=2,label='Optimized FOPDT')
plt.ylabel('Output')
plt.legend()
plt.subplot(2,1,2)
plt.plot(t,u,'b.-',lw=2,label='u')
plt.legend()
plt.ylabel('Input')
plt.show()