(Sympy) Why do I get an unsupported operand error when trying to solve a specific equation block using Newton's Method despite other blocks working?

Question

I have a system of 58 equations (in which some equations are non-linear and/or multivariate) which I have separated into blocks (via Trajan's algorithm) for easier solving and faster computational time. A solution for the whole system definitely exists as I've tried it on a software (Engineering Equation Solver (EES)) and I'm getting the right answers.

Now, I'm trying to recreate/convert my original system on EES to Python via (mostly) Sympy. Right now, I'm near completion and all blocks can be solved EXCEPT one (it's maddening!!!).

The code in question relates to "Block 2" (B2) and Newtons Method where:

The equation/function matrix is "B2FunctionMatrix" and consists of 21 equations.
The variable, guess (initial), next guess and Jacobian matrices are labelled accordingly in the code.
The tolerance is 0.05
The solved function matrix is "B2fsolve"
The solved jacobian matrix is "B2jsolve"
The solved function matrix with respect to the jacobian is "B2delta_x0"
The residuals matrix is "B2guessVariance"

Throughout the code I have placed numerous "placeholders" to see where the code gets up to upon execution. The code is as follows:

B2FunctionMatrix = Matrix([[fxdot_m], [fa_2], [fMValve], [fPValve], [fq_3c], [fC_vm], [fq_m], [fq_1], [fC_vp], [fq_2], [fmix], [fp1q_m], [fh_L1], [fCV1q_m], [fp2q_m], [fh_L2], [fp3q_m], [fh_L3], [fCV1q_m], [fp4q_m], [fh_L4]]).subs(([rho, guessrho], [g, guessg], [m_m, guessm_m], [m_p, guessm_p], [a_1, guessa_1], [a_p, guessa_p], [k_spr, guessk_spr], [theta, guesstheta], [P_sp, guessP_sp], [C_vfo, guessC_vfo], [C_valpha, guessC_valpha], [C_vbeta, guessC_vbeta], [a_alpha, guessa_alpha], [a_beta, guessa_beta], [h_pump, guessh_pump], [h_atm, guessh_atm], [f, guessf], [L1, guessL1], [L2, guessL2], [L3, guessL3], [L4, guessL4], [D1, guessD1], [D2, guessD2], [D3, guessD3], [D4, guessD4], [A1, guessA1], [A2, guessA2], [A3, guessA3], [A4, guessA4], [P_spscale, guessP_spscale]))

B2VariableMatrix = Matrix([[x_m], [a_2], [q_3c], [h_in], [h_out], [h_c], [q_m], [x_p], [C_vm], [C_vp], [h_t], [q_1], [q_2], [h_L1], [h_L2], [h_L3], [h_L4], [v1], [v2], [v3], [v4]])


print('Block 2 (Multiple Unknown Equations)')
dof2 = abs(B2VariableMatrix.shape[0] - B2FunctionMatrix.shape[0])
if dof2 == 0:
    B2guessmatrix = B2VariableMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq3_c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4]))
    B2iter_n = 0
    B2tolerance = Matrix.ones(B2FunctionMatrix.shape[0], 1) * tolerance
    B2guessVariance = Matrix.ones(B2FunctionMatrix.shape[0], 1) * 10
    B2JacobianMatrix = B2FunctionMatrix.jacobian(B2VariableMatrix)
    
    
    while B2guessVariance[0] > B2tolerance[0] or B2guessVariance[1] > B2tolerance[1]:
        print('placeholder 1')
        B2fsolve = B2FunctionMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq3_c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4]))
        print('placeholder 2')
        B2jsolve = B2JacobianMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq3_c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4]))
        print('placeholder 3')
        B2delta_x0, fv = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
        print('placeholder 4')
        B2guessmatrix = B2VariableMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq3_c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4]))
        B2nextguessmatrix = B2delta_x0 + B2guessmatrix
        guessx_m = B2nextguessmatrix[0]
        guessa_2 = B2nextguessmatrix[1]
        guessq_3c = B2nextguessmatrix[2]
        guessh_in = B2nextguessmatrix[3]
        guessh_out = B2nextguessmatrix[4]
        guessh_c = B2nextguessmatrix[5]
        guessq_m = B2nextguessmatrix[6]
        guessx_p = B2nextguessmatrix[7]
        guessC_vm = B2nextguessmatrix[8]
        guessC_vp = B2nextguessmatrix[9]
        guessh_t = B2nextguessmatrix[10]
        guessq_1 = B2nextguessmatrix[11]
        guessq_2 = B2nextguessmatrix[12]
        guessh_L1 = B2nextguessmatrix[13]
        guessh_L2 = B2nextguessmatrix[14]
        guessh_L3 = B2nextguessmatrix[15]
        guessh_L4 = B2nextguessmatrix[16]
        guessv1 = B2nextguessmatrix[17]
        guessv2 = B2nextguessmatrix[18]
        guessv3 = B2nextguessmatrix[19]
        guessv4 = B2nextguessmatrix[20]        
        B2guessVariance = abs(B2nextguessmatrix - B2guessmatrix)
        B2guessmatrix = B2nextguessmatrix
        print('placeholder 5')
        B2iter_n += 1
        print(f'{B2iter_n} iterations have been completed so far. Moving onto the next iteration...')
        if B2iter_n >= max_iter:
            print('The maximum Number of iterations has been reached')
            break
        
    if B2guessVariance[0] <= B2tolerance[0] or B2guessVariance[1] <= B2tolerance[1]:
        print('The solution Matrix is')
        display(B2nextguessmatrix)
        print(f'Which was achieved after {B2iter_n} iterations with a tolerance of {tolerance}.')
        print(f'Displayed as integers, the solutions for variable {x_m} and {a_2} converge at {sp.Float(B2nextguessmatrix[0])} and {sp.Float(B2nextguessmatrix[1])} as floats respectively.')
    else:
        print('The Equation set has no solution or the initial guesses are too far from the solution.')
        
