scipy 非线性系统中方程的优先级排序

我试图用最小二乘法解一个非线性方程组。我得到了一个结果，但不幸的是它不是很准确。主要问题是方程29至32必须被充满，但它们不是。有没有一种方法可以在非线性系统中对方程进行优先级排序？
下面是我的代码：

# Ethanol and 1_Propanol with 3 trays

from scipy.optimize import least_squares
import numpy as np

def equations(x):
    eq1 = x[18] - x[16] * x[19] * x[11] * (x[21] / (x[19] * x[11] * x[13])) ** (2.0 / 3.0)
    eq2 = x[19] - ((732.0 * x[10]) + (747.0 * x[9]))
    eq3 = x[13] - ((1.25E-5 * x[10]) + (1.6E-6 * x[9]))
    eq4 = x[2] - ((((((-6021.54) + (0.055134232180918 * (x[29] - 298.0))) + 2.61275786690816) * x[10])
                   + ((((-6883.827) + (0.058471271259418 * (x[29] - 298.0))) + 2.38446400137128) * x[9])))
    eq5 = x[11] - (((0.055134232180918 * (x[29])) + 2.61275786690816) * x[10] + ((0.058471271259418 * (x[29])) +
                                                                                 2.38446400137128) * x[9])
    eq6 = x[21] - ((((((0.001318 * (x[29])) + 0.16762) * x[10])
                     + (((-0.001075424413 * (x[29])) + 0.155701011727667) * x[9]))))
    eq7 = x[24] - (1.746 * x[22] + 2.76 * x[20])
    eq8 = x[15] - (1.0E-5 * x[22] + 6.5E-6 * x[20])
    eq9 = x[12] - ((((((((-6021.54) + (0.015220202121093 * (x[26] - 298.0))) + 1.60395800121061) + 922.194)
                      * x[22]) + ((((((-6883.827) + (0.01847750929725 * (x[26] - 298.0))) + 1.43843773750579)
                                    + 1067.523)) * x[20]))))
    eq10 = x[6] - ((((((((0.015220202121093 * (x[26])) + 1.60395800121061)
                        * x[22]) + (((0.01847750929725 * (x[26])) + 1.43843773750579)
                                    * x[20]))))))
    eq11 = x[5] - ((((((((5.07392141746E-4 * (x[29])) + 0.015070468544221)
                        * x[22]) + (((4.89248E-4 * (x[29])) + 0.014354384)
                                    * x[20]))))))
    eq12 = x[3] - (1.0 / (10.0 ** (5.24677 - 1598.673 / (x[28] + 46.424))))
    eq13 = x[0] - (1.0 / (10.0 ** (5.31384 - 1690.864 / (x[28] + 51.804))))
    eq14 = x[27] - (x[25] * x[24] * x[6] * ((x[5] / (x[24] * x[6] * x[15])) ** (2.0 / 3.0)))
    eq15 = 0.0 - ((2.0 * 28000.0) - (x[23] * x[2]) + (x[17] * 0.002))
    eq16 = 0.0 - ((0.2 * 30000.0) - (x[14] * x[12]) - (x[17] * 0.002))
    eq17 = 0.0 - ((2.0 * 0.5) - (x[23] * x[10]) + (x[8] * 0.002))
    eq18 = 0.0 - ((2.0 * 0.5) - (x[23] * x[9]) + (x[7] * 0.002))
    eq19 = 0.0 - ((0.2 * 0.6) - (x[14] * x[22]) - (x[8] * 0.002))
    eq20 = 0.0 - ((0.2 * 0.4) - (x[14] * x[20]) - (x[7] * 0.002))
    eq21 = 0.0 - (x[31] - (x[3] * x[4]))
    eq22 = 0.0 - (x[30] - (x[0] * x[1]))
    eq23 = x[17] - (x[18] * (x[28] - x[29]) + (x[8] * (((-6021.54 + (0.055134232180918 * (x[29] - 298.0))) +
                                                        2.61275786690816) + (x[7] * (
            (-6883.827 + (0.058471271259418 * (x[29] - 298.0))) + 2.38446400137128)))))
    eq24 = x[17] - (x[27] * (x[26] - x[28]) + (x[8] * ((((-6021.54 + (0.015220202121093 * (x[26] - 298.0))) +
                                                         1.60395800121061) + 922.194) + (x[7] * (
            ((-6883.827 + (0.01847750929725 * (x[26] - 298.0))) + 1.43843773750579) + 1067.523)))))
    eq25 = x[8] - (x[19] * x[16] * (x[4] - x[10]) + x[10] * (x[8] + x[7]))
    eq26 = x[7] - (x[19] * x[16] * (x[1] - x[9]) + x[9] * (x[8] + x[7]))
    eq27 = x[8] - (x[24] * x[25] * (x[22] - x[31]) + x[22] * (x[8] + x[7]))
    eq28 = x[7] - (x[24] * x[25] * (x[20] - x[30]) + x[20] * (x[8] + x[7]))
    eq29 = 1.0 - (x[10] + x[9])
    eq30 = 1.0 - (x[4] + x[1])
    eq31 = 1.0 - (x[22] + x[20])
    eq32 = 1.0 - (x[31] + x[30])

    return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16,
            eq17, eq18, eq19, eq20, eq21, eq22, eq23, eq24, eq25, eq26, eq27, eq28, eq29, eq30,
            eq31, eq32]

