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!
from scipy.optimize import least_squares
import numpy as np
def equations(x):
eq1 = x[7] - 1 / (np.exp((5.24677 - 1598.673 / (60.0 - -46.424)) * np.log(10.0)))
eq2 = x[4] - 1 / (np.exp((5.08354 - 1663.125 / (60.0 - -45.622)) * np.log(10.0)))
eq3 = 0.0 - (0.5 * 0.4) + x[2] * x[9] - x[14] * (0.02 * 0.1)
eq4 = 0.0 - (0.5 * 0.6) + x[2] * x[5] - x[12] * (0.02 * 0.1)
eq5 = 0.0 - (0.4 * 0.6) + x[3] * x[0] + x[14] * (0.02 * 0.1)
eq6 = 0.0 - (0.4 * 0.4) + x[3] * x[15] + x[12] * (0.02 * 0.1)
eq7 = 0.0 - (x[10] - x[7] * x[13])
eq8 = 0.0 - (x[8] - x[4] * x[11])
eq9 = x[14] - (1.0 * x[6] * (x[13] - x[9]) + x[9] * (x[14] + x[12]))
eq10 = x[12] - (1.0 * x[6] * (x[11] - x[5]) + x[5] * (x[14] + x[12]))
eq11 = x[14] - (0.01 * x[1] * (x[0] - x[10]) + x[0] * (x[14] + x[12]))
eq12 = x[12] - (0.01 * x[1] * (x[15] - x[8]) + x[15] * (x[14] + x[12]))
eq13 = 1.0 - (x[9] + x[5])
eq14 = 1.0 - (x[0] + x[15])
eq15 = 1.0 - (x[13] + x[11])
eq16 = 1.0 - (x[10] + x[8])
return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16]
# constraints
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]))
# inital
x0 = np.zeros(16)
# x0 = np.array([1, 1e-4, 1, 1, 1, 1, 1e-4, 1, 1, 1, 1, 1, 1, 1, 1, 1])
res = least_squares(equations, x0, bounds=bounds)
print(res.x)
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.
from scipy.optimize import least_squares
import numpy as np
def sides(x: np.ndarray) -> np.ndarray:
# x order as original, y reordered so that x indices are closer to their use
return np.array((
(x[14], 0.01*x[1]*(x[ 0] - x[10]) + x[ 0]*(x[14] + x[12])),
(x[12], 0.01*x[1]*(x[15] - x[ 8]) + x[15]*(x[14] + x[12])),
(0, -0.5*0.6 + x[2]* x[ 5] - x[12]*(0.02 * 0.1)),
(0, -0.5*0.4 + x[2]* x[ 9] - x[14]*(0.02 * 0.1)),
(x[4], 1 / np.exp((5.08354 - 1663.125/(60.0 - -45.622)) * np.log(10))),
(x[12], x[6]*(x[11] - x[5]) + x[5]*(x[14] + x[12])),
(x[14], x[6]*(x[13] - x[9]) + x[9]*(x[14] + x[12])),
(x[7], 1 / np.exp((5.24677 - 1598.673/(60.0 - -46.424)) * np.log(10))),
(1, x[10] + x[ 8]),
(1, x[ 9] + x[ 5]),
(0, x[10] - x[7]*x[13]),
(0, x[ 8] - x[4]*x[11]),
(0, -0.4*0.4 + x[3]* x[15] + x[12]*(0.02 * 0.1)),
(1, x[13] + x[11]),
(0, -0.4*0.6 + x[3]* x[ 0] + x[14]*(0.02 * 0.1)),
(1, x[ 0] + x[15]),
))
def equations(x: np.ndarray) -> np.ndarray:
left, right = sides(x)[[
0, 1, 2, 3,
5, 6,
8, 9, 10, 11, 12, 13, 14, 15,
], :].T
return left - right
# Ordered as original
x0 = np.array((
1, 1e-4, 1, 1,
10, # 1 / np.exp((5.08354 - 1663.125/(60.0 - -45.622)) * np.log(10)),
1, 1e-4,
10, # 1 / np.exp((5.24677 - 1598.673/(60.0 - -46.424)) * np.log(10)),
1, 1, 1, 1, 1, 1, 1, 1,
))
# Ordered as original, with x[4] and x[7] bounds changed to inf (not fixed)
bounds = np.array((
np.zeros(16),
(
1, np.inf, np.inf, np.inf,
np.inf, # 1,
np.inf, np.inf,
np.inf, # 1,
1, 1, 1, 1, np.inf, 1, np.inf, 1,
),
))
res = least_squares(fun=equations, x0=x0, bounds=bounds)
assert res.success, res.message
for i, (x, (left, right)) in enumerate(zip(res.x, sides(res.x))):
print(f'y{i:<2} = {left:8.3g} ~ {right:9.3g}, '
f'err={right-left:8.1e}, '
f'x{i:<2} = {x:.3g}')
y0 = 6.25 ~ 6.25, err= 5.3e-15, x0 = 0.591
y1 = 1.41 ~ 1.41, err= 0.0e+00, x1 = 361
y2 = 0 ~ -1.25e-13, err=-1.3e-13, x2 = 0.515
y3 = 0 ~ 1.78e-13, err= 1.8e-13, x3 = 0.385
y4 = 1.56 ~ 4.6e+10, err= 4.6e+10, x4 = 1.56
y5 = 1.41 ~ 1.41, err=-6.9e-15, x5 = 0.588
y6 = 6.25 ~ 6.25, err= 1.2e-14, x6 = 139
y7 = 0.266 ~ 5.96e+09, err= 6.0e+09, x7 = 0.266
y8 = 1 ~ 1, err=-2.2e-15, x8 = 0.884
y9 = 1 ~ 1, err= 1.9e-13, x9 = 0.412
y10 = 0 ~ 4.16e-17, err= 4.2e-17, x10 = 0.116
y11 = 0 ~ 2.22e-16, err= 2.2e-16, x11 = 0.565
y12 = 0 ~ 1.07e-13, err= 1.1e-13, x12 = 1.41
y13 = 1 ~ 1, err= 1.8e-13, x13 = 0.435
y14 = 0 ~ -7.43e-14, err=-7.4e-14, x14 = 6.25
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:
import numpy as np
from sympy import Symbol, Eq, solveset
for y, values in (
(1.690, (5.08354, 1663.125, 60, -45.622)),
(0.138, (5.24677, 1598.673, 60, -46.424)),
):
a = Symbol('a', real=True, infinite=False)
b = Symbol('b', real=True, infinite=False)
c = Symbol('c', real=True, infinite=False)
d = Symbol('d', real=True, infinite=False)
variables = a, b, c, d
eq = Eq(-np.log10(y), a - b/(c - d))
print('Do one of:')
for target, old_val in zip(variables, values):
seq = eq.subs({
s: v for s, v in zip(variables, values)
if s is not target
})
replacement, = solveset(seq, target)
print(f'Replace {target}={old_val} with {replacement:.1f}')
print()
Do one of:
Replace a=5.08354 with 15.5
Replace b=1663.125 with 561.0
Replace c=60 with 267.5
Replace d=-45.622 with -253.1
Do one of:
Replace a=5.24677 with 15.9
Replace b=1598.673 with 466.8
Replace c=60 with 318.0
Replace d=-46.424 with -304.4