elif dof2 !=0:
    print(f'This system has {B2FunctionMatrix.shape[0]} equations and {B2VariableMatrix.shape[0]} variables which represents a d.o.f value of {dof2} which ≠ 0. Therefore, the system cannot be solved.')

I only get to "placeholder 3" in the code. To be clear, the output I get is:

Block 2 (Multiple Unknown Equations)
placeholder 1
placeholder 2
placeholder 3

Until I receive the following error which has me very confused:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\cache.py in wrapper(*args, **kwargs)
     71             try:
---> 72                 retval = cfunc(*args, **kwargs)
     73             except TypeError:

D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\operations.py in __new__(cls, evaluate, _sympify, *args)
     84 
---> 85         c_part, nc_part, order_symbols = cls.flatten(args)
     86         is_commutative = not nc_part

D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\mul.py in flatten(cls, seq)
    690         # order commutative part canonically
--> 691         _mulsort(c_part)
    692 

D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\mul.py in _mulsort(args)
     30     # in-place sorting of args
---> 31     args.sort(key=_args_sortkey)
     32 

D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\basic.py in compare(self, other)
    229             else:
--> 230                 c = (l > r) - (l < r)
    231             if c:

TypeError: unsupported operand type(s) for -: 'StrictGreaterThan' and 'StrictLessThan'

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
<ipython-input-22-9ad901ca6789> in <module>
     20         B2jsolve = B2JacobianMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq3_c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4]))
     21         print('placeholder 3')
---> 22         B2delta_x0, fv = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
     23         print('placeholder 4')
     24         B2guessmatrix = B2VariableMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq3_c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4]))

D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\matrices.py in gauss_jordan_solve(self, B, freevar)
   2158 
   2159     def gauss_jordan_solve(self, B, freevar=False):
-> 2160         return _gauss_jordan_solve(self, B, freevar=freevar)
   2161 
   2162     def pinv_solve(self, B, arbitrary_matrix=None):

D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\solvers.py in _gauss_jordan_solve(M, B, freevar)
    563 
    564     # solve by reduced row echelon form
--> 565     A, pivots = aug.rref(simplify=True)
    566     A, v      = A[:, :-B_cols], A[:, -B_cols:]
    567     pivots    = list(filter(lambda p: p < col, pivots))

D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\matrices.py in rref(self, iszerofunc, simplify, pivots, normalize_last)
    169             normalize_last=True):
    170         return _rref(self, iszerofunc=iszerofunc, simplify=simplify,
--> 171             pivots=pivots, normalize_last=normalize_last)
    172 
    173     echelon_form.__doc__ = _echelon_form.__doc__

D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\reductions.py in _rref(M, iszerofunc, simplify, pivots, normalize_last)
    304 
    305     mat, pivot_cols, _ = _row_reduce(M, iszerofunc, simpfunc,
--> 306             normalize_last, normalize=True, zero_above=True)
    307 
    308     if pivots:

