FreeMat和NumPy中的矩陣求逆給出了不同的輸出

Question

我完全困惑。我試着在FreeMat如下：

Y1 = [0.1637 - 0.4288i, -0.0460 + 0.1117i, -0.1165 + 0.2828i,  0; 
    -0.0460 + 0.1117i, 0.1090 - 0.3526i, -0.0583 + 0.1414i,  0; 
    -0.1165 + 0.2828i, -0.0583 + 0.1414i, 0.5422 - 1.0004i, -0.3663 + 0.5495i; 
     0,     0,   -0.3663 + 0.5495i, 0.3663 - 0.5495i] 

Z1 = inv(Y1) 

輸出是：

Z1 = 
0.9263 + 7.9980i -0.0021 + 5.4516i 0.4569 + 6.7866i 0.4569 + 6.7866i 
-0.0021 + 5.4516i 0.4434 + 6.6566i 0.0342 + 5.5501i 0.0342 + 5.5501i 
0.4569 + 6.7866i 0.0342 + 5.5501i 0.9363 + 7.9718i 0.9363 + 7.9718i 
0.4569 + 6.7866i 0.0342 + 5.5501i 0.9363 + 7.9718i 1.7763 + 9.2318i 

我試着在NumPy同樣的事情：

Y1 = np.matrix( ((0.1637 - 0.4288j, -0.0460 + 0.1117j, -0.1165 + 0.2828j,  0), 
    (-0.0460 + 0.1117j, 0.1090 - 0.3526j, -0.0583 + 0.1414j,  0), 
    (-0.1165 + 0.2828j, -0.0583 + 0.1414j, 0.5422 - 1.0004j, -0.3663 + 0.5495j), 
    (  0,     0,   -0.3663 + 0.5495j, 0.3663 - 0.5495j))) 
print np.linalg.inv(Y1) 

輸出是：

[[ 0.86953959+6.88397793j 0.00715465+4.54141578j 0.32473444+5.47695118j 
    0.32473444+5.47695118j] 
[ 0.00715465+4.54141578j 0.49628324+5.91573859j -0.00822597+4.47684505j 
    -0.00822597+4.47684505j] 
[ 0.32473444+5.47695118j -0.00822597+4.47684505j 0.70375023+6.43596174j 
    0.70375023+6.43596174j] 
[ 0.32473444+5.47695118j -0.00822597+4.47684505j 0.70375023+6.43596174j 
    1.54375023+7.69596174j]] 

爲什麼輸出不同？

關於第二個想法，在我看來，這個問題可能是在源代碼中的其他地方，所以這裏是Python程序：

#ROEPER SHORT CIRCUIT CALCULATION EXAMPLE 
import numpy as np 

#Constants 
a = np.exp(2*np.pi/3); 
A = np.matrix(((1., 1., 1.), (1., a**2., a), (1., a, a**2.)), dtype=np.complex64) 

#System description at Bus 1 
S_NG1 = 100. #generator rated apparent power in MVA 
U_NG1 = 10.5 #generator rated voltage in kV 
xst_dG1 = 10.5 #relative initial subtransient reactance in percent 
r_G1 = 0.3 #generator resistance in percent 
S_NT1 = 100. # transformer rated capacity in MVA 
U_NT1LV = 10.5 #in kV 
U_NT1HV = 115. #in kV 
u_kNT1 = 11.5 #in percent 
u_RNT1 = 0.5 #in percent 
n_ZT1 = 0.8 #ratio of zero to positive sequence impedance 
#Sequence impedance calculations ----------------------- 
Z1_G1 = complex(r_G1, xst_dG1)*(1./100.)*((U_NG1**2)/S_NG1)*(U_NT1HV/U_NT1LV)**2 
#print Z1_G1 
Z2_G1 = Z1_G1 
Z1_T1 = complex(u_RNT1/100., u_kNT1/100.)*((U_NT1HV**2)/S_NT1) 
#print Z1_T1 
Z2_T1 = Z1_T1 
Z0_T1 = n_ZT1*Z1_T1 

#System description at Bus 2 
S_NT2 = 200. # transformer rated capacity in MVA 
U_NT2LV = 115. #transformer LV side nominal voltage in kV 
U_NT2HV = 230. #in kV 
u_kNT2 = 12. #in percent 
u_RNT2 = 0.3 #in percent 
n_ZT2 = 2.4 #ratio of zero to positive sequence impedance 
#Sequence impedance calculations ------------------ 
Z1_T2 = complex(u_RNT2/100, u_kNT2/100.)*((U_NT2LV**2)/S_NT2) 
#print Z1_T2 
Z2_T2 = Z1_T2 
Z0_T2 = n_ZT2*Z1_T2 

