Python - 求解非线性方程组（sqrt 问题）

Question

大家早上好

我正在尝试使用 python 求解非线性方程组（有损失的欠膨胀射流 - 感兴趣的物理学家的莫尔科夫）

方程组如下图所示：

为了做到这一点，我尝试使用 Scipy 库和 fsolve 模块，如下所示。

def equations(p):

    P2, T2, u2, rho2, P3, T3, u3, rho3 = p;

    return (P2 - TANK_PRES + rho2 * (u2**2) * ((K / 4) + 1),
            Cp_H2 * TANK_TEMP - Cp_H2 * T2 - ((K + 1) * (u2**2) / 2),
            P2 - ((rho2 * R_H2 * T2) / (1 - ABEL_NOBLE_COEF * rho2)),
            P3 - P2 + rho2 * (u2**2) * ((F / 4) - 1) + rho3 * (u3**2) * ((F / 4) + 1),
            Cp_H2 * T2 + ((u2**2) / 2) - Cp_H2 * T3 - ((F / 4) + 1) * ((u3**2) / 2),
            rho3 - (P3 / (ABEL_NOBLE_COEF * P3 + R_H2 * T3)),
            rho2 * u2 - rho3 * u3,
            u3 - np.sqrt((H2_GAMMA * R_H2 * T3)) / (1 - ABEL_NOBLE_COEF * rho3))


P2, T2, u2, rho2, P3, T3, u3, rho3 =  fsolve(equations, (10**7, 1000, 1500, 1, 10**7, 1000, 1500, 1))
print(equations((P2, T2, u2, rho2, P3, T3, u3, rho3)))
print("*****")
print(P2, T2, u2, rho2, P3, T3, u3, rho3)

但是，由于最后一个方程中的

np.sqrt((H2_GAMMA * R_H2 * T3))

部分，我遇到了一个问题。 Python 正在尝试解决系统尝试为 T3 取负值（因此由于平方根为负项而导致的数学错误）。 下面是收到的错误消息

(-27468750.0, -1967406.7266, -31562027.65265894, 55125.0, -13781.25, 0.7579822540902079, 0.0, -6157.533689399985)
*****
10000000.0 10000.0 1500.0 1.0 10000000.0 10000.0 1500.0 1.0
c:\users\guillaume.thiriet\downloads\test_h2jet.py:89: RuntimeWarning: invalid value encountered in sqrt
  u3 - np.sqrt((H2_GAMMA * R_H2 * T3)) / (1 - ABEL_NOBLE_COEF * rho3))
c:\users\guillaume.thiriet\downloads\test_h2jet.py:92: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.

我不知道如何解决这个问题。这个想法是能够设置边界（最小/最大值），但似乎 fsolve 不允许这样做（或者我没有找到任何相关内容）。

我在网上发现也许可以使用Scipy的最小二乘函数，但它没有给我一个合适的结果（也许我也用得不好）。

有人可以帮助我解决这个系统，或者提供一些提示来解决我的问题（其他库的其他功能也许......）？

提前感谢您的帮助

Answer 1

最简单的解决方案：

将 sqrt 中可能为负的项替换为其绝对值 x -> abs(x)，然后在优化解中取解的绝对值。

推荐解决方案：

寻找约束多变量优化求解器。

https://docs.scipy.org/doc/scipy/tutorial/optimize.html#constrained-minimization-of-multivariate-scalar-functions-minimize

或：

https://pygad.readthedocs.io/en/latest/