D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\reductions.py in _row_reduce(M, iszerofunc, simpfunc, normalize_last, normalize, zero_above)
    127     mat, pivot_cols, swaps = _row_reduce_list(list(M), M.rows, M.cols, M.one,
    128             iszerofunc, simpfunc, normalize_last=normalize_last,
--> 129             normalize=normalize, zero_above=zero_above)
    130 
    131     return M._new(M.rows, M.cols, mat), pivot_cols, swaps

D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\reductions.py in _row_reduce_list(mat, rows, cols, one, iszerofunc, simpfunc, normalize_last, normalize, zero_above)
     67         pivot_offset, pivot_val, \
     68         assumed_nonzero, newly_determined = _find_reasonable_pivot(
---> 69                 get_col(piv_col)[piv_row:], iszerofunc, simpfunc)
     70 
     71         # _find_reasonable_pivot may have simplified some things

D:\ProgramData\Anaconda3\lib\site-packages\sympy\matrices\determinant.py in _find_reasonable_pivot(col, iszerofunc, simpfunc)
     80         if possible_zeros[i] is not None:
     81             continue
---> 82         simped = simpfunc(x)
     83         is_zero = iszerofunc(simped)
     84         if is_zero == True or is_zero == False:

D:\ProgramData\Anaconda3\lib\site-packages\sympy\simplify\simplify.py in simplify(expr, ratio, measure, rational, inverse, doit, **kwargs)
    658     if expr.has(Piecewise):
    659         # Fold into a single Piecewise
--> 660         expr = piecewise_fold(expr)
    661         # Apply doit, if doit=True
    662         expr = done(expr)

D:\ProgramData\Anaconda3\lib\site-packages\sympy\functions\elementary\piecewise.py in piecewise_fold(expr)
   1127             args = expr.args
   1128         # fold
-> 1129         folded = list(map(piecewise_fold, args))
   1130         for ec in cartes(*[
   1131                 (i.args if isinstance(i, Piecewise) else

D:\ProgramData\Anaconda3\lib\site-packages\sympy\functions\elementary\piecewise.py in piecewise_fold(expr)
   1132                  [(i, true)]) for i in folded]):
   1133             e, c = zip(*ec)
-> 1134             new_args.append((expr.func(*e), And(*c)))
   1135 
   1136     return Piecewise(*new_args)

D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\cache.py in wrapper(*args, **kwargs)
     72                 retval = cfunc(*args, **kwargs)
     73             except TypeError:
---> 74                 retval = func(*args, **kwargs)
     75             return retval
     76 

D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\operations.py in __new__(cls, evaluate, _sympify, *args)
     83             return args[0]
     84 
---> 85         c_part, nc_part, order_symbols = cls.flatten(args)
     86         is_commutative = not nc_part
     87         obj = cls._from_args(c_part + nc_part, is_commutative)

D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\mul.py in flatten(cls, seq)
    689 
    690         # order commutative part canonically
--> 691         _mulsort(c_part)
    692 
    693         # current code expects coeff to be always in slot-0

D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\mul.py in _mulsort(args)
     29 def _mulsort(args):
     30     # in-place sorting of args
---> 31     args.sort(key=_args_sortkey)
     32 
     33 

D:\ProgramData\Anaconda3\lib\site-packages\sympy\core\basic.py in compare(self, other)
    228                 c = l.compare(r)
    229             else:
--> 230                 c = (l > r) - (l < r)
    231             if c:
    232                 return c

TypeError: unsupported operand type(s) for -: 'StrictGreaterThan' and 'StrictLessThan'

The other 3 blocks solve completely fine. I was thinking that maybe it has something to do with the big size of the sparse Jacobian matrix in the block ("B2Jsolve" which is a sparse 21 x 21 square matrix)? Any ideas would be very appreciated! Thank you for reading :)

The output for B2fsolve is as follows: B2fsolve

The output for B2Jsolve is as follows: B2jsolve (This is a merged jpg file only with good resolution when you zoom in to make it clearer!)

EDIT (Adding source code): Here is the simplified version of the code which I've added ##headings to so that it's a bit clearer. The solving block is at the bottom of the code. I think that the problems may be to do with taking the Jacobian (B2jsolve) of some equations which use Sympy's sign function (original equations on lines 33, 34, 36, 47 & 55).

## 1.0 CODE SETTINGS, LIBRARIES & PHYSICAL CONSTANTS ##
from sympy.interactive import printing 
printing.init_printing(use_latex = True)
import numpy as np
import sympy as sp
import math
import functools
from sympy import Matrix 
from sympy.functions import sign

