Python GEKKO：已知延迟的延迟系统的参数估计

使用Gekko中的m.delay()函数，如此Time-Delay example所示.它创建了一个离散状态空间模型来模拟延迟.由于它是混合方程和离散状态空间模型，因此NODES的数字必须是2.结果是-0.2034而不是-0.2，可能是由于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)

还存在初始化问题，因为延迟状态从[0:4]开始未知，因为没有有关状态的过go 信息.使用m.options.TIME_SHIFT=1(默认)重复解决该问题以模拟后退地平线.

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()

还可以估计时间延迟，如外部输入所示，如另一个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()