# Boundaries
# Upper bounds for specific variables
upper_bounds = np.ones(32) * np.inf
upper_bounds[1] = 1
upper_bounds[4] = 1
upper_bounds[9] = 1
upper_bounds[10] = 1
upper_bounds[20] = 1
upper_bounds[22] = 1
upper_bounds[30] = 1
upper_bounds[31] = 1

# Upper bounds for specific variables
bottom_bounds = np.zeros(32)
bottom_bounds[26] = 340
bottom_bounds[28] = 340
bottom_bounds[29] = 340

# Bounds array
bounds = (bottom_bounds, upper_bounds)
# bounds = (0, np.inf)

# Initial values
x0 = np.ones(32)
x0[1] = 0.5
x0[4] = 0.5
x0[9] = 0.5
x0[10] = 0.5
x0[20] = 0.5
x0[22] = 0.5
x0[30] = 0.5
x0[31] = 0.5

x0[26] = 370
x0[28] = 370
x0[29] = 370

# x0 = np.array([1, 1e-4, 1, 1, 1, 1, 1e-4, 1, 1, 1, 1, 1, 1, 1, 1, 1])

# Lösung des Gleichungssystems mit Variablengrenzen
res = least_squares(equations, x0, bounds=bounds, method='trf', ftol=1e-10, xtol=1e-10, gtol=1e-10)

print(res.x)
print('e0_K_tr1_c2 {0}'.format(res.x[0]))
print('e0_x_Lint_tr1_c2 {0}'.format(res.x[1]))
print('e0_h_L_tr1 {0}'.format(res.x[2]))
print('e0_K_tr1_c1 {0}'.format(res.x[3]))
print('e0_x_Lint_tr1_c1 {0}'.format(res.x[4]))
print('e0_greek_lambda_V_tr1 {0}'.format(res.x[5]))
print('e0_cp_V_tr1 {0}'.format(res.x[6]))
print('e0_f_int_tr1_c2 {0}'.format(res.x[7]))
print('e0_f_int_tr1_c1 {0}'.format(res.x[8]))
print('e0_x_L_tr1_c2 {0}'.format(res.x[9]))
print('e0_x_L_tr1_c1 {0}'.format(res.x[10]))
print('e0_cp_L_tr1 {0}'.format(res.x[11]))
print('e0_h_V_tr1 {0}'.format(res.x[12]))
print('e0_D_L_tr1 {0}'.format(res.x[13]))
print('e0_F_V_tr1 {0}'.format(res.x[14]))
print('e0_D_V_tr1 {0}'.format(res.x[15]))
print('e0_greek_beta_L_tr1 {0}'.format(res.x[16]))
print('e0_e_int_tr1 {0}'.format(res.x[17]))
print('e0_greek_alpha_L_tr1 {0}'.format(res.x[18]))
print('e0_greek_rho_L_tr1 {0}'.format(res.x[19]))
print('e0_x_V_tr1_c2 {0}'.format(res.x[20]))
print('e0_greek_lambda_L_tr1 {0}'.format(res.x[21]))
print('e0_x_V_tr1_c1 {0}'.format(res.x[22]))
print('e0_F_L_tr1 {0}'.format(res.x[23]))
print('e0_greek_rho_V_tr1 {0}'.format(res.x[24]))
print('e0_greek_beta_V_tr1 {0}'.format(res.x[25]))
print('e0_T_V_tr1 {0}'.format(res.x[26]))
print('e0_greek_alpha_V_tr1 {0}'.format(res.x[27]))
print('e0_T_int_tr1 {0}'.format(res.x[28]))
print('e0_T_L_tr1 {0}'.format(res.x[29]))
print('e0_x_Vint_tr1_c2 {0}'.format(res.x[30]))
print('e0_x_Vint_tr1_c1 {0}'.format(res.x[31]))