rho, g = sp.symbols('rho g')
cos = sp.cos
e = sp.exp
pi = np.pi
sqrt = sp.sqrt

## 2.0 PRV Matrix Constants ##
m_m, m_p, a_1, a_p, k_spr, theta = sp.symbols('m_m m_p a_1 a_p k_spr theta')

## 3.0 PRV FORCE BALANCE RELATIONSHIPS ##
xdot_m, a_2, q_3c, x_m, x_m0, t_0, t_f, t, P_sp = sp.symbols('xdot_m a_2 q_3c x_m x_m0 t_0 t_f t P_sp')
fa_2 = Matrix([a_2 - (1/(3700*(0.02732 - x_m)))])
fxdot_m = Matrix([0 - 3700*(0.02732 - x_m) * q_3c])

## 4.0 MAIN VALVE & PILOT VALVE FORCE BALANCES (m_m d^2m_m/dt^2) ##
h_in, h_out, h_c, q_m, x_p = sp.symbols('h_in h_out h_c q_m x_p')
fMValve = Matrix([rho * g * (h_in * a_1 + h_out * (a_2 - a_1) - h_c * a_2) - m_m * g * cos(theta) + (rho * ((q_m)**2)/a_1) - 0])
fPValve = Matrix([k_spr * (P_sp - x_p) - rho * g * h_out * a_p + m_p * g - 0])

## 5.0 FLOW EQUATIONS ##
C_vm, C_vfo, C_vp, C_vnvout, C_vnvin, h_t, q_1, q_2 = sp.symbols('C_vm C_vfo C_vp C_vnvout C_vnvin h_t q_1 q_2')
fC_vm = Matrix([0.02107 - 0.02962 * e(-51.1322 * x_m) + 0.0109 * e(-261 * x_m) - 0.00325 * e(-683.17 * x_m) + 0.0009 * e(-399.5 * x_m) - C_vm])
fq_m = Matrix([C_vm * (sqrt((abs(h_in - h_out)))) * sign(h_in - h_out) -  q_m])
fq_1 = Matrix([C_vfo * (sqrt(abs(h_in - h_t))) * sign(h_in - h_t) - q_1])
fC_vp = Matrix([0.0000753 * (1 - e(-1135 * x_p)) - C_vp])
fq_2 = Matrix([C_vp * (sqrt(abs(h_t - h_out))) * sign(h_t - h_out) - q_2])
fmix = Matrix([q_1 + q_3c - q_2])
fq_3c = Matrix([h_c - h_t])

## 6.0 NETWORK EQUATIONS FROM PUMP TO ATMOSPHERIC OUTLET ##
C_valpha, C_vbeta, a_alpha, a_beta, h_pump, h_L1, h_L2, h_L3, h_L4, h_atm, f, L1, L2, L3, L4, D1, D2, D3, D4, A1, A2, A3, A4, v1, v2, v3, v4 = sp.symbols('C_valpha C_vbeta a_alpha a_beta h_pump h_L1 h_L2 h_L3 h_L4 h_atm f L1 L2 L3 L4 D1 D2 D3 D4 A1 A2 A3 A4 v1 v2 v3 v4')

#Pipe1
fp1q_m = Matrix([q_m - (A1 * v1)])
fh_L1 = Matrix([h_L1 - (f*(L1/D1)*((v1**2)/(2*g)))])
#AlphaValve
fCV1q_m = Matrix([C_valpha * a_alpha * (sqrt(abs(h_pump - h_L1 - (h_in + h_L2)))) * sign(h_pump - h_L1 - (h_in + h_L2)) - q_m])
#Pipe2
fp2q_m = Matrix([q_m - (A2 * v2)])
fh_L2 = Matrix([h_L2 - (f*(L2/D2)*((v2**2)/(2*g)))])
#Pipe3
fp3q_m = Matrix([q_m - A3 * v3])
fh_L3 = Matrix([h_L3 - (f*(L3/D3)*((v3**2)/(2*g)))])
#BetaValve
fCV2q_m = Matrix([C_vbeta * a_beta * (sqrt(abs(h_out - h_L3 - (h_atm + h_L4)))) * sign(h_out - h_L3 - (h_atm + h_L4)) - q_m])
#Pipe4
fp4q_m = Matrix([q_m - (A4 * v4)])
fh_L4 = Matrix([h_L4 - (f*(L4/D4)*((v4**2)/(2*g)))])

