GEKKO中的FOPDT方程组-使用多个输入

如何解决GEKKO中的FOPDT方程组-使用多个输入

我想使用两个或多个输入来创建变量的更精确估计。我已经仅使用一个输入和一个FOPDT方程进行了估算，但是当我尝试添加一个以上的输入以及相应的k，tau和theta以及另一个方程时，出现“未找到解决方案”错误。我可以这样创建方程组吗？

有关以下求解器输出的更多详细信息。即使我添加了比方程式更多的变量来使用多个输入，这仍然使我的问题具有负的自由度。

from gekko import GEKKO
import numpy as np
import pandas as pd
import plotly.express as px 

d19jc = d19jc.dropna()

d19jcslice = d19jc.loc['2019-10-22 05:30:00':'2019-10-22 09:30:00'] #jc22

d19jcslice.index = pd.to_datetime(d19jcslice.index)
d19jcsliceGroupMin = d19jcslice.groupby(pd.Grouper(freq='T')).mean()
data = d19jcsliceGroupMin.dropna()

xdf1 = data['Cond_PID_SP']
xdf2 = data['Front_PID_SP']
ydf1 = data['Cond_Center_Top_TC']

xms1 = pd.Series(xdf1)
xm1 = np.array(xms1)
xms2 = pd.Series(xdf2)
xm2 = np.array(xms2)
yms = pd.Series(ydf1)
ym = np.array(yms)

xm_r = len(xm1)
tm = np.linspace(0,xm_r-1,xm_r)

m = GEKKO()
m.options.IMODE=5
m.time = tm; time = m.Var(0); m.Equation(time.dt()==1)

k1 = m.FV(lb=0.1,ub=5); k1.STATUS=1
tau1 = m.FV(lb=1,ub=300); tau1.STATUS=1
theta1 = m.FV(lb=0,ub=30); theta1.STATUS=1
k2 = m.FV(lb=0.1,ub=5); k2.STATUS=1
tau2 = m.FV(lb=1,ub=300); tau2.STATUS=1
theta2 = m.FV(lb=0,ub=30); theta2.STATUS=1

# create cubic spline with t versus u
uc1 = m.Var(xm1); tc1 = m.Var(tm); m.Equation(tc1==time-theta1)
m.cspline(tc1,uc1,tm,xm1,bound_x=False)
# create cubic spline with t versus u
uc2 = m.Var(xm2); tc2 = m.Var(tm); m.Equation(tc2==time-theta2)
m.cspline(tc2,uc2,xm2,bound_x=False)

x1 = m.Param(value=xm1)
x2 = m.Param(value=xm2)
y = m.Var(value=ym)
yObj = m.Param(value=ym)

m.Equation(tau1*y.dt()+(y-ym[0])==k1 * (uc1-xm1[0]))
m.Equation(tau2*y.dt()+(y-ym[0])==k2 * (uc2-xm2[0]))
m.Minimize((y-yObj)**2)
m.options.EV_TYPE=2

print('solve start')
m.solve(disp=True)

print('k1: ',k1.value[0])
print('tau1: ',tau1.value[0])
print('theta1: ',theta1.value[0])
print('k2: ',k2.value[0])
print('tau2: ',tau2.value[0])
print('theta2: ',theta2.value[0])

