issue with scipy.optimize.root returing non-zero errors

Question

The code below is a simplified version of a set of non-linear equations I am trying to solve using scipy.optimize.root. The if statements associated with variables e_eqn3 and e_eqn4 in function system_err try to prevent the solver from heading off into non-physical solution space. Everything seems to work fine for x12t >= 0.26 in the main program, but at lower values, I get back my initial guess, as the following output shows:

x12t = 0.25
           q12    p12
Initial: 19.00  1695.00
Final:   19.00  1695.00
Exited successfully: True

x12t = 0.26
           q12    p12
Initial: 19.00  1695.00
Final:   21.62  2187.55
Exited successfully: True

I checked this using EXCEL and I get a sensible solution for x12t = 0.25. I also checked the values being returned by system_err for x12t < 0.25 and an example of the final values before successful convergence is [1000000.0, 0.0, 1.378, 0.0, 0.0], which I thought meant that the overall error was still too high. Would appreciate any comments about either what I am doing wrong or a possible fix. Thanks in advance!

import math
import scipy.optimize
import numpy as np

# Calculate the system error

def system_err(x, pres, psgp, ptgp, x12t, x11t, x11s):
    
    # Unpack the unknowns
    
    p11 = x[0]
    p12 = x[1]
    q11 = x[2]
    q12 = x[3]
    r   = x[4]
    
    # Define constants
    
    a = 1.5152
    b = 1.7576
    c = 1.3433
    
    beta = 0.0129
    
    C11 = 0.010
    C12 = 0.018
    n11 = 0.595
    n12 = 0.533
    
    # Define the error associated with the constituent equations
    
    if p12 > pres:
        e_eqn3 = 1.0E06
    else:
        e_eqn3 = q12 - C12 * ((pres ** 2 - p12 ** 2) ** n12)
        
    if p11 > pres:
        e_eqn4 = 1.0E06
    else:
        e_eqn4 = q11 - C11 * ((pres ** 2 - p11 ** 2) ** n11)
        
    e_eqn5 = q12 - c * x12t * p12 * math.sqrt(beta * ( (ptgp / p12) ** a - (ptgp / p12) ** b) )
    e_eqn6 = r * q11 - c * x11t * p11 * math.sqrt(beta * ( (ptgp / p11) ** a - (ptgp / p11) ** b) )
    e_eqn7 = (1.0 - r) * q11 - c * x11s * p11 * math.sqrt(beta * ( (psgp / p11) ** a - (psgp / p11) ** b) )
    
    # Return the errors
    
    err = [e_eqn3, e_eqn4, e_eqn5, e_eqn6, e_eqn7]
    
    return err

# MAIN program

if __name__ == "__main__":
    
    # Initialise
    
    psgp = 1000.0           # plant 1 pressure (psia)
    ptgp = 1200.0           # plant 2 pressure (psia)

    pres = 2190.0           # inlet pressure (psia)
    
    x11s = 0.0              # FCV 1 setting
    x11t = 0.0              # FCV 2 setting
    
    # Vary FCV 3 setting
    
    for x12t in np.arange(0.20, 0.301, 0.01):
        
        # Initial guess for solution
        
        p11 = pres
        p12 = (ptgp + pres) / 2.0
        q11 = 0.0
        q12 = 50.0
        r = 0.0
        
        x0 = [p11, p12, q11, q12, r]
         
        print("\nx12t = %4.2f" % x12t)
        print("           q12    p12")
        print("Initial: %5.2f  %7.2f" % (q12, p12) )
        
        # Solve system
        
        sol = scipy.optimize.root(system_err, x0, args=(pres, psgp, ptgp, x12t, x11t, x11s))
        
        # Unpack the solution
        
        p11 = sol.x[0]
        p12 = sol.x[1]
        q11 = sol.x[2]
        q12 = sol.x[3]
        r   = sol.x[4]
        
        print("Final:   %5.2f  %7.2f" % (q12, p12))
        print("Exited successfully: %s" % sol.success)

Answer 1

Thought a bit harder about the physics of the problem and when eg p12 > pres then the flow direction reverses. This leads to raising a negative number to a fractional power which the original fix tried to prevent. Changed the code to represent the flow reversal explicitly

    if p12 > pres:
        e_eqn3 = q12 + C12 * ((p12 ** 2 - pres ** 2) ** n12)
    else:
        e_eqn3 = q12 - C12 * ((pres ** 2 - p12 ** 2) ** n12)

and the issue seems to have gone away.