非线性方程组带约束条件

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

Nonlinear equation system with constraints

问题

  1. from scipy.optimize import least_squares
  2. import numpy as np
  4. def equations(x):
  5.     eq1 = x[7] - 1 / (np.exp((5.24677 - 1598.673 / (60.0 - -46.424)) * np.log(10.0)))
  6.     eq2 = x[4] - 1 / (np.exp((5.08354 - 1663.125 / (60.0 - -45.622)) * np.log(10.0)))
  7.     eq3 = 0.0 - (0.5 * 0.4) + x[2] * x[9] - x[14] * (0.02 * 0.1)
  8.     eq4 = 0.0 - (0.5 * 0.6) + x[2] * x[5] - x[12] * (0.02 * 0.1)
  9.     eq5 = 0.0 - (0.4 * 0.6) + x[3] * x[0] + x[14] * (0.02 * 0.1)
  10.     eq6 = 0.0 - (0.4 * 0.4) + x[3] * x[15] + x[12] * (0.02 * 0.1)
  11.     eq7 = 0.0 - (x[10] - x[7] * x[13])
  12.     eq8 = 0.0 - (x[8] - x[4] * x[11])
  13.     eq9 = x[14] - (1.0 * x[6] * (x[13] - x[9]) + x[9] * (x[14] + x[12]))
  14.     eq10 = x[12] - (1.0 * x[6] * (x[11] - x[5]) + x[5] * (x[14] + x[12]))
  15.     eq11 = x[14] - (0.01 * x[1] * (x[0] - x[10]) + x[0] * (x[14] + x[12]))
  16.     eq12 = x[12] - (0.01 * x[1] * (x[15] - x[8]) + x[15] * (x[14] + x[12]))
  17.     eq13 = 1.0 - (x[9] + x[5])
  18.     eq14 = 1.0 - (x[0] + x[15])
  19.     eq15 = 1.0 - (x[13] + x[11])
  20.     eq16 = 1.0 - (x[10] + x[8])
  21.     return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16]
  23. # constraints
  24. bounds = (np.zeros(16), np.array([1, np.inf, np.inf, np.inf, 1, np.inf, np.inf, 1, 1, 1, 1, 1, np.inf, 1, np.inf, 1]))
  26. # inital 
  27. x0 = np.zeros(16)
  28. # x0 = np.array([1, 1e-4, 1, 1, 1, 1, 1e-4, 1, 1, 1, 1, 1, 1, 1, 1, 1])
  30. res = least_squares(equations, x0, bounds=bounds)
  32. print(res.x)
英文:

Hello I want to solve a thermodynamic equation system. Unfortunately, I am quite a beginner in Python and need some help. I tried to solve it with the least_squares method but I always got my initial vector as a result.

I hope someone can help!

  1. from scipy.optimize import least_squares
  2. import numpy as np
  3. def equations(x):
  4. eq1 = x[7] - 1 / (np.exp((5.24677 - 1598.673 / (60.0 - -46.424)) * np.log(10.0)))
  5. eq2 = x[4] - 1 / (np.exp((5.08354 - 1663.125 / (60.0 - -45.622)) * np.log(10.0)))
  6. eq3 = 0.0 - (0.5 * 0.4) + x[2] * x[9] - x[14] * (0.02 * 0.1)
  7. eq4 = 0.0 - (0.5 * 0.6) + x[2] * x[5] - x[12] * (0.02 * 0.1)
  8. eq5 = 0.0 - (0.4 * 0.6) + x[3] * x[0] + x[14] * (0.02 * 0.1)
  9. eq6 = 0.0 - (0.4 * 0.4) + x[3] * x[15] + x[12] * (0.02 * 0.1)
  10. eq7 = 0.0 - (x[10] - x[7] * x[13])
  11. eq8 = 0.0 - (x[8] - x[4] * x[11])
  12. eq9 = x[14] - (1.0 * x[6] * (x[13] - x[9]) + x[9] * (x[14] + x[12]))
  13. eq10 = x[12] - (1.0 * x[6] * (x[11] - x[5]) + x[5] * (x[14] + x[12]))
  14. eq11 = x[14] - (0.01 * x[1] * (x[0] - x[10]) + x[0] * (x[14] + x[12]))
  15. eq12 = x[12] - (0.01 * x[1] * (x[15] - x[8]) + x[15] * (x[14] + x[12]))
  16. eq13 = 1.0 - (x[9] + x[5])
  17. eq14 = 1.0 - (x[0] + x[15])
  18. eq15 = 1.0 - (x[13] + x[11])
  19. eq16 = 1.0 - (x[10] + x[8])
  20. return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16]
  21. # constraints
  22. bounds = (np.zeros(16), np.array([1, np.inf, np.inf, np.inf, 1, np.inf, np.inf, 1, 1, 1, 1, 1, np.inf, 1, np.inf, 1]))
  23. # inital 
  24. x0 = np.zeros(16)
  25. # x0 = np.array([1, 1e-4, 1, 1, 1, 1, 1e-4, 1, 1, 1, 1, 1, 1, 1, 1, 1])
  26. res = least_squares(equations, x0, bounds=bounds)
  27. print(res.x)