q_mscale, q_1scale, q_2scale, q_3cscale, x_mscale, x_pscale, P_spscale, a_2scale = sp.symbols('q_mscale q_1scale q_2scale q_3cscale x_mscale x_pscale P_spscale a_2scale')

## 7.1 INITIAL GUESSES & SOLVER CONDITIONS ##
guessrho, guessg, guessm_m, guessm_p, guessa_1, guessa_p, guessk_spr, guesstheta, guessP_sp, guessx_m, guessa_2, guessq_3c, guessh_in, guessh_out, guessh_c, guessq_m, guessx_p, guessC_vm, guessC_vfo, guessC_vp, guessh_t, guessq_1, guessq_2, guessC_valpha, guessC_vbeta, guessa_alpha, guessa_beta, guessh_pump, guessh_L1, guessh_L2, guessh_L3, guessh_L4, guessh_atm, guessf, guessL1, guessL2, guessL3, guessL4, guessD1, guessD2, guessD3, guessD4, guessA1, guessA2, guessA3, guessA4, guessv1, guessv2, guessv3, guessv4, guessq_mscale, guessq_1scale, guessq_2scale, guessq_3cscale, guessx_mscale, guessx_pscale, guessP_spscale, guessa_2scale = sp.symbols('guessrho guessg guessm_m guessm_p guessa_1 guessa_p guessk_spr guesstheta guessP_sp guessx_m guessa_2 guessq_3c guessh_in guessh_out guessh_c guessq_m guessx_p guessC_vm guessC_vfo guessC_vp guessh_t guessq_1 guessq_2 guessC_valpha guessC_vbeta guessa_alpha guessa_beta guessh_pump guessh_L1 guessh_L2 guessh_L3 guessh_L4 guessh_atm guessf guessL1 guessL2 guessL3 guessL4 guessD1 guessD2 guessD3 guessD4 guessA1 guessA2 guessA3 guessA4 guessv1 guessv2 guessv3 guessv4 guessq_mscale guessq_1scale guessq_2scale guessq_3cscale guessx_mscale guessx_pscale guessP_spscale guessa_2scale')

guessrho = 997 # physical constant
guessg = 9.81 # physical constant
guessm_m = 8 # physical constant
guessm_p = 0.1 # physical constant
guessa_1 = 0.0078 # physical constant
guessa_p = 0.00196 # physical constant
guessk_spr = 70000 # physical constant
guesstheta = 0 #0.000000000000000000000001 # physical constant 31/5/21 CHANGED SO IT'S NON-ZERO
guessP_sp = 0.00700000000 
guessx_m = 0.013647064
guessa_2 = 0.01977
guessq_3c = -2.578E-24 #0.000000000000000000000001 # physical constant 31/5/21 CHANGED SO IT'S NON-ZERO
guessh_in = 41.38
guessh_out = 23.465726
guessh_c = 30.65
guessq_m = 0.02811
guessx_p = 0.0005878
guessC_vm = 0.006642
guessC_vfo = 0.00003
guessC_vp = 0.00003666
guessh_t = 30.65
guessq_1 = 0.00009826
guessq_2 = 0.00009826
guessC_valpha = 1
guessC_vbeta = 0.6
guessa_alpha = 0.01
guessa_beta = 0.01
guessh_pump = 50
guessh_L1 = 0.1959
guessh_L2 = 0.5224
guessh_L3 = 0.6529
guessh_L4 = 0.8619
guessh_atm = 0 #0.000000000000000000000001 # 31/5/21 CHANGED SO IT'S NON-ZERO
guessf = 0.02
guessL1 = 1.5
guessL2 = 4
guessL3 = 5
guessL4 = 6.6
guessD1 = 0.1
guessD2 = 0.1
guessD3 = 0.1
guessD4 = 0.1
guessA1 = 0.007854
guessA2 = 0.007854
guessA3 = 0.007854
guessA4 = 0.007854
guessv1 = 3.579
guessv2 = 3.579
guessv3 = 3.579
guessv4 = 3.579
guessq_mscale = 101.2 #q_m * 3600
guessq_1scale = 0.3537 #q_1 * 3600
guessq_2scale = 0.3537 #q_2 * 3600
guessq_3cscale = -9.281E-21 #q_3c * 36000
guessx_mscale = 13.6470637 #x_m * 1000
guessx_pscale = 0.5878 #x_p * 1000
guessP_spscale = 7 #P_sp * 1000
guessa_2scale = 197.7 #a_2 * (100^2)

