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

Question

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

有关下面求解器输出的更多详细信息。尽管我添加的变量多于使用多个输入的方程式，但这使我的问题具有负自由度。

--------- APM Model Size ------------
 Each time step contains
   Objects      :            2
   Constants    :            0
   Variables    :           15
   Intermediates:            0
   Connections  :            4
   Equations    :            6
   Residuals    :            6
 
 Number of state variables:           1918
 Number of total equations: -         2151
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :           -233
 
 * Warning: DOF <= 0
 **********************************************
 Dynamic Estimation with Interior Point Solver
 **********************************************
  
  
 Info: Exact Hessian

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.10, running with linear solver ma57.

Number of nonzeros in equality constraint Jacobian...:     6451
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     1673

Exception of type: TOO_FEW_DOF in file "IpIpoptApplication.cpp" at line 891:
 Exception message: status != TOO_FEW_DEGREES_OF_FREEDOM evaluated false: Too few degrees of freedom (rethrown)!

EXIT: Problem has too few degrees of freedom.
 
 An error occured.
 The error code is          -10
 
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   2.279999999154825E-002 sec
 Objective      :   0.000000000000000E+000
 Unsuccessful with error code            0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found

这是代码

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

Answer 1

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

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 是未测量的状态。