答案1

得分: 1

以下是已经翻译好的部分：

"One probable problem is that the scale of the first and second RHS is enormous compared to the others. x4 and x7 appear to be fixed, but you do not reflect this in your initial values. And indeed, the upper bound for x4 and x7 is 1, which makes absolutely no sense given the first and second equations. So that needs to be fixed for your system of equations to be sensible.

I show below that disregarding the first two equations (reordered to be equations 4 and 7) produces a seemingly reasonable solution."

"这里可能存在的问题是，与其他方程相比，第一个和第二个右侧的比例巨大。x4和x7似乎是固定的，但您在初始值中没有反映这一点。而且，x4和x7的上限是1，这在考虑第一个和第二个方程时毫无意义。因此，必须修正这个问题，使您的方程系统变得合理。

我在下面展示，忽略了前两个方程（重新排列为方程4和7），得到了一个看似合理的解决方案。"

接下来的部分包括Python代码和输出，您可以根据需要研究它们。如果您需要更多翻译，请告诉我需要翻译的具体内容。

英文:

One probable problem is that the scale of the first and second RHS is enormous compared to the others. x4 and x7 appear to be fixed, but you do not reflect this in your initial values. And indeed, the upper bound for x4 and x7 is 1, which makes absolutely no sense given the first and second equations. So that needs to be fixed for your system of equations to be sensible.