print('Zusammensetztungen:')

print('e0_x_Vint_tr1_c1 {0}'.format(res.x[31]))
print('e0_x_Vint_tr1_c2 {0}'.format(res.x[30]))

print('e0_x_V_tr1_c1 {0}'.format(res.x[22]))
print('e0_x_V_tr1_c2 {0}'.format(res.x[20]))

print('e0_x_Lint_tr1_c1 {0}'.format(res.x[4]))
print('e0_x_Lint_tr1_c2 {0}'.format(res.x[1]))

print('e0_x_L_tr1_c1 {0}'.format(res.x[10]))
print('e0_x_L_tr1_c2 {0}'.format(res.x[9]))

字符串
以及x[31]、x[30]、x[22]、…为：

Zusammensetztungen:
e0_x_Vint_tr1_c1 0.5189962572539993
e0_x_Vint_tr1_c2 0.24082903753005624
e0_x_V_tr1_c1 2.3066615630447933e-32
e0_x_V_tr1_c2 1.246796204196686e-33
e0_x_Lint_tr1_c1 0.9999999999999999
e0_x_Lint_tr1_c2 0.01258253962013115
e0_x_L_tr1_c1 7.013309435002394e-33
e0_x_L_tr1_c2 3.0913699644636455e-33

型
e0_x_Vint_tr1_c1 + e0_x_Vint_tr1_c2应该是1，但不是。

删除多余的括号并重新设置公式的格式，以便更容易调试。
用命名变量替换索引变量。我展示了一种通过NamedTuple实现此目的的快速方法。
不幸的是，尽管方法lm似乎可以改善您的结果-它声称找到的解决方案的成本（30.7）比trf（399）低10倍-但它不可用，因为它不尊重您的边界。没有什么灵丹妙药可以提高准确性。了解领域问题是必要的。如果某些方程比其他方程对解决方案“更重要”，则可以简单地按其权重缩放这些方程。
添加回归测试。对于我展示的结果，这些方程等价于你原来的方程。
始终检查result.success。
第一遍看起来像

"""Ethanol and 1_Propanol with 3 trays"""

from pprint import pprint
from typing import NamedTuple

from numpy.random import default_rng
from scipy.optimize import least_squares
import numpy as np

class Vars(NamedTuple):
    e0_K_tr1_c2: float            # 0
    e0_x_Lint_tr1_c2: float       # 1
    e0_h_L_tr1: float             # 2
    e0_K_tr1_c1: float            # 3
    e0_x_Lint_tr1_c1: float       # 4
    e0_lambda_V_tr1: float        # 5
    e0_cp_V_tr1: float            # 6
    e0_f_int_tr1_c2: float        # 7
    e0_f_int_tr1_c1: float        # 8
    e0_x_L_tr1_c2: float          # 9
    e0_x_L_tr1_c1: float          # 10
    e0_cp_L_tr1: float            # 11
    e0_h_V_tr1: float             # 12
    e0_D_L_tr1: float             # 13
    e0_F_V_tr1: float             # 14
    e0_D_V_tr1: float             # 15
    e0_beta_L_tr1: float          # 16
    e0_e_int_tr1: float           # 17
    e0_alpha_L_tr1: float         # 18
    e0_rho_L_tr1: float           # 19
    e0_x_V_tr1_c2: float          # 20
    e0_lambda_L_tr1: float        # 21
    e0_x_V_tr1_c1: float          # 22
    e0_F_L_tr1: float             # 23
    e0_rho_V_tr1: float           # 24
    e0_beta_V_tr1: float          # 25
    e0_T_V_tr1: float             # 26
    e0_alpha_V_tr1: float         # 27
    e0_T_int_tr1: float           # 28
    e0_T_L_tr1: float             # 29
    e0_x_Vint_tr1_c2: float       # 30
    e0_x_Vint_tr1_c1: float       # 31