#System description at Bus 3 
S_NG3 = 75. #generator rated apparent power in MVA 
U_NG3 = 10.5 #generator rated voltage in kV 
xst_dG3 = 11.2 #relative initial subtransient reactance in percent 
r_G3 = 0.3 #generator resistance in percent 
S_NT3 = 75. # transformer rated capacity in MVA 
U_NT3LV = 10.5 #in kV 
U_NT3HV = 115. #in kV 
u_kNT3 = 10. #in percent 
u_RNT3 = 0.6 #in percent 
n_ZT3 = 0.8 
Sst_kQ = 12000. #initial a.c. short circuit power at 2Q in MVA 
U_NQ = 230. #rated voltage of the system at point of connection 2Q in kV 
U_Q = 220. 
U_N2 = 115. 
RX_Q = 0.1 
#Sequence impedance calculations ------------------ 
Z1_G3 = complex(r_G3, xst_dG3)*(1/100)*((U_NG3**2)/S_NG3)*(U_NT3HV/U_NT3LV)**2 
Z2_G3 = Z1_G3 
Z1_T3 = complex(u_RNT3/100, u_kNT3/100)*((U_NT3HV**2)/S_NT3) 
Z2_T3 = Z1_T3 
Z0_T3 = n_ZT3*Z1_T3 
Z1Q = complex(0.1, 1.0)*1.1*((U_Q**2)/Sst_kQ)*(U_N2/U_NQ)**2 
Z2Q = Z1Q 

#Bus interconnections 
Z1_12 = complex(3.154, 7.657) #in ohms 
Z2_12 = Z1_12 #in ohms 
Z0_12 = complex(10.431, 31.673) #in ohms 
G0_12 = complex(0, 1./(1.056*(10.**4.))) #in mhos 

Z1_23 = complex(2.490, 6.045) #in ohms 
Z2_23 = Z1_23 #in ohms 
Z0_23 = complex(8.235, 25.005) #in ohms 
G0_23 = complex(0, 1./(1.337*(10.**4.))) #in mhos 

Z1_13 = complex(1.245, 3.023) #in ohms 
Z2_13 = Z1_13 #in ohms 
Z0_13 = complex(4.118, 12.503) #in ohms 
G0_13 = complex(0, 1./(2.674*(10.**4.))) #in mhos 

Z1_2Q = complex(0.173, 0.980) #in ohms 
Z2_2Q = Z1_2Q 

Z1_34 = complex(0.840, 1.260) #in ohms 
Z2_34 = Z1_34 
Z0_34 = complex(5.440, 3.430) 
G0_34 = complex(0, 1./(1.098*(10.**3.))) #in mhos 

#Creation of Y matrix for positive & negative sequence components 
y1_11 = 1/(Z1_G1+Z1_T1)+1/Z1_12+1/Z1_13; 
y1_12 = -1/Z1_12; 
y1_13 = -1/Z1_13; 
y1_14 = 0+0j; 
y1_21 = y1_12; 
y1_22 = 1/(Z1_T2+Z1_2Q+Z1Q)+1/Z1_12+1/Z1_23; 
y1_23 = -1/Z1_23; 
y1_24 = 0+0j; 
y1_31 = y1_13; 
y1_32 = y1_23; 
y1_33 = 1/(Z1_G3+Z1_T3)+1/Z1_13+1/Z1_23+1/Z1_34; 
y1_34 = -1/Z1_34; 
y1_41 = 0+0j; 
y1_42 = 0+0j; 
y1_43 = y1_34; 
y1_44 = 1/Z1_34; 


Y1 = np.matrix(((y1_11, y1_12, y1_13, y1_14), 
     (y1_21, y1_22, y1_23, y1_24), 
     (y1_31, y1_32, y1_33, y1_34), 
     (y1_41, y1_42, y1_43, y1_44))) 

Z1 = np.linalg.inv(Y1); 
Ist_k3 = 1.1*110./(np.sqrt(3)*Z1[3,3]) 
If3ph = abs(Ist_k3) 
print If3ph 
phi_3ph = np.angle(Ist_k3, deg=True) 
print phi_3ph 

Answer 1

有些矩陣是不穩定的數值運算，喜歡你的矩陣Y1。矩陣的數字行爲的特點是其condition number。條件數字接近1是好的，cond(Y1)給出45.

例如參見numpy.linalg.cond。