#Broad Loop Solver Conditions
tolerance = 0.05
max_iter = 10

#!-!#! 7.2.2 BLOCK 2 (MULTIPLE UNKNOWN EQUATIONS) (THE PROBLEMS ARE WITH THIS BLOCK!!!!) #!-!#!-!#!-

print('*PROBLEM* SECTION 7.2.2 BLOCK 2 (MULTIPLE UNKNOWN EQUATIONS)')
B2FunctionMatrix = Matrix([[fxdot_m], [fa_2], [fMValve], [fPValve], [fq_3c], [fC_vm], [fq_m], [fq_1], [fC_vp], [fq_2], [fmix], [fp1q_m], [fh_L1], [fCV1q_m], [fp2q_m], [fh_L2], [fp3q_m], [fh_L3], [fCV1q_m], [fp4q_m], [fh_L4]]).subs(([rho, guessrho], [g, guessg], [m_m, guessm_m], [m_p, guessm_p], [a_1, guessa_1], [a_p, guessa_p], [k_spr, guessk_spr], [theta, guesstheta], [P_sp, guessP_sp], [C_vfo, guessC_vfo], [C_valpha, guessC_valpha], [C_vbeta, guessC_vbeta], [a_alpha, guessa_alpha], [a_beta, guessa_beta], [h_pump, guessh_pump], [h_atm, guessh_atm], [f, guessf], [L1, guessL1], [L2, guessL2], [L3, guessL3], [L4, guessL4], [D1, guessD1], [D2, guessD2], [D3, guessD3], [D4, guessD4], [A1, guessA1], [A2, guessA2], [A3, guessA3], [A4, guessA4], [P_spscale, guessP_spscale]))

B2VariableMatrix = Matrix([[x_m], [a_2], [q_3c], [h_in], [h_out], [h_c], [q_m], [x_p], [C_vm], [C_vp], [h_t], [q_1], [q_2], [h_L1], [h_L2], [h_L3], [h_L4], [v1], [v2], [v3], [v4]])