def equations(x: np.ndarray) -> tuple[float, ...]:
    v = Vars(*x)

    eq1 = v.e0_alpha_L_tr1 - v.e0_beta_L_tr1 * v.e0_rho_L_tr1 * v.e0_cp_L_tr1 * (
        v.e0_lambda_L_tr1 / (v.e0_rho_L_tr1 * v.e0_cp_L_tr1 * v.e0_D_L_tr1)
    ) ** (2 / 3)

    eq2 = v.e0_rho_L_tr1 - (732 * v.e0_x_L_tr1_c1 + 747 * v.e0_x_L_tr1_c2)
    eq3 = v.e0_D_L_tr1 - (1.25E-5 * v.e0_x_L_tr1_c1 + 1.6E-6 * v.e0_x_L_tr1_c2)

    eq4 = v.e0_h_L_tr1 - (
        (
            -6021.54
            + 0.055134232180918 * (v.e0_T_L_tr1 - 298)
            + 2.61275786690816
        ) * v.e0_x_L_tr1_c1
        + (
            -6883.827
            + 0.058471271259418 * (v.e0_T_L_tr1 - 298)
            + 2.38446400137128
        ) * v.e0_x_L_tr1_c2
    )

    eq5 = v.e0_cp_L_tr1 - (
        (
            0.055134232180918 * v.e0_T_L_tr1 + 2.61275786690816
        ) * v.e0_x_L_tr1_c1
        + (
            0.058471271259418 * v.e0_T_L_tr1 + 2.38446400137128
        ) * v.e0_x_L_tr1_c2
    )
    eq6 = v.e0_lambda_L_tr1 - (
        (
            0.001318 * v.e0_T_L_tr1 + 0.16762
        ) * v.e0_x_L_tr1_c1
        + (
            -0.001075424413 * v.e0_T_L_tr1 + 0.155701011727667
        ) * v.e0_x_L_tr1_c2
    )

    eq7 = v.e0_rho_V_tr1 - (1.746 * v.e0_x_V_tr1_c1 + 2.76 * v.e0_x_V_tr1_c2)
    eq8 = v.e0_D_V_tr1 - (1E-5 * v.e0_x_V_tr1_c1 + 6.5E-6 * v.e0_x_V_tr1_c2)

    eq9 = v.e0_h_V_tr1 - (
        (
            -6021.54
            + 0.015220202121093 * (v.e0_T_V_tr1 - 298)
            + 1.60395800121061
            + 922.194
        ) * v.e0_x_V_tr1_c1
        + (
            -6883.827
            + 0.01847750929725 * (v.e0_T_V_tr1 - 298)
            + 1.43843773750579
            + 1067.523
        ) * v.e0_x_V_tr1_c2
    )
    eq10 = v.e0_cp_V_tr1 - (
        (
            0.015220202121093 * v.e0_T_V_tr1 + 1.60395800121061
        ) * v.e0_x_V_tr1_c1
        + (
            0.01847750929725 * v.e0_T_V_tr1 + 1.43843773750579
        ) * v.e0_x_V_tr1_c2
    )
    eq11 = v.e0_lambda_V_tr1 - (
        (
            5.07392141746E-4 * v.e0_T_L_tr1 + 0.015070468544221
        ) * v.e0_x_V_tr1_c1
        + (
            4.89248E-4 * v.e0_T_L_tr1 + 0.014354384
        ) * v.e0_x_V_tr1_c2
    )

    eq12 = v.e0_K_tr1_c1 - 1 / 10**(5.24677 - 1598.673 / (v.e0_T_int_tr1 + 46.424))
    eq13 = v.e0_K_tr1_c2 - 1 / 10**(5.31384 - 1690.864 / (v.e0_T_int_tr1 + 51.804))

    eq14 = v.e0_alpha_V_tr1 - (
        v.e0_beta_V_tr1 * v.e0_rho_V_tr1 * v.e0_cp_V_tr1 * (
            v.e0_lambda_V_tr1 / (
                v.e0_rho_V_tr1 * v.e0_cp_V_tr1 * v.e0_D_V_tr1
            )
        ) ** (2 / 3)
    )

    eq15 = -(2 * 28000 - v.e0_F_L_tr1 * v.e0_h_L_tr1 + v.e0_e_int_tr1 * 0.002)
    eq16 = -(0.2 * 30000 - v.e0_F_V_tr1 * v.e0_h_V_tr1 - v.e0_e_int_tr1 * 0.002)
    eq17 = -(2 * 0.5 - v.e0_F_L_tr1 * v.e0_x_L_tr1_c1 + v.e0_f_int_tr1_c1 * 0.002)
    eq18 = -(2 * 0.5 - v.e0_F_L_tr1 * v.e0_x_L_tr1_c2 + v.e0_f_int_tr1_c2 * 0.002)
    eq19 = -(0.2 * 0.6 - v.e0_F_V_tr1 * v.e0_x_V_tr1_c1 - v.e0_f_int_tr1_c1 * 0.002)
    eq20 = -(0.2 * 0.4 - v.e0_F_V_tr1 * v.e0_x_V_tr1_c2 - v.e0_f_int_tr1_c2 * 0.002)

    eq21 = -(v.e0_x_Vint_tr1_c1 - v.e0_K_tr1_c1 * v.e0_x_Lint_tr1_c1)
    eq22 = -(v.e0_x_Vint_tr1_c2 - v.e0_K_tr1_c2 * v.e0_x_Lint_tr1_c2)

    eq23 = v.e0_e_int_tr1 - (
        v.e0_alpha_L_tr1 * (v.e0_T_int_tr1 - v.e0_T_L_tr1)
        + v.e0_f_int_tr1_c1 * (
            -6021.54
            + 0.055134232180918 * (v.e0_T_L_tr1 - 298)
            + 2.61275786690816
            # really?
        # )
            + v.e0_f_int_tr1_c2 * (
                -6883.827 +
                0.058471271259418 * (v.e0_T_L_tr1 - 298)
                + 2.38446400137128
            )
        )
    )

    eq24 = v.e0_e_int_tr1 - (
        v.e0_alpha_V_tr1 * (v.e0_T_V_tr1 - v.e0_T_int_tr1)
        + v.e0_f_int_tr1_c1 * (
            -6021.54
            + 0.015220202121093 * (v.e0_T_V_tr1 - 298)
            + 1.60395800121061
            + 922.194
            # really?
        # )
            + v.e0_f_int_tr1_c2 * (
                -6883.827
                + 0.01847750929725 * (v.e0_T_V_tr1 - 298)
                + 1.43843773750579
                + 1067.523
            )
        )
    )

    eq25 = v.e0_f_int_tr1_c1 - (
        v.e0_rho_L_tr1 * v.e0_beta_L_tr1 * (
            v.e0_x_Lint_tr1_c1 - v.e0_x_L_tr1_c1
        )
        + v.e0_x_L_tr1_c1 * (v.e0_f_int_tr1_c1 + v.e0_f_int_tr1_c2)
    )
    eq26 = v.e0_f_int_tr1_c2 - (
        v.e0_rho_L_tr1 * v.e0_beta_L_tr1 * (
            v.e0_x_Lint_tr1_c2 - v.e0_x_L_tr1_c2
        )
        + v.e0_x_L_tr1_c2 * (v.e0_f_int_tr1_c1 + v.e0_f_int_tr1_c2)
    )
    eq27 = v.e0_f_int_tr1_c1 - (
        v.e0_rho_V_tr1 * v.e0_beta_V_tr1 * (
            v.e0_x_V_tr1_c1 - v.e0_x_Vint_tr1_c1
        )
        + v.e0_x_V_tr1_c1 * (v.e0_f_int_tr1_c1 + v.e0_f_int_tr1_c2)
    )
    eq28 = v.e0_f_int_tr1_c2 - (
        v.e0_rho_V_tr1 * v.e0_beta_V_tr1 * (
            v.e0_x_V_tr1_c2 - v.e0_x_Vint_tr1_c2
        )
        + v.e0_x_V_tr1_c2 * (v.e0_f_int_tr1_c1 + v.e0_f_int_tr1_c2)
    )

    eq29 = 1 - (v.e0_x_L_tr1_c1 + v.e0_x_L_tr1_c2)
    eq30 = 1 - (v.e0_x_Lint_tr1_c1 + v.e0_x_Lint_tr1_c2)
    eq31 = 1 - (v.e0_x_V_tr1_c1 + v.e0_x_V_tr1_c2)
    eq32 = 1 - (v.e0_x_Vint_tr1_c1 + v.e0_x_Vint_tr1_c2)

    return (eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10,
            eq11, eq12, eq13, eq14, eq15, eq16, eq17, eq18, eq19, eq20,
            eq21, eq22, eq23, eq24, eq25, eq26, eq27, eq28, eq29, eq30,
            eq31, eq32)

