提高 Python gekko 中非线性回归的质量

Question

我正在尝试使用 python GEKKO 非线性回归工具来执行使用阶跃响应的二阶过阻尼系统的系统识别。

我的代码如下：

m = GEKKO()
m_input = m.Param(value=input)
m_time=m.Param(value=time)
m_T1 = m.FV(value=initT1, lb=T1bounds[0], ub=T1bounds[1])
m_T1.STATUS = 1
m_k = m.FV(value=initk,lb=100)
m_k.STATUS = 1

m_T2 = m.FV(value=initT2, lb=T2bounds[0], ub=T2bounds[1])
m_T2.STATUS = 1


m_output = m.CV(value=output)
m_output.FSTATUS=1

m.Equation(m_output==(m_k/(m_T1+m_T2))*(1+((m_T1/(m_T2-m_T1))*m.exp(-m_time/m_T2))-((m_T2/(m_T2-m_T1))*m.exp(-m_time/m_T1)))*m_input)
m.options.IMODE = 2
m.options.MAX_ITER = 10000
m.options.OTOL = 1e-8
m.options.RTOL = 1e-8
m.solve(disp=True)

结果并不乐观。似乎优化器似乎陷入了目标函数的局部最小值，从而使目标函数过高

求解器的输出为：

The final value of the objective function is    160453.282142838     
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :    7.60390000000189      sec
 Objective      :    160453.282605857     
 Successful solution
 ---------------------------------------------------

我可以做些什么来提高合身的质量？我可以对目标函数值进行限制吗？

Answer 1

尝试将方程用于欠阻尼二阶系统，而不是过阻尼二阶系统。在Process Dynamics and Control website 上有更多关于二阶系统显式解决方案的信息。更好的方法是不使用显式解决方案，其中必须预先确定它是欠阻尼、临界阻尼还是过阻尼。如果使用原始的二阶微分方程，那么求解器可以决定它是欠阻尼（振荡过冲）还是过阻尼。这是拟合二阶系统的示例代码。

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 (y)
url = 'http://apmonitor.com/pdc/uploads/Main/data_sopdt.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
zeta = m.FV(1,lb=0.1,ub=3);  zeta.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); xp = m.Var(y)
m.Equation(xp==yp.dt())
m.Equation((tau**2)*xp.dt()+2*zeta*tau*yp.dt()+yp==K*(uc-u[0]))

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

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

print('Kp: ', K.value[0])
print('taup: ',  tau.value[0])
print('thetap: ', theta.value[0])
print('zetap: ', zeta.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 SOPDT')
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()

请用你的数据试试这个，看看它是否能提供更好的解决方案。