optimizationscipyconstraintsminimize

Nonlinear equation system with constraints


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)

Solution

  • 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