solve_ivp求解水库充泄问题时步长异常过小的故障排查

核心函数与无降雨场景求解 #

拥有三个核心函数：

• 水位与蓄水量的关系函数 V(h)
• 坝体出水口流量与水位的关系函数 q(h)
• 降雨时长对应的入流函数 Vr(tr)

忽略降雨入流时，假设初始蓄水量V0瞬间进入水库，以下代码可成功求解放空时间：

def dtdh(h, y):
    t = y[0]
    dVdh = misc.derivative(V_h_fun, h, dx=0.1, args=(aCV, bCV, cCV))
    return (- dVdh * 1e3 / Q_h_fun(h, aCQ, bCQ, cCQ)) # V_h_fun relates h[m] to V[dam^3], so 1e3 is to get [m^3]
                                                      # as Q_h_fun gives q[m^3/s] from h[m]

elevs = np.linspace(235,231,100) # initial elev. 235, to final elev. 231 m
sol = integrate.solve_ivp(dtdh, t_span=(235,231), y0=[0], t_eval=elevs) # y0 = t0

其中V_h_fun和Q_h_fun对应上述前两个函数，aCV, bCV, cCV, aCQ, bCQ, cCQ为曲线拟合得到的系数。

加入降雨入流后的微分方程实现 #

加入降雨入流影响后，微分方程实现如下：

def dtdh2(h, y):
    t = y[0]
    dVdh = misc.derivative(V_h_fun, h, dx=0.1, args=(aCV, bCV, cCV))
    dVrdt = misc.derivative(initial_volume, t, dx=0.1, args=(V0, t_rain))
    dVrdt.clip(min=0, out=dVrdt) # the derivative goes to a large negative value before going back to 0
    return (- dVdh * 1e3 / (Q_h_fun(h, aCQ, bCQ, cCQ)  - dVrdt))

入流函数定义：

def initial_volume(t, V0, t_rain):
    # V0 in [m3]
    # t_rain in [s]
    if type(t) != np.ndarray:
        t = np.array([t])
    
    f = np.zeros_like(t)
    c1 = (t < 0)
    c2 = (t >= 0) & (t <= t_rain)
    c3 = (t > t_rain)
    f[c1] = 0
    f[c2] = V0 / t_rain * t[c2] # straight line within the duration of the rain, and 0 otherwise
    f[c3] = 0
    return f

问题现象 #

当降雨时长t_rain设置为大于9500秒时，求解器能正常收敛并给出放空时间；但当时长小于该值时，求解器在初始水位附近的步长变得极小（如234.99999996789282），无法终止，需手动中断。

问题定位与解决方案计划 #

已定位问题——降雨时蓄水量增加会导致水位上升，但水位作为自变量只能单向变化，计划分上升阶段和下降阶段分别求解。

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