优先考虑非线性系统中的方程。

huangapple
go评论127阅读模式
英文:

Prioritize equations in a nonlinear system

问题

我正在尝试使用最小二乘法解决非线性等式系统（EQS）的问题。我得到了一个结果，但不幸的是它并不非常准确。主要问题是方程29到32必须被满足，但它们并没有被满足。有没有一种方法可以在非线性系统中优先考虑方程？

以下是我的代码：

  1. # 乙醇和1-丙醇与3个塔板
  3. from scipy.optimize import least_squares
  4. import numpy as np
  7. def equations(x):
  8.     # 你的方程式在这里
  10. # 界限
  11. # 特定变量的上界
  12. upper_bounds = np.ones(32) * np.inf
  13. # 你的上界在这里
  15. # 特定变量的下界
  16. bottom_bounds = np.zeros(32)
  17. # 你的下界在这里
  19. # 边界数组
  20. bounds = (bottom_bounds, upper_bounds)
  22. # 初始值
  23. x0 = np.ones(32)
  24. # 你的初始值在这里
  26. # 使用变量边界解决方程组
  27. res = least_squares(equations, x0, bounds=bounds, method='trf', ftol=1e-10, xtol=1e-10, gtol=1e-10)
  29. # 打印结果
  30. print(res.x)
  31. # 你的结果打印在这里

e0_x_Vint_tr1_c1 和 e0_x_Vint_tr1_c2 的结果似乎不正确，因为它们的总和不等于1。您可能需要检查方程29和32以确保它们正确地表示了变量之间的关系。如果这些方程正确，但结果仍然不准确，您可能需要进一步调整初始值或尝试不同的优化方法。

英文:

I am trying to solve a nonlinear EQS using the least squares method. I get a result, but unfortunately it is not very accurate. The main problem is that the equations 29 to 32 must be fullfilled, but they are not. Is there a way to maybe prioritize equations in a nonlinear system?

