GEKKO无不等式模型抛出Inequality Definition错误的解决求助

错误原因 #

1. 生成器表达式解析错误：代码中用m.Equation(...) for i in range(nc)生成器表达式直接传递给m.Equation，GEKKO无法正确解析该结构，误判为无效不等式。
2. Numpy与GEKKO变量不兼容：np.sum(x * cpmol)试图用Numpy直接处理GEKKO变量数组，Numpy无法识别GEKKO变量的运算逻辑，导致形状错误。
3. 变量赋值冲突：x[3] = m.CV(value=ymeas)直接覆盖数组中原本的Var对象，后续数组运算时易引发类型不匹配问题。

修正步骤 #

• 循环添加组分方程：用for循环逐个添加4个组分的物料衡算方程，替代生成器表达式。
• 使用GEKKO内置求和函数：用m.sum()替代np.sum()，处理GEKKO变量的求和运算。
• 正确初始化CV变量：创建x数组时直接将x[3]定义为CV，避免后续赋值覆盖，确保测量数据正确关联。

完整修正代码 #

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)    # 创建GEKKO模型
m.time = t                 # 时间点

# 物性参数
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

nc = 4
# 创建GEKKO常数
Mo = m.Const(value=M[0])
Mm = m.Const(value=M[1])
VR = m.Const(20)
dHr = m.Const(-6309)

# 待辨识参数
r = m.FV(value=par0)
r.STATUS = 1  # 允许参数被调整

# 输入变量
u0 = m.Param(value=u[0, :])
u1 = m.Param(value=u[1, :])
To = m.Param(value=u[2, :])  # 直接用u2作为温度参数，简化代码
Tm = m.Param(value=u[3, :])

# 创建变量数组，x[3]直接定义为CV
x = m.Array(m.Var, nc)
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(value=333.5500)

# 反应速率项
rx = [-r, -3*r, r, 3*r]

# 中间变量
No = m.Intermediate(u0 / Mo / 3600)
Nm = m.Intermediate(u1 / Mm / 3600)
# 用GEKKO求和替代Numpy求和
cpmolR = m.Intermediate(m.sum([x[i] * cpmol[i] for i in range(nc)]))
nR = m.Intermediate(VR / m.sum([x[i] * vmol[i] for i in range(nc)]))

# 循环添加物料衡算方程
for i in range(nc):
    m.Equation(nR * x[i].dt() == Nm*(xm[i] - x[i]) + No*(xo[i] - x[i]) + rx[i]*VR)

# 能量衡算方程
m.Equation(nR * cpmolR * T.dt() == Nm*(cp[1]*M[1])*(Tm - T) + No*(cp[0]*M[0])*(To - T) + VR*(-dHr)*r)

# 求解配置
m.options.IMODE = 5  # 动态参数估计模式
m.options.NODES = 3  # 配置节点数
m.solve(disp=True)  # 显示求解过程，方便调试

print(f'辨识得到的反应速率r: {r.value[0]:.6f}')

# 可选：绘制拟合结果
plt.figure(figsize=(10,6))
plt.plot(t/3600, ymeas, 'bo', label='测量值')
plt.plot(t/3600, x[3].value, 'r-', label='拟合值')
plt.xlabel('时间(h)')
plt.ylabel('生物柴油组分')
plt.legend()
plt.show()

关键说明 #

• 用for循环逐个添加方程，避免生成器表达式导致的解析错误。
• 所有涉及GEKKO变量的求和运算都使用m.sum()，确保求解器能正确处理变量关系。
• 直接将x[3]定义为CV并关联测量数据，避免数组元素类型冲突。
• 简化了温度参数的定义，直接将u2、u3作为To、Tm的Param变量，代码更清晰。

内容的提问来源于stack exchange，提问作者Joel Sansana