I show below that disregarding the first two equations (reordered to be equations 4 and 7) produces a seemingly reasonable solution.

  1. from scipy.optimize import least_squares
  2. import numpy as np
  5. def sides(x: np.ndarray) -> np.ndarray:
  6.     # x order as original, y reordered so that x indices are closer to their use
  7.     return np.array((
  8.         (x[14], 0.01*x[1]*(x[ 0] - x[10]) + x[ 0]*(x[14] + x[12])),
  9.         (x[12], 0.01*x[1]*(x[15] - x[ 8]) + x[15]*(x[14] + x[12])),
  10.         (0,     -0.5*0.6 + x[2]* x[ 5] - x[12]*(0.02 * 0.1)),
  11.         (0,     -0.5*0.4 + x[2]* x[ 9] - x[14]*(0.02 * 0.1)),
  12.         (x[4],  1 / np.exp((5.08354 - 1663.125/(60.0 - -45.622)) * np.log(10))),
  13.         (x[12], x[6]*(x[11] - x[5]) + x[5]*(x[14] + x[12])),
  14.         (x[14], x[6]*(x[13] - x[9]) + x[9]*(x[14] + x[12])),
  15.         (x[7],  1 / np.exp((5.24677 - 1598.673/(60.0 - -46.424)) * np.log(10))),
  16.         (1,     x[10] + x[ 8]),
  17.         (1,     x[ 9] + x[ 5]),
  18.         (0,     x[10] - x[7]*x[13]),
  19.         (0,     x[ 8] - x[4]*x[11]),
  20.         (0,     -0.4*0.4 + x[3]* x[15] + x[12]*(0.02 * 0.1)),
  21.         (1,     x[13] + x[11]),
  22.         (0,     -0.4*0.6 + x[3]* x[ 0] + x[14]*(0.02 * 0.1)),
  23.         (1,     x[ 0] + x[15]),
  24.     ))
  27. def equations(x: np.ndarray) -> np.ndarray:
  28.     left, right = sides(x)[[
  29.         0, 1, 2, 3,
  30.         5, 6,
  31.         8, 9, 10, 11, 12, 13, 14, 15,
  32.     ], :].T
  33.     return left - right
  36. # Ordered as original
  37. x0 = np.array((
  38.     1, 1e-4, 1, 1,
  39.     10,  # 1 / np.exp((5.08354 - 1663.125/(60.0 - -45.622)) * np.log(10)),
  40.     1, 1e-4,
  41.     10,  # 1 / np.exp((5.24677 - 1598.673/(60.0 - -46.424)) * np.log(10)),
  42.     1, 1, 1, 1, 1, 1, 1, 1,
  43. ))
  45. # Ordered as original, with x[4] and x[7] bounds changed to inf (not fixed)
  46. bounds = np.array((
  47.     np.zeros(16),
  48.     (
  49.         1, np.inf, np.inf, np.inf,
  50.         np.inf,  # 1,
  51.         np.inf, np.inf,
  52.         np.inf,  # 1,
  53.         1, 1, 1, 1, np.inf, 1, np.inf, 1,
  54.     ),
  55. ))
  57. res = least_squares(fun=equations, x0=x0, bounds=bounds)
  58. assert res.success, res.message
  60. for i, (x, (left, right)) in enumerate(zip(res.x, sides(res.x))):
  61.     print(f'y{i:<2} = {left:8.3g} ~ {right:9.3g}, '
  62.           f'err={right-left:8.1e}, '
  63.           f'x{i:<2} = {x:.3g}')
  1. y0  =     6.25 ~      6.25, err= 5.3e-15, x0  = 0.591
  2. y1  =     1.41 ~      1.41, err= 0.0e+00, x1  = 361
  3. y2  =        0 ~ -1.25e-13, err=-1.3e-13, x2  = 0.515
  4. y3  =        0 ~  1.78e-13, err= 1.8e-13, x3  = 0.385
  5. y4  =     1.56 ~   4.6e+10, err= 4.6e+10, x4  = 1.56
  6. y5  =     1.41 ~      1.41, err=-6.9e-15, x5  = 0.588
  7. y6  =     6.25 ~      6.25, err= 1.2e-14, x6  = 139
  8. y7  =    0.266 ~  5.96e+09, err= 6.0e+09, x7  = 0.266
  9. y8  =        1 ~         1, err=-2.2e-15, x8  = 0.884
  10. y9  =        1 ~         1, err= 1.9e-13, x9  = 0.412
  11. y10 =        0 ~  4.16e-17, err= 4.2e-17, x10 = 0.116
  12. y11 =        0 ~  2.22e-16, err= 2.2e-16, x11 = 0.565
  13. y12 =        0 ~  1.07e-13, err= 1.1e-13, x12 = 1.41
  14. y13 =        1 ~         1, err= 1.8e-13, x13 = 0.435
  15. y14 =        0 ~ -7.43e-14, err=-7.4e-14, x14 = 6.25
  16. y15 =        1 ~         1, err=-2.2e-16, x15 = 0.409

This sidekick program will show you some options where values in the two problematic equations could be substituted to produce a reasonable solution in the original problem:

  1. import numpy as np
  2. from sympy import Symbol, Eq, solveset
  4. for y, values in (
  5.     (1.690, (5.08354, 1663.125, 60, -45.622)),
  6.     (0.138, (5.24677, 1598.673, 60, -46.424)),
  7. ):
  8.     a = Symbol('a', real=True, infinite=False)
  9.     b = Symbol('b', real=True, infinite=False)
  10.     c = Symbol('c', real=True, infinite=False)
  11.     d = Symbol('d', real=True, infinite=False)
  12.     variables = a, b, c, d
  14.     eq = Eq(-np.log10(y), a - b/(c - d))
  16.     print('Do one of:')
  17.     for target, old_val in zip(variables, values):
  18.         seq = eq.subs({
  19.             s: v for s, v in zip(variables, values)
  20.             if s is not target
  21.         })
  22.         replacement, = solveset(seq, target)
  24.         print(f'Replace {target}={old_val} with {replacement:.1f}')
  25.     print()
  1. Do one of:
  2. Replace a=5.08354 with 15.5
  3. Replace b=1663.125 with 561.0
  4. Replace c=60 with 267.5
  5. Replace d=-45.622 with -253.1
  6. Do one of:
  7. Replace a=5.24677 with 15.9
  8. Replace b=1598.673 with 466.8
  9. Replace c=60 with 318.0
  10. Replace d=-46.424 with -304.4

huangapple
  • 本文由 发表于 2023年7月3日 20:13:18
  • 转载请务必保留本文链接：https://go.coder-hub.com/76604636.html
  • constraints
  • minimize
  • optimization
  • scipy
匿名

发表评论

匿名网友

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

确定