Here is the my code:

  1. # Ethanol and 1_Propanol with 3 trays
  3. from scipy.optimize import least_squares
  4. import numpy as np
  7. def equations(x):
  8.     eq1 = x[18] - x[16] * x[19] * x[11] * (x[21] / (x[19] * x[11] * x[13])) ** (2.0 / 3.0)
  9.     eq2 = x[19] - ((732.0 * x[10]) + (747.0 * x[9]))
  10.     eq3 = x[13] - ((1.25E-5 * x[10]) + (1.6E-6 * x[9]))
  11.     eq4 = x[2] - ((((((-6021.54) + (0.055134232180918 * (x[29] - 298.0))) + 2.61275786690816) * x[10])
  12.                    + ((((-6883.827) + (0.058471271259418 * (x[29] - 298.0))) + 2.38446400137128) * x[9])))
  13.     eq5 = x[11] - (((0.055134232180918 * (x[29])) + 2.61275786690816) * x[10] + ((0.058471271259418 * (x[29])) +
  14.                                                                                  2.38446400137128) * x[9])
  15.     eq6 = x[21] - ((((((0.001318 * (x[29])) + 0.16762) * x[10])
  16.                      + (((-0.001075424413 * (x[29])) + 0.155701011727667) * x[9]))))
  17.     eq7 = x[24] - (1.746 * x[22] + 2.76 * x[20])
  18.     eq8 = x[15] - (1.0E-5 * x[22] + 6.5E-6 * x[20])
  19.     eq9 = x[12] - ((((((((-6021.54) + (0.015220202121093 * (x[26] - 298.0))) + 1.60395800121061) + 922.194)
  20.                       * x[22]) + ((((((-6883.827) + (0.01847750929725 * (x[26] - 298.0))) + 1.43843773750579)
  21.                                     + 1067.523)) * x[20]))))
  22.     eq10 = x[6] - ((((((((0.015220202121093 * (x[26])) + 1.60395800121061)
  23.                         * x[22]) + (((0.01847750929725 * (x[26])) + 1.43843773750579)
  24.                                     * x[20]))))))
  25.     eq11 = x[5] - ((((((((5.07392141746E-4 * (x[29])) + 0.015070468544221)
  26.                         * x[22]) + (((4.89248E-4 * (x[29])) + 0.014354384)
  27.                                     * x[20]))))))
  28.     eq12 = x[3] - (1.0 / (10.0 ** (5.24677 - 1598.673 / (x[28] + 46.424))))
  29.     eq13 = x[0] - (1.0 / (10.0 ** (5.31384 - 1690.864 / (x[28] + 51.804))))
  30.     eq14 = x[27] - (x[25] * x[24] * x[6] * ((x[5] / (x[24] * x[6] * x[15])) ** (2.0 / 3.0)))
  31.     eq15 = 0.0 - ((2.0 * 28000.0) - (x[23] * x[2]) + (x[17] * 0.002))
  32.     eq16 = 0.0 - ((0.2 * 30000.0) - (x[14] * x[12]) - (x[17] * 0.002))
  33.     eq17 = 0.0 - ((2.0 * 0.5) - (x[23] * x[10]) + (x[8] * 0.002))
  34.     eq18 = 0.0 - ((2.0 * 0.5) - (x[23] * x[9]) + (x[7] * 0.002))
  35.     eq19 = 0.0 - ((0.2 * 0.6) - (x[14] * x[22]) - (x[8] * 0.002))
  36.     eq20 = 0.0 - ((0.2 * 0.4) - (x[14] * x[20]) - (x[7] * 0.002))
  37.     eq21 = 0.0 - (x[31] - (x[3] * x[4]))
  38.     eq22 = 0.0 - (x[30] - (x[0] * x[1]))
  39.     eq23 = x[17] - (x[18] * (x[28] - x[29]) + (x[8] * (((-6021.54 + (0.055134232180918 * (x[29] - 298.0))) +
  40.                                                         2.61275786690816) + (x[7] * (
  41.             (-6883.827 + (0.058471271259418 * (x[29] - 298.0))) + 2.38446400137128)))))
  42.     eq24 = x[17] - (x[27] * (x[26] - x[28]) + (x[8] * ((((-6021.54 + (0.015220202121093 * (x[26] - 298.0))) +
  43.                                                          1.60395800121061) + 922.194) + (x[7] * (
  44.             ((-6883.827 + (0.01847750929725 * (x[26] - 298.0))) + 1.43843773750579) + 1067.523)))))
  45.     eq25 = x[8] - (x[19] * x[16] * (x[4] - x[10]) + x[10] * (x[8] + x[7]))
  46.     eq26 = x[7] - (x[19] * x[16] * (x[1] - x[9]) + x[9] * (x[8] + x[7]))
  47.     eq27 = x[8] - (x[24] * x[25] * (x[22] - x[31]) + x[22] * (x[8] + x[7]))
  48.     eq28 = x[7] - (x[24] * x[25] * (x[20] - x[30]) + x[20] * (x[8] + x[7]))
  49.     eq29 = 1.0 - (x[10] + x[9])
  50.     eq30 = 1.0 - (x[4] + x[1])
  51.     eq31 = 1.0 - (x[22] + x[20])
  52.     eq32 = 1.0 - (x[31] + x[30])
  54.     return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16,
  55.             eq17, eq18, eq19, eq20, eq21, eq22, eq23, eq24, eq25, eq26, eq27, eq28, eq29, eq30,
  56.             eq31, eq32]
  59. # Boundaries
  60. # Upper bounds for specific variables
  61. upper_bounds = np.ones(32) * np.inf
  62. upper_bounds[1] = 1
  63. upper_bounds[4] = 1
  64. upper_bounds[9] = 1
  65. upper_bounds[10] = 1
  66. upper_bounds[20] = 1
  67. upper_bounds[22] = 1
  68. upper_bounds[30] = 1
  69. upper_bounds[31] = 1
  71. # Upper bounds for specific variables
  72. bottom_bounds = np.zeros(32)
  73. bottom_bounds[26] = 340
  74. bottom_bounds[28] = 340
  75. bottom_bounds[29] = 340
  77. # Bounds array
  78. bounds = (bottom_bounds, upper_bounds)
  79. # bounds = (0, np.inf)
  82. # Initial values
  83. x0 = np.ones(32)
  84. x0[1] = 0.5
  85. x0[4] = 0.5
  86. x0[9] = 0.5
  87. x0[10] = 0.5
  88. x0[20] = 0.5
  89. x0[22] = 0.5
  90. x0[30] = 0.5
  91. x0[31] = 0.5
  93. x0[26] = 370
  94. x0[28] = 370
  95. x0[29] = 370
  98. # x0 = np.array([1, 1e-4, 1, 1, 1, 1, 1e-4, 1, 1, 1, 1, 1, 1, 1, 1, 1])
  101. # Lösung des Gleichungssystems mit Variablengrenzen
  102. res = least_squares(equations, x0, bounds=bounds, method='trf', ftol=1e-10, xtol=1e-10, gtol=1e-10)
  104. print(res.x)
  105. print('e0_K_tr1_c2 {0}'.format(res.x[0]))
  106. print('e0_x_Lint_tr1_c2 {0}'.format(res.x[1]))
  107. print('e0_h_L_tr1 {0}'.format(res.x[2]))
  108. print('e0_K_tr1_c1 {0}'.format(res.x[3]))
  109. print('e0_x_Lint_tr1_c1 {0}'.format(res.x[4]))
  110. print('e0_greek_lambda_V_tr1 {0}'.format(res.x[5]))
  111. print('e0_cp_V_tr1 {0}'.format(res.x[6]))
  112. print('e0_f_int_tr1_c2 {0}'.format(res.x[7]))
  113. print('e0_f_int_tr1_c1 {0}'.format(res.x[8]))
  114. print('e0_x_L_tr1_c2 {0}'.format(res.x[9]))
  115. print('e0_x_L_tr1_c1 {0}'.format(res.x[10]))
  116. print('e0_cp_L_tr1 {0}'.format(res.x[11]))
  117. print('e0_h_V_tr1 {0}'.format(res.x[12]))
  118. print('e0_D_L_tr1 {0}'.format(res.x[13]))
  119. print('e0_F_V_tr1 {0}'.format(res.x[14]))
  120. print('e0_D_V_tr1 {0}'.format(res.x[15]))
  121. print('e0_greek_beta_L_tr1 {0}'.format(res.x[16]))
  122. print('e0_e_int_tr1 {0}'.format(res.x[17]))
  123. print('e0_greek_alpha_L_tr1 {0}'.format(res.x[18]))
  124. print('e0_greek_rho_L_tr1 {0}'.format(res.x[19]))
  125. print('e0_x_V_tr1_c2 {0}'.format(res.x[20]))
  126. print('e0_greek_lambda_L_tr1 {0}'.format(res.x[21]))
  127. print('e0_x_V_tr1_c1 {0}'.format(res.x[22]))
  128. print('e0_F_L_tr1 {0}'.format(res.x[23]))
  129. print('e0_greek_rho_V_tr1 {0}'.format(res.x[24]))
  130. print('e0_greek_beta_V_tr1 {0}'.format(res.x[25]))
  131. print('e0_T_V_tr1 {0}'.format(res.x[26]))
  132. print('e0_greek_alpha_V_tr1 {0}'.format(res.x[27]))
  133. print('e0_T_int_tr1 {0}'.format(res.x[28]))
  134. print('e0_T_L_tr1 {0}'.format(res.x[29]))
  135. print('e0_x_Vint_tr1_c2 {0}'.format(res.x[30]))
  136. print('e0_x_Vint_tr1_c1 {0}'.format(res.x[31]))
  138. print('Zusammensetztungen:')
  140. print('e0_x_Vint_tr1_c1 {0}'.format(res.x[31]))
  141. print('e0_x_Vint_tr1_c2 {0}'.format(res.x[30]))
  143. print('e0_x_V_tr1_c1 {0}'.format(res.x[22]))
  144. print('e0_x_V_tr1_c2 {0}'.format(res.x[20]))
  146. print('e0_x_Lint_tr1_c1 {0}'.format(res.x[4]))
  147. print('e0_x_Lint_tr1_c2 {0}'.format(res.x[1]))
  149. print('e0_x_L_tr1_c1 {0}'.format(res.x[10]))
  150. print('e0_x_L_tr1_c2 {0}'.format(res.x[9]))