print('Block 2 (Multiple Unknown Equations)')
dof2 = abs(B2VariableMatrix.shape[0] - B2FunctionMatrix.shape[0])
if dof2 == 0:
    B2guessmatrix = B2VariableMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq_3c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4]))
    #display(B2guessmatrix)
    #display(B2FunctionMatrix)
    B2iter_n = 0
    B2tolerance = Matrix.ones(B2FunctionMatrix.shape[0], 1) * tolerance
    B2guessVariance = Matrix.ones(B2FunctionMatrix.shape[0], 1) * 10
    B2JacobianMatrix = B2FunctionMatrix.jacobian(B2VariableMatrix)
    #display(B2JacobianMatrix)
    
    while (B2guessVariance[0,0] >= B2tolerance[0,0]) and (B2guessVariance[1,0] >= B2tolerance[1,0]) and (B2guessVariance[2,0] >= B2tolerance[2,0]) and (B2guessVariance[3,0] >= B2tolerance[3,0]) and (B2guessVariance[4,0] >= B2tolerance[4,0]) and (B2guessVariance[5,0] >= B2tolerance[5,0]) and (B2guessVariance[6,0] >= B2tolerance[6,0]) and (B2guessVariance[7,0] >= B2tolerance[7,0]) and (B2guessVariance[8,0] >= B2tolerance[8,0]) and (B2guessVariance[9,0] >= B2tolerance[9,0]) and (B2guessVariance[10,0] >= B2tolerance[10,0]) and (B2guessVariance[11,0] >= B2tolerance[11,0]) and (B2guessVariance[12,0] >= B2tolerance[12,0]) and (B2guessVariance[13,0] >= B2tolerance[13,0]) and (B2guessVariance[14,0] >= B2tolerance[14,0]) and (B2guessVariance[15,0] >= B2tolerance[15,0]) and (B2guessVariance[16,0] >= B2tolerance[16,0]) and (B2guessVariance[17,0] >= B2tolerance[17,0]) and (B2guessVariance[18,0] >= B2tolerance[18,0]) and (B2guessVariance[19,0] >= B2tolerance[19,0]) and (B2guessVariance[20,0] >= B2tolerance[20,0]):
        print('placeholder 1 (displays B2fsolve, the function matrix with guess values substituted)')
        B2fsolve = B2FunctionMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq_3c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4]))
        print(repr(B2fsolve))
        #display(B2fsolve)
        print('placeholder 2 (displays B2jsolve, the jacobian matrix with guess values substituted)')
        B2jsolve = B2JacobianMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq_3c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4]))
        print(repr(B2jsolve))
        #display(B2jsolve)
        print('placeholder 3')
        B2delta_x0, fv = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
        #The code does not run past this point, please help!
        print(repr(B2delta_x0))
        #display(B2delta_x0)
        print('placeholder 4')
        B2guessmatrix = B2VariableMatrix.subs(([x_m, guessx_m], [a_2, guessa_2], [q_3c, guessq_3c], [h_in, guessh_in], [h_out, guessh_out], [h_c, guessh_c], [q_m, guessq_m], [x_p, guessx_p], [C_vm, guessC_vm], [C_vp, guessC_vp], [h_t, guessh_t], [q_1, guessq_1], [q_2, guessq_2], [h_L1, guessh_L1], [h_L2, guessh_L2], [h_L3, guessh_L3], [h_L4, guessh_L4], [v1, guessv1], [v2, guessv2], [v3, guessv3], [v4, guessv4]))
        B2nextguessmatrix = B2delta_x0 + B2guessmatrix
        guessx_m = B2nextguessmatrix[0]
        guessa_2 = B2nextguessmatrix[1]
        guessq_3c = B2nextguessmatrix[2]
        guessh_in = B2nextguessmatrix[3]
        guessh_out = B2nextguessmatrix[4]
        guessh_c = B2nextguessmatrix[5]
        guessq_m = B2nextguessmatrix[6]
        guessx_p = B2nextguessmatrix[7]
        guessC_vm = B2nextguessmatrix[8]
        guessC_vp = B2nextguessmatrix[9]
        guessh_t = B2nextguessmatrix[10]
        guessq_1 = B2nextguessmatrix[11]
        guessq_2 = B2nextguessmatrix[12]
        guessh_L1 = B2nextguessmatrix[13]
        guessh_L2 = B2nextguessmatrix[14]
        guessh_L3 = B2nextguessmatrix[15]
        guessh_L4 = B2nextguessmatrix[16]
        guessv1 = B2nextguessmatrix[17]
        guessv2 = B2nextguessmatrix[18]
        guessv3 = B2nextguessmatrix[19]
        guessv4 = B2nextguessmatrix[20]        
        B2guessVariance = abs(B2nextguessmatrix - B2guessmatrix)
        B2guessmatrix = B2nextguessmatrix
        print('placeholder 5')
        B2iter_n += 1
        print(f'{B2iter_n} iterations have been completed so far. Moving onto the next iteration...')
        if B2iter_n >= max_iter:
            print('The maximum Number of iterations has been reached')
            break
        
    if (B2guessVariance[0,0] < B2tolerance[0,0]) and (B2guessVariance[1,0] < B2tolerance[1,0]) and (B2guessVariance[2,0] < B2tolerance[2,0]) and (B2guessVariance[3,0] < B2tolerance[3,0]) and (B2guessVariance[4,0] < B2tolerance[4,0]) and (B2guessVariance[5,0] < B2tolerance[5,0]) and (B2guessVariance[6,0] < B2tolerance[6,0]) and (B2guessVariance[7,0] < B2tolerance[7,0]) and (B2guessVariance[8,0] < B2tolerance[8,0]) and (B2guessVariance[9,0] < B2tolerance[9,0]) and (B2guessVariance[10,0] < B2tolerance[10,0]) and (B2guessVariance[11,0] < B2tolerance[11,0]) and (B2guessVariance[12,0] < B2tolerance[12,0]) and (B2guessVariance[13,0] < B2tolerance[13,0]) and (B2guessVariance[14,0] < B2tolerance[14,0]) and (B2guessVariance[15,0] < B2tolerance[15,0]) and (B2guessVariance[16,0] < B2tolerance[16,0]) and (B2guessVariance[17,0] < B2tolerance[17,0]) and (B2guessVariance[18,0] < B2tolerance[18,0]) and (B2guessVariance[19,0] < B2tolerance[19,0]) and (B2guessVariance[20,0] < B2tolerance[20,0]):
        print('The solution Matrix is')
        print(repr(B2nextguessmatrix))
        #display(B2nextguessmatrix)
        print(f'Which was achieved after {B2iter_n} iterations with a tolerance of {tolerance}.')
        print(f'Displayed as integers, the solutions for variable {x_m} and {a_2} converge at {sp.Float(B2nextguessmatrix[0])} and {sp.Float(B2nextguessmatrix[1])} as floats respectively.')
    else:
        print('The Equation set has no solution or the initial guesses are too far from the solution.')
        
