英文:
Prioritize equations in a nonlinear system
问题
我正在尝试使用最小二乘法解决非线性等式系统（EQS）的问题。我得到了一个结果，但不幸的是它并不非常准确。主要问题是方程29到32必须被满足，但它们并没有被满足。有没有一种方法可以在非线性系统中优先考虑方程？
以下是我的代码：
# 乙醇和1-丙醇与3个塔板
from scipy.optimize import least_squares
import numpy as np
def equations(x):
# 你的方程式在这里
# 界限
# 特定变量的上界
upper_bounds = np.ones(32) * np.inf
# 你的上界在这里
# 特定变量的下界
bottom_bounds = np.zeros(32)
# 你的下界在这里
# 边界数组
bounds = (bottom_bounds, upper_bounds)
# 初始值
x0 = np.ones(32)
# 你的初始值在这里
# 使用变量边界解决方程组
res = least_squares(equations, x0, bounds=bounds, method='trf', ftol=1e-10, xtol=1e-10, gtol=1e-10)
# 打印结果
print(res.x)
# 你的结果打印在这里
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:
# 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]))
and the results for x[31], x[30], x[22],... are:
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
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-丙醇，3个托盘"""
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
<details>
<summary>英文:</summary>
Remove your redundant parentheses and reformat your equations so that they can be more easily debugged.
Replace your indexed variables with named variables. I show a quick way to do this via `NamedTuple`.
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's not usable because it doesn't respect your bounds. There is no magic bullet to improve the accuracy. Knowledge of the domain problem is necessary. If some equations "matter more" to the solution than others, you can simply scale those equations by their weights.
Add regression tests. For the results that I show, the equations are equivalent to your original equations.
Always check `result.success`.
A first pass can look like
```python
"""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()
{'e0_D_L_tr1': 0.9984442691049119,
'e0_D_V_tr1': 0.9986041251126907,
'e0_F_L_tr1': 128.59184097299763,
'e0_F_V_tr1': 194.05957359631154,
'e0_K_tr1_c1': 0.9994676910341592,
'e0_K_tr1_c2': 0.999729237718534,
'e0_T_L_tr1': 340.0,
'e0_T_V_tr1': 559.936902192616,
'e0_T_int_tr1': 369.2496832286959,
'e0_alpha_L_tr1': 7.5810026415896,
'e0_alpha_V_tr1': 30.124493018891847,
'e0_beta_L_tr1': 0.0,
'e0_beta_V_tr1': 1.2479650000565898,
'e0_cp_L_tr1': 5.468269302698259,
'e0_cp_V_tr1': 1.0196495285046383,
'e0_e_int_tr1': 1856.5186325404206,
'e0_f_int_tr1_c1': 0.11670635205679133,
'e0_f_int_tr1_c2': 0.39727951871442857,
'e0_h_L_tr1': 190.08447743399068,
'e0_h_V_tr1': 22.650731758245296,
'e0_lambda_L_tr1': 1.0083185179925418,
'e0_lambda_V_tr1': 1.0368032854477796,
'e0_rho_L_tr1': 1.0430730156966466,
'e0_rho_V_tr1': 0.9986707757704435,
'e0_x_L_tr1_c1': 0.0,
'e0_x_L_tr1_c2': 0.0,
'e0_x_Lint_tr1_c1': 0.5001045740608061,
'e0_x_Lint_tr1_c2': 0.49992448934879186,
'e0_x_V_tr1_c1': 1.4663484981775705e-16,
'e0_x_V_tr1_c2': 0.00041819098423356144,
'e0_x_Vint_tr1_c1': 0.4998465855173379,
'e0_x_Vint_tr1_c2': 0.4995202978525341}
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论