使用 Sympy 求解方程组时,如何删除输出中的自动 tau0 变量?
How do I remove the automated tau0 variables in my output when solving a system of equations using Sympy?
我制作了一个求解器,它适用于我尝试过的其他几个方程组。但是,我当前的系统有 21 个方程(一些非线性方程和多变量方程),当放入同一个求解器时似乎不起作用。基于牛顿法的矩阵求解器根据 while 循环运行,该循环将猜测的输出残差(称为 B2guessVariance
)与指定的容差进行比较。
出于某种原因,并且仅在本示例中,我得到了一个奇怪的输出。在一次循环迭代后,我的 B2guessVariance
矩阵有一个嵌套在矩阵中的 tau0 变量。我的代码中没有使用 tau 符号的变量。我怀疑它可能与我在 Sympy 中使用的 gauss_jordan_solve
运算符有关?
关于系统的一些注意事项,以防系统在任何评论中被抢购一空:
- 我的 21 方程系统的解决方案肯定存在,因为我已经在软件(工程方程求解器 (EES))上尝试过,并且我得到了正确的答案。
- 求解的函数矩阵为
B2fsolve
- 求解的雅可比矩阵为
B2jsolve
- 关于雅可比的求解函数矩阵是
B2delta_x0
- 残差矩阵为
B2guessVariance
- 无论输入的方程是线性的、非线性的、and/or 多变量等,我制作的求解器以前都有效
在整个代码中,我放置了许多“占位符”以查看代码在执行时到达的位置。我在我的代码中使用 sympy displays
但知道 print(repr())
在大多数 programs/terminals 中更可取;因此,displays
被注释掉并使用 print(repr())
。缩短后的代码如下:
## 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)))]).doit()
fxdot_m = Matrix([0 - 3700*(0.02732 - x_m) * q_3c]).doit()
## 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]).doit()
fPValve = Matrix([k_spr * (P_sp - x_p) - rho * g * h_out * a_p + m_p * g - 0]).doit()
## 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]).doit()
fq_m = Matrix([C_vm * (sqrt(((h_in - h_out)))) - q_m]).doit()
fq_1 = Matrix([C_vfo * (sqrt((h_in - h_t))) - q_1]).doit()
fC_vp = Matrix([0.0000753 * (1 - e(-1135 * x_p)) - C_vp]).doit()
fq_2 = Matrix([C_vp * (sqrt((h_t - h_out))) - q_2]).doit()
fmix = Matrix([q_1 + q_3c - q_2]).doit()
fq_3c = Matrix([h_c - h_t]).doit()
## 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)]).doit()
fh_L1 = Matrix([h_L1 - (f*(L1/D1)*((v1**2)/(2*g)))]).doit()
#AlphaValve
fCV1q_m = Matrix([C_valpha * a_alpha * (sqrt((h_pump - h_L1 - (h_in + h_L2)))) - q_m]).doit()
#Pipe2
fp2q_m = Matrix([q_m - (A2 * v2)]).doit()
fh_L2 = Matrix([h_L2 - (f*(L2/D2)*((v2**2)/(2*g)))]).doit()
#Pipe3
fp3q_m = Matrix([q_m - A3 * v3]).doit()
fh_L3 = Matrix([h_L3 - (f*(L3/D3)*((v3**2)/(2*g)))]).doit()
#BetaValve
fCV2q_m = Matrix([C_vbeta * a_beta * (sqrt((h_out - h_L3 - (h_atm + h_L4)))) - q_m]).doit()
#Pipe4
fp4q_m = Matrix([q_m - (A4 * v4)]).doit()
fh_L4 = Matrix([h_L4 - (f*(L4/D4)*((v4**2)/(2*g)))]).doit()
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])).doit()
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]]).doit()
#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])).doit()
#display(B2guessmatrix)
#display(B2FunctionMatrix)
B2iter_n = 0
B2tolerance = (Matrix.ones(B2FunctionMatrix.shape[0], 1) * tolerance).doit()
B2guessVariance = (Matrix.ones(B2FunctionMatrix.shape[0], 1) * 10).doit()
B2JacobianMatrix = (B2FunctionMatrix.jacobian(B2VariableMatrix)).doit()
#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])).doit()
#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])).doit()
#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])).doit()
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')
print('These are the two matrices which are compared in the while loop')
#display(B2guessVariance)
#display(B2tolerance)
print(repr(B2guessVariance))
print(repr(B2tolerance))
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.')
显示while循环中用于迭代的关系矩阵的输出如下:
*PROBLEM* SECTION 7.2.2 BLOCK 2 (MULTIPLE UNKNOWN EQUATIONS)
These are the two matrices which are compared in the while loop
Matrix([
[ Abs(0.00392183668209074*tau0 - 8.02180713197279e-7)],
[ Abs(0.00566975402880878*tau0 - 4.35392954981073e-6)],
[ Abs(7.3945310403193e-25*tau0 - 2.5781512492912e-24)],
[ Abs(4.81682791725627*tau0 - 0.000426049277400864)],
[ Abs(0.617506604286677*tau0 - 0.000126658199850027)],
[ Abs(4.04275751354693*tau0 - 0.00184684018329051)],
[ Abs(0.007854*tau0 - 5.33999999999951e-7)],
[Abs(0.000169107863923268*tau0 - 4.93212155862604e-8)],
[ Abs(0.00263411194406529*tau0 - 7.57176496023193e-7)],
[ Abs(7.41675230154456e-6*tau0 - 3.67497300027552e-9)],
[ Abs(4.04275751354693*tau0 - 0.00184684018329051)],
[ Abs(3.54463776960747e-6*tau0 - 3.53001248327642e-9)],
[ Abs(3.54463776960747e-6*tau0 - 3.53001248327642e-9)],
[ Abs(0.109449541284404*tau0 - 4.00458715596186e-5)],
[ Abs(0.291865443425076*tau0 - 0.000106788990825613)],
[ Abs(0.364831804281346*tau0 - 3.34862385320545e-5)],
[ Abs(0.481577981651376*tau0 - 0.000116201834862384)],
[ 1.0*Abs(tau0)],
[ 1.0*Abs(tau0)],
[ 1.0*Abs(tau0)],
[ Abs(tau0)]])
Matrix([
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05]])
1 iterations have been completed so far. Moving onto the next iteration...
这给出了错误:
TypeError Traceback (most recent call last)
<ipython-input-1-e583c0741f11> in <module>
145 #display(B2JacobianMatrix)
146
--> 147 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]):
148 #print('placeholder 1 (displays B2fsolve, the function matrix with guess values substituted)')
149 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])).doit()
C:\ProgramData\Anaconda3\lib\site-packages\sympy\core\relational.py in __bool__(self)
396
397 def __bool__(self):
--> 398 raise TypeError("cannot determine truth value of Relational")
399
400 def _eval_as_set(self):
TypeError: cannot determine the truth value of Relational
显然,由于tau0符号在输出的B2guessVariance
矩阵中,因此while循环继续迭代时,不能与容差矩阵相关联。有没有一种方法可以删除 tau0 变量,以便求解器可以正常执行并且我可以解决这个问题?
After some research,这似乎可以通过用 params
return 行预先指定 gauss_jordan_solve
然后用零替换 taus 来解决。
因此,通过更改以下内容解决了这个问题:
B2delta_x0, fv = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
至:
B2delta_x0, params = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
B2delta_x0 = B2delta_x0.xreplace({ tau:0 for tau in params })
不幸的是,我并不真正理解为什么 这会起作用,这令人沮丧。但它目前正在运行,并为我提供了很好的答案,这太棒了!
我制作了一个求解器,它适用于我尝试过的其他几个方程组。但是,我当前的系统有 21 个方程(一些非线性方程和多变量方程),当放入同一个求解器时似乎不起作用。基于牛顿法的矩阵求解器根据 while 循环运行,该循环将猜测的输出残差(称为 B2guessVariance
)与指定的容差进行比较。
出于某种原因,并且仅在本示例中,我得到了一个奇怪的输出。在一次循环迭代后,我的 B2guessVariance
矩阵有一个嵌套在矩阵中的 tau0 变量。我的代码中没有使用 tau 符号的变量。我怀疑它可能与我在 Sympy 中使用的 gauss_jordan_solve
运算符有关?
关于系统的一些注意事项,以防系统在任何评论中被抢购一空:
- 我的 21 方程系统的解决方案肯定存在,因为我已经在软件(工程方程求解器 (EES))上尝试过,并且我得到了正确的答案。
- 求解的函数矩阵为
B2fsolve
- 求解的雅可比矩阵为
B2jsolve
- 关于雅可比的求解函数矩阵是
B2delta_x0
- 残差矩阵为
B2guessVariance
- 无论输入的方程是线性的、非线性的、and/or 多变量等,我制作的求解器以前都有效
在整个代码中,我放置了许多“占位符”以查看代码在执行时到达的位置。我在我的代码中使用 sympy displays
但知道 print(repr())
在大多数 programs/terminals 中更可取;因此,displays
被注释掉并使用 print(repr())
。缩短后的代码如下:
## 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)))]).doit()
fxdot_m = Matrix([0 - 3700*(0.02732 - x_m) * q_3c]).doit()
## 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]).doit()
fPValve = Matrix([k_spr * (P_sp - x_p) - rho * g * h_out * a_p + m_p * g - 0]).doit()
## 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]).doit()
fq_m = Matrix([C_vm * (sqrt(((h_in - h_out)))) - q_m]).doit()
fq_1 = Matrix([C_vfo * (sqrt((h_in - h_t))) - q_1]).doit()
fC_vp = Matrix([0.0000753 * (1 - e(-1135 * x_p)) - C_vp]).doit()
fq_2 = Matrix([C_vp * (sqrt((h_t - h_out))) - q_2]).doit()
fmix = Matrix([q_1 + q_3c - q_2]).doit()
fq_3c = Matrix([h_c - h_t]).doit()
## 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)]).doit()
fh_L1 = Matrix([h_L1 - (f*(L1/D1)*((v1**2)/(2*g)))]).doit()
#AlphaValve
fCV1q_m = Matrix([C_valpha * a_alpha * (sqrt((h_pump - h_L1 - (h_in + h_L2)))) - q_m]).doit()
#Pipe2
fp2q_m = Matrix([q_m - (A2 * v2)]).doit()
fh_L2 = Matrix([h_L2 - (f*(L2/D2)*((v2**2)/(2*g)))]).doit()
#Pipe3
fp3q_m = Matrix([q_m - A3 * v3]).doit()
fh_L3 = Matrix([h_L3 - (f*(L3/D3)*((v3**2)/(2*g)))]).doit()
#BetaValve
fCV2q_m = Matrix([C_vbeta * a_beta * (sqrt((h_out - h_L3 - (h_atm + h_L4)))) - q_m]).doit()
#Pipe4
fp4q_m = Matrix([q_m - (A4 * v4)]).doit()
fh_L4 = Matrix([h_L4 - (f*(L4/D4)*((v4**2)/(2*g)))]).doit()
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])).doit()
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]]).doit()
#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])).doit()
#display(B2guessmatrix)
#display(B2FunctionMatrix)
B2iter_n = 0
B2tolerance = (Matrix.ones(B2FunctionMatrix.shape[0], 1) * tolerance).doit()
B2guessVariance = (Matrix.ones(B2FunctionMatrix.shape[0], 1) * 10).doit()
B2JacobianMatrix = (B2FunctionMatrix.jacobian(B2VariableMatrix)).doit()
#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])).doit()
#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])).doit()
#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])).doit()
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')
print('These are the two matrices which are compared in the while loop')
#display(B2guessVariance)
#display(B2tolerance)
print(repr(B2guessVariance))
print(repr(B2tolerance))
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.')
显示while循环中用于迭代的关系矩阵的输出如下:
*PROBLEM* SECTION 7.2.2 BLOCK 2 (MULTIPLE UNKNOWN EQUATIONS)
These are the two matrices which are compared in the while loop
Matrix([
[ Abs(0.00392183668209074*tau0 - 8.02180713197279e-7)],
[ Abs(0.00566975402880878*tau0 - 4.35392954981073e-6)],
[ Abs(7.3945310403193e-25*tau0 - 2.5781512492912e-24)],
[ Abs(4.81682791725627*tau0 - 0.000426049277400864)],
[ Abs(0.617506604286677*tau0 - 0.000126658199850027)],
[ Abs(4.04275751354693*tau0 - 0.00184684018329051)],
[ Abs(0.007854*tau0 - 5.33999999999951e-7)],
[Abs(0.000169107863923268*tau0 - 4.93212155862604e-8)],
[ Abs(0.00263411194406529*tau0 - 7.57176496023193e-7)],
[ Abs(7.41675230154456e-6*tau0 - 3.67497300027552e-9)],
[ Abs(4.04275751354693*tau0 - 0.00184684018329051)],
[ Abs(3.54463776960747e-6*tau0 - 3.53001248327642e-9)],
[ Abs(3.54463776960747e-6*tau0 - 3.53001248327642e-9)],
[ Abs(0.109449541284404*tau0 - 4.00458715596186e-5)],
[ Abs(0.291865443425076*tau0 - 0.000106788990825613)],
[ Abs(0.364831804281346*tau0 - 3.34862385320545e-5)],
[ Abs(0.481577981651376*tau0 - 0.000116201834862384)],
[ 1.0*Abs(tau0)],
[ 1.0*Abs(tau0)],
[ 1.0*Abs(tau0)],
[ Abs(tau0)]])
Matrix([
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05],
[0.05]])
1 iterations have been completed so far. Moving onto the next iteration...
这给出了错误:
TypeError Traceback (most recent call last)
<ipython-input-1-e583c0741f11> in <module>
145 #display(B2JacobianMatrix)
146
--> 147 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]):
148 #print('placeholder 1 (displays B2fsolve, the function matrix with guess values substituted)')
149 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])).doit()
C:\ProgramData\Anaconda3\lib\site-packages\sympy\core\relational.py in __bool__(self)
396
397 def __bool__(self):
--> 398 raise TypeError("cannot determine truth value of Relational")
399
400 def _eval_as_set(self):
TypeError: cannot determine the truth value of Relational
显然,由于tau0符号在输出的B2guessVariance
矩阵中,因此while循环继续迭代时,不能与容差矩阵相关联。有没有一种方法可以删除 tau0 变量,以便求解器可以正常执行并且我可以解决这个问题?
After some research,这似乎可以通过用 params
return 行预先指定 gauss_jordan_solve
然后用零替换 taus 来解决。
因此,通过更改以下内容解决了这个问题:
B2delta_x0, fv = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
至:
B2delta_x0, params = B2jsolve.gauss_jordan_solve(-1*B2fsolve)
B2delta_x0 = B2delta_x0.xreplace({ tau:0 for tau in params })
不幸的是,我并不真正理解为什么 这会起作用,这令人沮丧。但它目前正在运行,并为我提供了很好的答案,这太棒了!