elif dof2 !=0:
    print(f'This system has {B2FunctionMatrix.shape[0]} equations and {B2VariableMatrix.shape[0]} variables which represents a d.o.f value of {dof2} which ≠ 0. Therefore, the system cannot be solved.')

You don't need to add placeholders if you look at the stacktrace output. It shows you which line of your code the error happens on. The error comes from having a Boolean in place of an Expr but I don't know whether that's a bug in your code or in SymPy itself. The code you have shown is incomplete so I can't tell what's actually in the expression. Rather than completing the code you should just show the repr of `B2jsolve` and `B2fsolve` at the point where the exception is raised (i.e. make a minimal complete code example). — Oscar Benjamin, Jul 07 '21 at 12:40
Thanks for your comment @OscarBenjamin, I have added the photos in the main question now at the bottom. Apologies I didn't do it before. Very keen on your insight! :) — Hendrix13, Jul 13 '21 at 03:03
Images of the output are not what I meant by a code example. My guess is that this is a bug in sympy but you need to give runnable code for anyone else to check that. — Oscar Benjamin, Jul 13 '21 at 12:19
@OscarBenjamin I tried adding my complete runnable code in the body text but it made the whole question exceed the body count 30,000 character limit. Would you recommend me opening up a "Part 2" of this question so that I can post my runnable code? I've shortened it as much as I can. :/ I'm also thinking that the issues may be with the use of Sympy's signum functions as I use them in five of my 58 equations. Do you have any experience using SymPy's signum functions ("sign" from sympy.functions)? — Hendrix13, Jul 22 '21 at 09:09
Don't post all of your code. Extract only the part that is necessary. As I said you can use repr to get the code needed to create the expressions that are being passed into the function call that fails. http://www.sscce.org/ — Oscar Benjamin, Jul 22 '21 at 11:37
@OscarBenjamin I've just added the simplified code example to the main question. You should be able to copy/paste it and run to see the outputs I refer to in other parts of the quesion. Thank you so much for your help and patience so far. Very keen on your reply! — Hendrix13, Jul 23 '21 at 07:04
Thanks. This is a bug: https://github.com/sympy/sympy/issues/21773 — Oscar Benjamin, Jul 23 '21 at 23:51
@OscarBenjamin wicked, thank you so much for that! Hopefully there'll be a solution soon! Do you think maybe trying to embed a signum function made from scratch (i.e. a for loop) for those 5 equations might work? I'll keep an eye on github, thanks again for this :) — Hendrix13, Jul 25 '21 at 11:08
You might find that just calling `.doit()` in the right place fixes this. The `Subs` objects shown in the sympy issue are equivalent to zero. — Oscar Benjamin, Jul 25 '21 at 20:36
@OscarBenjamin could you expand on this please? Do you mean putting the equations that use the `sign` function in the loop before the Jacobian is calculated with a `.doit()` so that they're evaluated before being put into the Jacobian Matrix? I'm not sure what your last statement is referring to as most of the initial guesses (like the ones used in the `sign` part of some of the equations) are non-zero. The `Subs` seems to work fine for all variables from what I've seen so far? — Hendrix13, Jul 26 '21 at 02:27
You could call `doit` to evaluate the `Subs` just before the call that fails. I haven't traced through to find where the `Subs` are coming from so there might be a better place. — Oscar Benjamin, Jul 26 '21 at 15:35
@OscarBenjamin `.doit()` was added to all subs matrices and then to all equations and loop operations and the output was no different than before in all examples. However, the error code was slightly different referring to the `doit()`. Unfortunately, due to the body word limit I can't display this new error message in the main body text. One thing I also found when looking at the `display` (not `repr`) of `B2JacobianMatrix` was that there were some squared `sign` terms which doesn't seem right to me. Looking at `B2jsolve` there were also some substitutions notated but not actually substituted — Hendrix13, Aug 03 '21 at 03:26
You'll find this easier to figure out I think if you distil your problem down to the minimum that is needed to set up the call that fails. — Oscar Benjamin, Aug 04 '21 at 08:47
@OscarBenjamin That sounds like a good approach. Because I can't add anymore body text I might make a new question related to this but with minimal code. I will make a reference to this question in the new post too. Cheers! — Hendrix13, Aug 04 '21 at 09:10

(Sympy) Why do I get an unsupported operand error when trying to solve a specific equation block using Newton's Method despite other blocks working?

0 Answers0

Linked