我正在try 使用Gekko来解决一个动态参数辨识问题. 作为背景，这是一个生产生物柴油的CSTR，我试图在其中找到平均react 速率.
x[3]是react 堆输出中的生物柴油分数，我正在try 将其与来自实验室测量的ymeas相匹配.u是具有四个负载变量的测量值的矩阵.
给出了以下公式:
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
t = np.arange(0, 36000, 3600)
ymeas = np.array([0.429, 0.429, 0.427, 0.429, 0.428, 0.429, 0.429, 0.428, 0.428, 0.429])
u = np.array([
[2996.9, 2989.4, 2997.0, 3005.6, 3001.7, 2988.8, 2986.6, 3005.5, 2993.9, 2994.0,],
[656.9, 656.8, 657.7, 658.8, 657.3, 660.2, 659.5, 662.0, 659.7, 659.0,],
[333.1, 334.1, 332.6, 333.6, 331.0, 330.9, 331.0, 330.8, 331.7, 331.3,],
[323.2, 324.2, 324.9, 325.9, 326.9, 327.11, 327.6, 327.6, 327.4, 326.8],
])
par0=0.046
m = GEKKO(remote=False) # create GEKKO model
m.time = t # time points
M = np.array([0.853, 0.032, 0.092, 0.286,])
ro = np.array([954.0, 757.0, 1340.0, 844.0,])
cp = np.array([2110.0, 2785.0, 2556.0, 2146.0,])
xo = np.array([1, 0, 0, 0,])
xm = np.array([0, 1, 0, 0,])
vmol = M / ro
cpmol = cp * M
# create GEKKO constants
nc = 4
Mo = m.Const(value=0.853)
Mm = m.Const(value=0.032)
cpmolo = m.Const(value=cpmol[0])
cpmolm = m.Const(value=cpmol[1])
VR = m.Const(20)
dHr = m.Const(-6309)
# create GEKKO parameter
r = m.FV(value=par0)
u0 = m.Param(value=u[0, :])
u1 = m.Param(value=u[1, :])
u2 = m.Param(value=u[2, :])
u3 = m.Param(value=u[3, :])
# create GEKKO variables
x = m.Array(m.Var, 4)
x[0].value = 0.0031
x[1].value = 0.4235
x[2].value = 0.1432
x[3] = m.CV(value=ymeas)
x[3].FSTATUS = 1
T = m.Var(333.5500)
rx = [-r, -3*r, r, 3*r,]
# create GEKKO equations
To = u2
Tm = u3
No = m.Intermediate(u0/Mo/3600)
Nm = m.Intermediate(u1/Mm/3600)
cpmolR = m.Intermediate(np.sum(x * cpmol))
nR = m.Intermediate(VR / np.sum(vmol * x))
m.Equation(nR*x[i].dt() == Nm*(xm[i] - x[i]) + No*(xo[i] - x[i]) + rx[i]*VR for i in range(nc))
m.Equation(nR*cpmolR*T.dt() == Nm*cpmolm*(Tm - T) + No*cpmolo*(To - T) + VR*(-dHr)*r)
# solve ODE
r.STATUS = 1
m.options.IMODE = 5 # dynamic estimation
m.options.NODES = 3 # collocation nodes
m.solve(disp=False)
print('Parameter identification: done!')
我得到以下错误:
Exception: @error: Inequality Definition
invalid inequalities: z > x < y
<generatorobject<genexpr>at0x7f88ecb7cf90>
STOPPING . . .
我想不出解决这个问题的办法，因为我没有定义不平等.
有谁有办法解决这个问题，让解算器找到r？
我try 将x[3]创建为独立的y，并将其附加到x数组中，但在中间公式cpmolR = m.Intermediate(np.sum(x * cpmol))中遇到形状错误.