def regression_test(x0: np.ndarray) -> None:
    y = equations(x0)
    assert np.allclose(
        y,
        [ 0.00000000e+00, -7.38500000e+02,  9.99992950e-01,  6.44709509e+03,
         -2.25156291e+01,  7.93463011e-01, -1.25300000e+00,  9.99991750e-01,
          5.45609068e+03, -6.75527448e+00,  8.00909148e-01,  9.60890608e-01,
          9.50476815e-01,  0.00000000e+00, -5.59990020e+04, -5.99899800e+03,
         -5.02000000e-01, -5.02000000e-01,  3.82000000e-01,  4.22000000e-01,
          0.00000000e+00,  0.00000000e+00,  1.28931902e+04,  1.09111814e+04,
          0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
          0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
          ],
        rtol=0, atol=1e-4,
    )

    rand = default_rng(seed=0)
    x = rand.uniform(low=0, high=2, size=32)
    y = equations(x)
    assert np.allclose(
        y,
        [-8.89664807e-02, -2.59056243e+03,  6.71477620e-02,  2.27496385e+04,
         -8.97636119e+00, -3.16304973e-01, -1.26737482e+00,  3.51297461e-01,
          7.17475798e+03, -1.06230892e+00,  1.80356351e+00, -1.59115277e+28,
         -3.04937193e+26, -6.69625188e-01, -5.59998961e+04, -5.99749540e+03,
          1.10987295e+00,  1.41775827e+00,  1.83947353e+00,  5.57239541e-03,
         -7.24077112e-01, -6.89518259e-01,  1.75063897e+04,  1.47800191e+04,
         -3.05993271e+00, -1.36098936e+00, -2.85999752e+00,  2.56167351e+00,
         -2.50185196e+00, -1.16611391e+00, -3.97888172e-01, -1.15473631e+00,
         ],
        rtol=1e-8, atol=0,
    )

def make_x0() -> np.ndarray:
    x0 = np.ones(32)
    x0[[1,4,9,10,20,22,30,31]] = 0.5
    x0[[26,28,29]] = 370
    return x0

def solve(x0: np.ndarray) -> Vars:
    upper_bounds = np.full(shape=32, fill_value=np.inf)
    upper_bounds[[1,4,9,10,20,22,30,31]] = 1

    bottom_bounds = np.zeros(32)
    bottom_bounds[[26,28,29]] = 340

    res = least_squares(
        fun=equations, x0=x0,
        bounds=(bottom_bounds, upper_bounds),  # not supported in lm
        method='trf',
        ftol=1e-10, xtol=1e-10, gtol=1e-10,
    )
    assert res.success
    return Vars(*res.x)

def dump(x: Vars) -> None:
    pprint(x._asdict())

def main() -> None:
    x0 = make_x0()
    regression_test(x0)
    x = solve(x0)
    dump(x)

if __name__ == '__main__':
    main()

个字符

scipy 非线性系统中方程的优先级排序

1条答案

相关问题

热门标签

最新问答