and the results for x[31], x[30], x[22],... are:

  1. Zusammensetztungen:
  2. e0_x_Vint_tr1_c1 0.5189962572539993
  3. e0_x_Vint_tr1_c2 0.24082903753005624
  4. e0_x_V_tr1_c1 2.3066615630447933e-32
  5. e0_x_V_tr1_c2 1.246796204196686e-33
  6. e0_x_Lint_tr1_c1 0.9999999999999999
  7. e0_x_Lint_tr1_c2 0.01258253962013115
  8. e0_x_L_tr1_c1 7.013309435002394e-33
  9. e0_x_L_tr1_c2 3.0913699644636455e-33

but e0_x_Vint_tr1_c1 + e0_x_Vint_tr1_c2 should be 1, but it's not.

答案1

得分: 1

去掉多余的括号并重新格式化你的方程，以便更容易调试。

用命名变量替换你的索引变量。我通过NamedTuple展示了一种快速的方法。

不幸的是，尽管方法lm似乎可以改善你的结果 - 它声称找到的解决方案的成本比trf低约10倍（30.7比399） - 但它不能使用，因为它不遵守你的边界。没有一个神奇的方法可以提高准确性。了解领域问题是必要的。如果某些方程“比其他方程更重要”，你可以简单地按照它们的权重来缩放这些方程。

添加回归测试。我展示的结果中，方程与你的原始方程等效。

始终检查result.success。

