使用GEKKO建模水箱液位时出现无解问题求助

问题诊断 #

1. 求解失败（添加lb=0时）：
   当给液位变量h1设置下界lb=0后，模拟过程中可能出现净流量为负的时段，导致液位有下降至0以下的趋势，GEKKO的联立方程求解器无法找到满足约束的可行解，因此抛出Solution Not Found错误。

2. 结果偏差（移除lb=0后）：
   GEKKO默认采用BDF积分方法，而scipy odeint默认使用LSODA方法，两者数值积分策略存在差异；同时求解器精度设置、参数插值方式的不同，也会导致模拟结果出现偏差。

解决方案 #

1. 修复求解失败问题（保留lb=0） #

通过条件约束避免液位违反下界限制，有两种实现方式：

• 方式一：条件导数方程
  当液位大于0时正常计算导数，当液位接近0时强制导数不小于0，防止液位为负：

  q_net = m.Intermediate(qin1 + qin2 - qout)
  m.Equation(h1.dt() == m.if3(h1, q_net/Ac, m.max2(q_net/Ac, 0)))
  

• 方式二：软约束辅助
  在保留硬约束的同时，添加一个宽松的导数约束，避免求解器陷入不可行域：

  h1 = m.Var(value=d2_h0, lb=0)
  m.Equation(h1.dt() >= -0.05)  # 根据实际工况调整允许的最大下降速率
  

2. 修正结果偏差问题 #

调整GEKKO的配置，匹配scipy odeint的数值行为：

# 设置参数采用线性插值（与scipy odeint对输入数组的处理一致）
qin1 = m.Param(value=inf1, integer=False)
qin2 = m.Param(value=inf2, integer=False)
qout = m.Param(value=eff, integer=False)

# 调整模型选项提升精度
m.options.CV_TYPE = 2  # 使用二阶连续变量，提升积分平滑度
m.options.RTOL = 1e-6  # 相对精度设置
m.options.ATOL = 1e-6  # 绝对精度设置

3. 预验证液位趋势 #

先手动计算累积液位变化，确认是否存在液位为负的风险：

import numpy as np
dt = np.diff(t)
q_net = inf1 + inf2 - eff
h_calc = d2_h0 + np.cumsum(q_net[:-1]/Ac * dt)
print(f"模拟过程中最小液位预测值: {min(h_calc):.4f} m")

如果结果小于0，说明确实需要添加约束处理液位下限问题。

完整修正代码示例 #

import numpy as np
from gekko import GEKKO

# 设置时间轴
tmax = 60*6
t = np.linspace(0, tmax, int(tmax/10)+1)

# 初始化参数
d2_h0 = 13.377
inf1 = np.array([32.6354599 , 32.41882451, 32.08460871, 32.11487071, 32.71570587,
       32.59923999, 31.66669464, 30.11240896, 29.31222725, 29.35761197,
       29.62183634, 29.67505582, 29.24057325, 29.13853518, 29.48321724,
       29.61703173, 29.49874306, 28.99679947, 29.24003156, 29.40070153,
       29.70169004, 29.2913545 , 29.47371801, 29.91566467, 31.31636302,
       31.6771698 , 31.65268326, 31.06637255, 31.39147377, 31.88083331,
       32.59566625, 32.70952861, 32.78859075, 32.87391027, 32.97800064,
       32.99872208, 33.02946218])
inf2 = np.array([66.91262309, 67.16797638, 67.77143351, 66.85663605, 67.43820954,
       67.96041107, 68.7215627 , 68.91900635, 69.20062764, 68.29413096,
       68.56461334, 67.67184957, 68.84806824, 67.61451467, 69.58069102,
       71.284935  , 75.60562642, 74.83906555, 74.06419373, 71.20425161,
       69.60981496, 69.45553589, 70.35860697, 71.17754873, 72.16390737,
       72.0528005 , 72.49635569, 73.09021505, 72.7195816 , 71.9975001 ,
       70.13828532, 71.11123403, 72.16157023, 73.27675883, 71.9024353 ,
       71.17524719, 70.34394582])
eff = np.array([110.97348786, 108.6726354 , 109.4272232 , 110.57080078,
       114.20512136, 114.84948222, 113.96173604, 110.81165822,
       110.4366506 , 111.61210887, 112.75804393, 111.23046112,
       108.35852305, 108.21724955, 110.47168223, 112.10458374,
       109.28511048, 107.31727092, 108.55026245, 111.30213165,
       111.88119253, 110.62695313, 111.76373037, 115.09386699,
       115.75547282, 113.47773488, 107.95795441, 106.46175893,
       105.83562978, 109.9902064 , 110.59869131, 110.49962108,
       109.35623678, 108.35690053, 107.0867513 , 104.34462484,
       103.1198527 ])

# 创建GEKKO模型
m = GEKKO(remote=False)
m.time = t
qin1 = m.Param(value=inf1, integer=False)
qin2 = m.Param(value=inf2, integer=False)
Ac = m.Const(value=226.98)
qout = m.Param(value=eff, integer=False)
h1 = m.Var(value=d2_h0, lb=0)

# 条件导数方程避免液位为负
q_net = m.Intermediate(qin1 + qin2 - qout)
m.Equation(h1.dt() == m.if3(h1, q_net/Ac, m.max2(q_net/Ac, 0)))

# 模型配置
m.options.IMODE = 4
m.options.SOLVER = 3
m.options.CV_TYPE = 2
m.options.RTOL = 1e-6
m.options.ATOL = 1e-6

# 求解模型
m.solve(disp=True)

# 输出结果
print("模拟液位结果：", h1.value)

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