python - （Sympy）为什么在尝试使用牛顿法求解特定方程块时，尽管其他块工作，但我得到一个不受支持的操作数错误？
问题描述
我有一个由 58 个方程组成的系统（其中一些方程是非线性的和/或多元的），我将它们分成块（通过 Trajan 算法）以便于求解和加快计算时间。整个系统的解决方案肯定存在，因为我已经在软件（工程方程求解器（EES））上尝试过，我得到了正确的答案。
现在，我正在尝试通过（主要是）Sympy 将我在 EES 上的原始系统重新创建/转换为 Python。现在，我快完成了，所有的块都可以解决，除了一个（这太疯狂了！！！）。
有问题的代码与“Block 2”（B2）和牛顿法有关，其中：
- 方程/函数矩阵是“B2FunctionMatrix”，由 21 个方程组成。
- 变量、猜测（初始）、下一个猜测和雅可比矩阵在代码中进行了相应标记。
- 公差为 0.05
- 求解的函数矩阵为“B2fsolve”
- 求解的雅可比矩阵是“B2jsolve”
- 关于雅可比的求解函数矩阵是“B2delta_x0”
- 残差矩阵是“B2guessVariance”
在整个代码中，我放置了许多“占位符”，以查看代码在执行时的位置。代码如下：
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.')
我只能在代码中找到“占位符 3”。需要明确的是，我得到的输出是：
Block 2 (Multiple Unknown Equations)
placeholder 1
placeholder 2
placeholder 3
直到我收到以下让我非常困惑的错误：
---------------------------------------------------------------------------
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'
其他 3 块解决完全没问题。我在想这可能与块中稀疏雅可比矩阵的大尺寸有关（“B2Jsolve”是一个稀疏的 21 x 21 方阵）？任何想法将不胜感激！感谢您阅读:)
B2fsolve 的输出如下： B2fsolve
B2Jsolve 的输出如下： B2jsolve （这是一个合并的 jpg 文件，只有当您放大以使其更清晰时才具有良好的分辨率！）
编辑（添加源代码）：这是我添加 ##headings 的代码的简化版本，以便更清晰一些。求解块位于代码的底部。我认为问题可能与采用某些使用 Sympy 符号函数的方程的雅可比 (B2jsolve) 有关（第 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.')
解决方案
推荐阅读
- apache-spark - 禁用或强制 Spark 哈希分区
- python - GUI、无限循环进程和线程
- python - Python将列表列表转换为具有公共列和索引的数据框
- javascript - 当过滤一个
- 子弹消失。如何让他们留下来？
- python - Discord.py 加入/连接到语音通道命令不起作用。（在齿轮中）
- html - 具有固定定位的 Z 索引在 CodePen 中不起作用
- sql - Postgres 按子查询/多个子查询排序
- amazon-web-services - AWS Amplify GraphQL - 一对多连接在查询时返回空列表
- r - 更改 tibble 中变量的值
- filter - 如何将某些项目值转换为新项目值？