第一次尝试可以如下所示：

  1. """乙醇和1-丙醇，3个托盘"""
  3. from pprint import pprint
  4. from typing import NamedTuple
  6. from numpy.random import default_rng
  7. from scipy.optimize import least_squares
  8. import numpy as np
  11. class Vars(NamedTuple):
  12.     e0_K_tr1_c2: float            # 0
  13.     e0_x_Lint_tr1_c2: float       # 1
  14.     e0_h_L_tr1: float             # 2
  15.     e0_K_tr1_c1: float            # 3
  16.     e0_x_Lint_tr1_c1: float       # 4
  17.     e0_lambda_V_tr1: float        # 5
  18.     e0_cp_V_tr1: float            # 6
  19.     e0_f_int_tr1_c2: float        # 7
  20.     e0_f_int_tr1_c1: float        # 8
  21.     e0_x_L_tr1_c2: float          # 9
  22.     e0_x_L_tr1_c1: float          # 10
  23.     e0_cp_L_tr1: float            # 11
  24.     e0_h_V_tr1: float             # 12
  25.     e0_D_L_tr1: float             # 13
  26.     e0_F_V_tr1: float             # 14
  27.     e0_D_V_tr1: float             # 15
  28.     e0_beta_L_tr1: float          # 16
  29.     e0_e_int_tr1: float           # 17
  30.     e0_alpha_L_tr1: float         # 18
  31.     e0_rho_L_tr1: float           # 19
  32.     e0_x_V_tr1_c2: float          # 20
  33.     e0_lambda_L_tr1: float        # 21
  34.     e0_x_V_tr1_c1: float          # 22
  35.     e0_F_L_tr1: float             # 23
  36.     e0_rho_V_tr1: float           # 24
  37.     e0_beta_V_tr1: float          # 25
  38.     e0_T_V_tr1: float             # 26
  39.     e0_alpha_V_tr1: float         # 27
  40.     e0_T_int_tr1: float           # 28
  41.     e0_T_L_tr1: float             # 29
  42.     e0_x_Vint_tr1_c2: float       # 30
  43.     e0_x_Vint_tr1_c1: float       # 31
  46. def equations(x: np.ndarray) -> tuple[float, ...]:
  47.     v = Vars(*x)
  49.     eq1 = v.e0_alpha_L_tr1 - v.e0_beta_L_tr1 * v.e0_rho_L_tr1 * v.e0_cp_L_tr1 * (
  50.         v.e0_lambda_L_tr1 / (v.e0_rho_L_tr1 * v.e0_cp_L_tr1 * v.e0_D_L_tr1)
  51.     ) ** (2 / 3)
  53.     eq2 = v.e0_rho_L_tr1 - (732 * v.e0_x_L_tr1_c1 + 747 * v.e0_x_L_tr1_c2)
  54.     eq3 = v.e0_D_L_tr1 - (1.25E-5 * v.e0_x_L_tr1_c1 + 1.6E-6 * v.e0_x_L_tr1_c2)
  56.     eq4 = v.e0_h_L_tr1 - (
  57.         (
  58.             -6021.54
  59.             + 0.055134232180918 * (v.e0_T_L_tr1 - 298)
  60.             + 2.61275786690816
  61.         ) * v.e0_x_L_tr1_c1
  62.         + (
  63.             -6883.827
  64.             + 0.058471271259418 * (v.e0_T_L_tr1 - 298)
  65.             + 2.38446400137128
  66.         ) * v.e0_x_L_tr1_c2
  67.     )
  69.     eq5 = v.e0_cp_L_tr1 - (
  70.         (
  71.             0.055134232180918 * v.e0_T_L_tr1 + 2.61275786690816
  72.         ) * v.e0_x_L_tr1_c1
  73.         + (
  74.             0.058471271259418 * v.e0_T_L_tr1 + 2.38446400137128
  75.         ) * v.e0_x_L_tr1_c2
  76.     )
  77.     eq6 = v.e0_lambda_L_tr1 - (
  78.         (
  79.             0.001318 * v.e0_T_L_tr1 + 0.16762
  80.         ) * v.e0_x_L_tr1_c1
  81.         + (
  82.             -0.001075424413 * v.e0_T_L_tr1 + 0.155701011727667
  83.         ) * v.e0_x_L_tr1_c2
  84.     )
  86.     eq7 = v.e0_rho_V_tr1 - (1.746 * v.e0_x_V_tr1_c1 + 2.76 * v.e0_x_V_tr1_c2)
  87.     eq8 = v.e0_D_V_tr1 - (1E-5 * v.e0_x_V_tr1_c1 + 6.5E-6 * v.e0_x_V_tr1_c2)
  89.     eq9 = v.e0_h_V_tr1 - (
  90.         (
  91.             -6021.54
  92.             + 0.015220202121093 * (v.e0_T_V_tr1 - 298)
  93.             + 1.60395800121061
  94.             + 922.194
  95.         ) * v.e0_x_V_tr1_c1
  96.         + (
  97.             -6883.827
  98.             + 0.01847750929725 * (v.e0_T_V_tr1 - 298)
  99.             + 1.43843773750579
  100.             + 1067.523
  101.         ) * v.e0_x_V_tr1_c2
  102.     )
  103.     eq10 = v.e0_cp
  105. <details>
  106. <summary>英文:</summary>
  108. Remove your redundant parentheses and reformat your equations so that they can be more easily debugged.
  110. Replace your indexed variables with named variables. I show a quick way to do this via `NamedTuple`.
  112. Unfortunately, whereas method `lm` may seem to improve your results - it claims to find a solution with ~ 10 times less cost (30.7) than that from `trf` (399) - it&#39;s not usable because it doesn&#39;t respect your bounds. There is no magic bullet to improve the accuracy. Knowledge of the domain problem is necessary. If some equations &quot;matter more&quot; to the solution than others, you can simply scale those equations by their weights.
  114. Add regression tests. For the results that I show, the equations are equivalent to your original equations.
  116. Always check `result.success`.
  118. A first pass can look like
  120. ```python
  121. &quot;&quot;&quot;Ethanol and 1_Propanol with 3 trays&quot;&quot;&quot;
  123. from pprint import pprint
  124. from typing import NamedTuple
  126. from numpy.random import default_rng
  127. from scipy.optimize import least_squares
  128. import numpy as np
  131. class Vars(NamedTuple):
  132.     e0_K_tr1_c2: float            # 0
  133.     e0_x_Lint_tr1_c2: float       # 1
  134.     e0_h_L_tr1: float             # 2
  135.     e0_K_tr1_c1: float            # 3
  136.     e0_x_Lint_tr1_c1: float       # 4
  137.     e0_lambda_V_tr1: float        # 5
  138.     e0_cp_V_tr1: float            # 6
  139.     e0_f_int_tr1_c2: float        # 7
  140.     e0_f_int_tr1_c1: float        # 8
  141.     e0_x_L_tr1_c2: float          # 9
  142.     e0_x_L_tr1_c1: float          # 10
  143.     e0_cp_L_tr1: float            # 11
  144.     e0_h_V_tr1: float             # 12
  145.     e0_D_L_tr1: float             # 13
  146.     e0_F_V_tr1: float             # 14
  147.     e0_D_V_tr1: float             # 15
  148.     e0_beta_L_tr1: float          # 16
  149.     e0_e_int_tr1: float           # 17
  150.     e0_alpha_L_tr1: float         # 18
  151.     e0_rho_L_tr1: float           # 19
  152.     e0_x_V_tr1_c2: float          # 20
  153.     e0_lambda_L_tr1: float        # 21
  154.     e0_x_V_tr1_c1: float          # 22
  155.     e0_F_L_tr1: float             # 23
  156.     e0_rho_V_tr1: float           # 24
  157.     e0_beta_V_tr1: float          # 25
  158.     e0_T_V_tr1: float             # 26
  159.     e0_alpha_V_tr1: float         # 27
  160.     e0_T_int_tr1: float           # 28
  161.     e0_T_L_tr1: float             # 29
  162.     e0_x_Vint_tr1_c2: float       # 30
  163.     e0_x_Vint_tr1_c1: float       # 31
  166. def equations(x: np.ndarray) -&gt; tuple[float, ...]:
  167.     v = Vars(*x)
  169.     eq1 = v.e0_alpha_L_tr1 - v.e0_beta_L_tr1 * v.e0_rho_L_tr1 * v.e0_cp_L_tr1 * (
  170.         v.e0_lambda_L_tr1 / (v.e0_rho_L_tr1 * v.e0_cp_L_tr1 * v.e0_D_L_tr1)
  171.     ) ** (2 / 3)
  173.     eq2 = v.e0_rho_L_tr1 - (732 * v.e0_x_L_tr1_c1 + 747 * v.e0_x_L_tr1_c2)
  174.     eq3 = v.e0_D_L_tr1 - (1.25E-5 * v.e0_x_L_tr1_c1 + 1.6E-6 * v.e0_x_L_tr1_c2)
  176.     eq4 = v.e0_h_L_tr1 - (
  177.         (
  178.             -6021.54
  179.             + 0.055134232180918 * (v.e0_T_L_tr1 - 298)
  180.             + 2.61275786690816
  181.         ) * v.e0_x_L_tr1_c1
  182.         + (
  183.             -6883.827
  184.             + 0.058471271259418 * (v.e0_T_L_tr1 - 298)
  185.             + 2.38446400137128
  186.         ) * v.e0_x_L_tr1_c2
  187.     )
  189.     eq5 = v.e0_cp_L_tr1 - (
  190.         (
  191.             0.055134232180918 * v.e0_T_L_tr1 + 2.61275786690816
  192.         ) * v.e0_x_L_tr1_c1
  193.         + (
  194.             0.058471271259418 * v.e0_T_L_tr1 + 2.38446400137128
  195.         ) * v.e0_x_L_tr1_c2
  196.     )
  197.     eq6 = v.e0_lambda_L_tr1 - (
  198.         (
  199.             0.001318 * v.e0_T_L_tr1 + 0.16762
  200.         ) * v.e0_x_L_tr1_c1
  201.         + (
  202.             -0.001075424413 * v.e0_T_L_tr1 + 0.155701011727667
  203.         ) * v.e0_x_L_tr1_c2
  204.     )
  206.     eq7 = v.e0_rho_V_tr1 - (1.746 * v.e0_x_V_tr1_c1 + 2.76 * v.e0_x_V_tr1_c2)
  207.     eq8 = v.e0_D_V_tr1 - (1E-5 * v.e0_x_V_tr1_c1 + 6.5E-6 * v.e0_x_V_tr1_c2)
  209.     eq9 = v.e0_h_V_tr1 - (
  210.         (
  211.             -6021.54
  212.             + 0.015220202121093 * (v.e0_T_V_tr1 - 298)
  213.             + 1.60395800121061
  214.             + 922.194
  215.         ) * v.e0_x_V_tr1_c1
  216.         + (
  217.             -6883.827
  218.             + 0.01847750929725 * (v.e0_T_V_tr1 - 298)
  219.             + 1.43843773750579
  220.             + 1067.523
  221.         ) * v.e0_x_V_tr1_c2
  222.     )
  223.     eq10 = v.e0_cp_V_tr1 - (
  224.         (
  225.             0.015220202121093 * v.e0_T_V_tr1 + 1.60395800121061
  226.         ) * v.e0_x_V_tr1_c1
  227.         + (
  228.             0.01847750929725 * v.e0_T_V_tr1 + 1.43843773750579
  229.         ) * v.e0_x_V_tr1_c2
  230.     )
  231.     eq11 = v.e0_lambda_V_tr1 - (
  232.         (
  233.             5.07392141746E-4 * v.e0_T_L_tr1 + 0.015070468544221
  234.         ) * v.e0_x_V_tr1_c1
  235.         + (
  236.             4.89248E-4 * v.e0_T_L_tr1 + 0.014354384
  237.         ) * v.e0_x_V_tr1_c2
  238.     )
  240.     eq12 = v.e0_K_tr1_c1 - 1 / 10**(5.24677 - 1598.673 / (v.e0_T_int_tr1 + 46.424))
  241.     eq13 = v.e0_K_tr1_c2 - 1 / 10**(5.31384 - 1690.864 / (v.e0_T_int_tr1 + 51.804))
  243.     eq14 = v.e0_alpha_V_tr1 - (
  244.         v.e0_beta_V_tr1 * v.e0_rho_V_tr1 * v.e0_cp_V_tr1 * (
  245.             v.e0_lambda_V_tr1 / (
  246.                 v.e0_rho_V_tr1 * v.e0_cp_V_tr1 * v.e0_D_V_tr1
  247.             )
  248.         ) ** (2 / 3)
  249.     )
  251.     eq15 = -(2 * 28000 - v.e0_F_L_tr1 * v.e0_h_L_tr1 + v.e0_e_int_tr1 * 0.002)
  252.     eq16 = -(0.2 * 30000 - v.e0_F_V_tr1 * v.e0_h_V_tr1 - v.e0_e_int_tr1 * 0.002)
  253.     eq17 = -(2 * 0.5 - v.e0_F_L_tr1 * v.e0_x_L_tr1_c1 + v.e0_f_int_tr1_c1 * 0.002)
  254.     eq18 = -(2 * 0.5 - v.e0_F_L_tr1 * v.e0_x_L_tr1_c2 + v.e0_f_int_tr1_c2 * 0.002)
  255.     eq19 = -(0.2 * 0.6 - v.e0_F_V_tr1 * v.e0_x_V_tr1_c1 - v.e0_f_int_tr1_c1 * 0.002)
  256.     eq20 = -(0.2 * 0.4 - v.e0_F_V_tr1 * v.e0_x_V_tr1_c2 - v.e0_f_int_tr1_c2 * 0.002)
  258.     eq21 = -(v.e0_x_Vint_tr1_c1 - v.e0_K_tr1_c1 * v.e0_x_Lint_tr1_c1)
  259.     eq22 = -(v.e0_x_Vint_tr1_c2 - v.e0_K_tr1_c2 * v.e0_x_Lint_tr1_c2)
  261.     eq23 = v.e0_e_int_tr1 - (
  262.         v.e0_alpha_L_tr1 * (v.e0_T_int_tr1 - v.e0_T_L_tr1)
  263.         + v.e0_f_int_tr1_c1 * (
  264.             -6021.54
  265.             + 0.055134232180918 * (v.e0_T_L_tr1 - 298)
  266.             + 2.61275786690816
  267.             # really?
  268.         # )
  269.             + v.e0_f_int_tr1_c2 * (
  270.                 -6883.827 +
  271.                 0.058471271259418 * (v.e0_T_L_tr1 - 298)
  272.                 + 2.38446400137128
  273.             )
  274.         )
  275.     )
  277.     eq24 = v.e0_e_int_tr1 - (
  278.         v.e0_alpha_V_tr1 * (v.e0_T_V_tr1 - v.e0_T_int_tr1)
  279.         + v.e0_f_int_tr1_c1 * (
  280.             -6021.54
  281.             + 0.015220202121093 * (v.e0_T_V_tr1 - 298)
  282.             + 1.60395800121061
  283.             + 922.194
  284.             # really?
  285.         # )
  286.             + v.e0_f_int_tr1_c2 * (
  287.                 -6883.827
  288.                 + 0.01847750929725 * (v.e0_T_V_tr1 - 298)
  289.                 + 1.43843773750579
  290.                 + 1067.523
  291.             )
  292.         )
  293.     )
  295.     eq25 = v.e0_f_int_tr1_c1 - (
  296.         v.e0_rho_L_tr1 * v.e0_beta_L_tr1 * (
  297.             v.e0_x_Lint_tr1_c1 - v.e0_x_L_tr1_c1
  298.         )
  299.         + v.e0_x_L_tr1_c1 * (v.e0_f_int_tr1_c1 + v.e0_f_int_tr1_c2)
  300.     )
  301.     eq26 = v.e0_f_int_tr1_c2 - (
  302.         v.e0_rho_L_tr1 * v.e0_beta_L_tr1 * (
  303.             v.e0_x_Lint_tr1_c2 - v.e0_x_L_tr1_c2
  304.         )
  305.         + v.e0_x_L_tr1_c2 * (v.e0_f_int_tr1_c1 + v.e0_f_int_tr1_c2)
  306.     )
  307.     eq27 = v.e0_f_int_tr1_c1 - (
  308.         v.e0_rho_V_tr1 * v.e0_beta_V_tr1 * (
  309.             v.e0_x_V_tr1_c1 - v.e0_x_Vint_tr1_c1
  310.         )
  311.         + v.e0_x_V_tr1_c1 * (v.e0_f_int_tr1_c1 + v.e0_f_int_tr1_c2)
  312.     )
  313.     eq28 = v.e0_f_int_tr1_c2 - (
  314.         v.e0_rho_V_tr1 * v.e0_beta_V_tr1 * (
  315.             v.e0_x_V_tr1_c2 - v.e0_x_Vint_tr1_c2
  316.         )
  317.         + v.e0_x_V_tr1_c2 * (v.e0_f_int_tr1_c1 + v.e0_f_int_tr1_c2)
  318.     )
  320.     eq29 = 1 - (v.e0_x_L_tr1_c1 + v.e0_x_L_tr1_c2)
  321.     eq30 = 1 - (v.e0_x_Lint_tr1_c1 + v.e0_x_Lint_tr1_c2)
  322.     eq31 = 1 - (v.e0_x_V_tr1_c1 + v.e0_x_V_tr1_c2)
  323.     eq32 = 1 - (v.e0_x_Vint_tr1_c1 + v.e0_x_Vint_tr1_c2)
  325.     return (eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10,
  326.             eq11, eq12, eq13, eq14, eq15, eq16, eq17, eq18, eq19, eq20,
  327.             eq21, eq22, eq23, eq24, eq25, eq26, eq27, eq28, eq29, eq30,
  328.             eq31, eq32)
  331. def regression_test(x0: np.ndarray) -&gt; None:
  332.     y = equations(x0)
  333.     assert np.allclose(
  334.         y,
  335.         [ 0.00000000e+00, -7.38500000e+02,  9.99992950e-01,  6.44709509e+03,
  336.          -2.25156291e+01,  7.93463011e-01, -1.25300000e+00,  9.99991750e-01,
  337.           5.45609068e+03, -6.75527448e+00,  8.00909148e-01,  9.60890608e-01,
  338.           9.50476815e-01,  0.00000000e+00, -5.59990020e+04, -5.99899800e+03,
  339.          -5.02000000e-01, -5.02000000e-01,  3.82000000e-01,  4.22000000e-01,
  340.           0.00000000e+00,  0.00000000e+00,  1.28931902e+04,  1.09111814e+04,
  341.           0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
  342.           0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
  343.           ],
  344.         rtol=0, atol=1e-4,
  345.     )
  347.     rand = default_rng(seed=0)
  348.     x = rand.uniform(low=0, high=2, size=32)
  349.     y = equations(x)
  350.     assert np.allclose(
  351.         y,
  352.         [-8.89664807e-02, -2.59056243e+03,  6.71477620e-02,  2.27496385e+04,
  353.          -8.97636119e+00, -3.16304973e-01, -1.26737482e+00,  3.51297461e-01,
  354.           7.17475798e+03, -1.06230892e+00,  1.80356351e+00, -1.59115277e+28,
  355.          -3.04937193e+26, -6.69625188e-01, -5.59998961e+04, -5.99749540e+03,
  356.           1.10987295e+00,  1.41775827e+00,  1.83947353e+00,  5.57239541e-03,
  357.          -7.24077112e-01, -6.89518259e-01,  1.75063897e+04,  1.47800191e+04,
  358.          -3.05993271e+00, -1.36098936e+00, -2.85999752e+00,  2.56167351e+00,
  359.          -2.50185196e+00, -1.16611391e+00, -3.97888172e-01, -1.15473631e+00,
  360.          ],
  361.         rtol=1e-8, atol=0,
  362.     )
  365. def make_x0() -&gt; np.ndarray:
  366.     x0 = np.ones(32)
  367.     x0[[1,4,9,10,20,22,30,31]] = 0.5
  368.     x0[[26,28,29]] = 370
  369.     return x0
  372. def solve(x0: np.ndarray) -&gt; Vars:
  373.     upper_bounds = np.full(shape=32, fill_value=np.inf)
  374.     upper_bounds[[1,4,9,10,20,22,30,31]] = 1
  376.     bottom_bounds = np.zeros(32)
  377.     bottom_bounds[[26,28,29]] = 340
  379.     res = least_squares(
  380.         fun=equations, x0=x0,
  381.         bounds=(bottom_bounds, upper_bounds),  # not supported in lm
  382.         method=&#39;trf&#39;,
  383.         ftol=1e-10, xtol=1e-10, gtol=1e-10,
  384.     )
  385.     assert res.success
  386.     return Vars(*res.x)
  389. def dump(x: Vars) -&gt; None:
  390.     pprint(x._asdict())
  393. def main() -&gt; None:
  394.     x0 = make_x0()
  395.     regression_test(x0)
  396.     x = solve(x0)
  397.     dump(x)
  400. if __name__ == &#39;__main__&#39;:
  401.     main()
  1. {&#39;e0_D_L_tr1&#39;: 0.9984442691049119,
  2. &#39;e0_D_V_tr1&#39;: 0.9986041251126907,
  3. &#39;e0_F_L_tr1&#39;: 128.59184097299763,
  4. &#39;e0_F_V_tr1&#39;: 194.05957359631154,
  5. &#39;e0_K_tr1_c1&#39;: 0.9994676910341592,
  6. &#39;e0_K_tr1_c2&#39;: 0.999729237718534,
  7. &#39;e0_T_L_tr1&#39;: 340.0,
  8. &#39;e0_T_V_tr1&#39;: 559.936902192616,
  9. &#39;e0_T_int_tr1&#39;: 369.2496832286959,
  10. &#39;e0_alpha_L_tr1&#39;: 7.5810026415896,
  11. &#39;e0_alpha_V_tr1&#39;: 30.124493018891847,
  12. &#39;e0_beta_L_tr1&#39;: 0.0,
  13. &#39;e0_beta_V_tr1&#39;: 1.2479650000565898,
  14. &#39;e0_cp_L_tr1&#39;: 5.468269302698259,
  15. &#39;e0_cp_V_tr1&#39;: 1.0196495285046383,
  16. &#39;e0_e_int_tr1&#39;: 1856.5186325404206,
  17. &#39;e0_f_int_tr1_c1&#39;: 0.11670635205679133,
  18. &#39;e0_f_int_tr1_c2&#39;: 0.39727951871442857,
  19. &#39;e0_h_L_tr1&#39;: 190.08447743399068,
  20. &#39;e0_h_V_tr1&#39;: 22.650731758245296,
  21. &#39;e0_lambda_L_tr1&#39;: 1.0083185179925418,
  22. &#39;e0_lambda_V_tr1&#39;: 1.0368032854477796,
  23. &#39;e0_rho_L_tr1&#39;: 1.0430730156966466,
  24. &#39;e0_rho_V_tr1&#39;: 0.9986707757704435,
  25. &#39;e0_x_L_tr1_c1&#39;: 0.0,
  26. &#39;e0_x_L_tr1_c2&#39;: 0.0,
  27. &#39;e0_x_Lint_tr1_c1&#39;: 0.5001045740608061,
  28. &#39;e0_x_Lint_tr1_c2&#39;: 0.49992448934879186,
  29. &#39;e0_x_V_tr1_c1&#39;: 1.4663484981775705e-16,
  30. &#39;e0_x_V_tr1_c2&#39;: 0.00041819098423356144,
  31. &#39;e0_x_Vint_tr1_c1&#39;: 0.4998465855173379,
  32. &#39;e0_x_Vint_tr1_c2&#39;: 0.4995202978525341}

huangapple
  • 本文由 发表于 2023年7月27日 21:44:41
  • 转载请务必保留本文链接：https://go.coder-hub.com/76780385.html
  • nonlinear-optimization
  • optimization
  • python
  • scipy
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定