df_plot = pd.DataFrame({'DateTime' : data.index,'Cond_Center_Top_TC' : np.array(yObj),'Fit Cond_Center_Top_TC - Train' : np.array(y),figGekko = px.line(df_plot,x='DateTime',y=['Cond_Center_Top_TC','Fit Cond_Center_Top_TC - Train'],labels={"value": "degrees Celsius"},title = "(Cond_PID_SP & Front_PID_SP) -> Cond_Center_Top_TC JC only - Train")
figGekko.update_layout(legend_title_text='')
figGekko.show()

这是代码

$ repo --version
repo version v2.4.1-9-g910dfe8
       (from C:\Users\KNG1LUD\bin\git-repo\.git)
repo launcher version 2.4
       (from C:\Users\KNG1LUD\bin\git-repo\repo)
       (currently at 2.4.1-9-g910dfe8)
repo User-Agent git-repo/2.4.1-9-g910dfe8 (Win32) git/2.29.2.windows.1 Python/3.6.0
git 2.29.2.windows.1
git User-Agent git/2.29.2.windows.1 (Win32) git-repo/2.4.1-9-g910dfe8
Python 3.6.0 (v3.6.0:41df79263a11,Dec 23 2016,08:06:12) [MSC v.1900 64 bit (AMD64)]
OS Windows 10 (10.0.17763)
cpu AMD64 (Intel64 Family 6 Model 142 Stepping 10,GenuineIntel)

解决方法

您目前只有一个变量和两个方程式。

y = m.Var(value=ym)
yObj = m.Param(value=ym)

m.Equation(tau1*y.dt()+(y-ym[0])==k1 * (uc1-xm1[0]))
m.Equation(tau2*y.dt()+(y-ym[0])==k2 * (uc2-xm2[0]))

您需要为两个方程式创建两个单独的变量y1和y2。

y1 = m.Var(value=ym)
y2 = m.Var(value=ym)
yObj = m.Param(value=ym)

m.Equation(tau1*y1.dt()+(y1-ym[0])==k1 * (uc1-xm1[0]))
m.Equation(tau2*y2.dt()+(y2-ym[0])==k2 * (uc2-xm2[0]))

如果每个数据都有不同的数据，则可能还需要为ym1和ym2创建两个单独的测量向量。

编辑：多输入单输出（MISO）系统

对于MISO系统，您需要按照Process Dynamics and Control course所示调整方程式。

y = m.Var(value=ym)
yObj = m.Param(value=ym)

m.Equation(tau*y.dt()+(y-ym[0])==k1 * (uc1-xm1[0]) + k2 * (uc2-xm2[0]))

这种形式只有一个时间常数。如果它们需要不同的时间常数，则还可以将两个输出与以下内容相加：

y1 = m.Var(value=ym[0])
y2 = m.Var(value=ym[0])
y  = m.Var(value=ym)
yObj = m.Param(value=ym)

m.Equation(tau1*y1.dt()+(y1-ym[0])==k1 * (uc1-xm1[0]))
m.Equation(tau2*y2.dt()+(y2-ym[0])==k2 * (uc2-xm2[0]))
m.Equation(y==y1+y2)

在这种情况下，y是具有数据的变量，而y1和y2是未测量状态。

,

请参阅以下示例代码，了解具有两个输入和一个输出的FOPDT模型。两种传递函数模型都具有不同的FOPDT参数，包括不同的死区时间。
它仍然处于动态仿真模式（IMODE = 4），但是您可以从此修改开始一点，之后再更改为动态估计模式（IMODE = 5）和MPC模式（IMODE = 6）。

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

tf = 100 
npt = 101 
t = np.linspace(0,tf,npt)
u1 = np.zeros(npt)
u2 = np.zeros(npt)
u1[10:] = 5
u2[40:] = -5

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

K1 = m.FV(1,lb=0,ub=1);      K1.STATUS=1
tau1 = m.FV(5,lb=1,ub=300);  tau1.STATUS=1
theta1 = m.FV(10,lb=2,ub=30); theta1.STATUS=1

K2 = m.FV(0.5,ub=1);      K2.STATUS=1
tau2 = m.FV(10,ub=300);  tau2.STATUS=1
theta2 = m.FV(20,ub=30); theta2.STATUS=1

uc1 = m.Var(u1) 
uc2 = m.Var(u2) 
tc1 = m.Var(t) 
tc2 = m.Var(t)
m.Equation(tc1==time-theta1)
m.Equation(tc2==time-theta2)
m.cspline(tc1,uc1,t,u1,bound_x=False)
m.cspline(tc2,uc2,u2,bound_x=False)

yp = m.Var() 
yp1 = m.Var()
yp2 = m.Var()
m.Equation(yp1.dt() == -yp1/tau1 + K1*uc1/tau1) 
m.Equation(yp2.dt() == -yp2/tau2 + K2*uc2/tau2)
m.Equation(yp == yp1 + yp2)

m.options.IMODE=4
m.solve()

print('K1: ',K1.value[0])
print('tau1: ',tau1.value[0])
print('theta1: ',theta1.value[0])
print('')

print('K2: ',K2.value[0])
print('tau2: ',tau2.value[0])
print('theta2: ',theta2.value[0])

plt.figure()
plt.subplot(2,1,1)
plt.plot(t,u1)
plt.plot(t,u2)
plt.legend([r'u1',r'u2'])
plt.ylabel('Inputs 1 & 2')
plt.subplot(2,2)
plt.plot(t,yp)
plt.legend([r'y1'])
plt.ylabel('Output')
plt.xlabel('Time')
plt.savefig('sysid.png')
plt.show()

K1:  1.0
tau1:  5.0
theta1:  10.0

K2:  0.5
tau2:  10.0
theta2:  20.0