以下是在模型中实现可变时间延迟的一些策略,例如当优化器调整一阶加死区时间 (FOPDT) 模型中的时间延迟时。
- 创建时间之间关系的三次样条(连续近似)
t
和输入u
。这允许不限于采样间隔的整数倍的分数时间延迟。
- Create
time
作为导数等于 1 的变量。
- Define
tc
用一个方程tc==time-theta
以获得时移值。这将查找样条曲线uc
与此对应的值tc
value.
你也可以拟合 FOPDT 模型使用 Excel 或其他工具获取数据。
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
解决这个问题的另一种方法是估计高阶ARX模型然后确定统计显着性beta
条款。这是使用 Gekko 的示例sysid
功能。
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
项几乎为零,但也可以将它们保留在模型中以获得系统的高阶表示。