如何用gekko估计FOPDT方程中的θ值？

Question

我试图用GEKKO来拟合某个数据集，用FOPDT优化方法来估计k，tau和θ

我在https://apmonitor.com/pdc/index.php/Main/FirstOrderOptimization上看到了使用odeint的例子，并尝试对GEKKO做同样的事情，但我不能在方程中使用θ的值

我看到了这个问题{a2}和文档{a3}，但在这个例子中，我想估计θ的值，而不是使用初始猜测值。我尝试使用gekko的延迟，但我得到一个错误，即只有当延迟是int值（而不是gekko FV）时，它才起作用。 我也试着在方程中直接使用时间，但我不知道如何将x（tθ）放在那里，因为我不能用gekko变量来做语法

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

下面是一些在模型中实现可变时延的策略，例如优化器在一阶加死时间（FOPDT）模型中调整时延时

• 创建时间t和输入u之间关系的三次样条曲线（连续近似）。这允许不限于采样间隔整数倍的分数延时
• 将time创建为导数等于1的变量
• 用方程tc==time-theta定义tc以获得时移值。这将查找与此tc值对应的样条曲线uc值

您还可以使用Excel或其他工具fit the FOPDT model转换数据

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

另一种方法是估计higher-order ARX model，然后确定beta项的统计显著性。下面是一个使用Gekkosysid函数的示例

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

结果如下：

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]

前两个beta项几乎为零，但它们也可以留在模型中，用于系统的高阶表示