回答此问题可获得 20 贡献值,回答如果被采纳可获得 50 分。
<p>我正在尝试将我的程序从MATLAB转换成Python。在我尝试反转矩阵(22x22)之前,一切似乎都很正常。我得到逆矩阵的不同值。在MATLAB中,我使用了函数<code>inv(J)</code>和{<cd2>}来反转矩阵,而在python中我使用了<code>J.I</code>和{<cd4>}。虽然python中的两个函数返回相同的结果,但它们与我在MATLAB中得到的结果不同。两种情况下的行列式值也不同。MATLAB返回10^23的值,python则返回10^20的值。一个矩阵有可能有多个逆矩阵吗?还是跟计算精度有关?在</p>
<p>本帖中也有类似的讨论:
<a href="https://stackoverflow.com/questions/39434302/discrepancies-between-matlab-and-numpy-matrix-inverse">Discrepancies between Matlab and Numpy matrix inverse</a>。
然而,我在讨论中没有找到任何结论性的答案,因此我再次发帖。在</p>
<p>编辑:矩阵条件不是很好。在</p>
<p>下面是python程序:</p>
<pre><code>import numpy as np
import pandas as pd
import cmath
import math
from decimal import *
from pandas import DataFrame
ldata = '''\
1 2 1 1 1 0 0.01938 0.05917 0.0528 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1 5 1 1 1 0 0.05403 0.22304 0.0492 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 3 1 1 1 0 0.04699 0.19797 0.0438 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 4 1 1 1 0 0.05811 0.17632 0.0340 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 5 1 1 1 0 0.05695 0.17388 0.0346 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
3 4 1 1 1 0 0.06701 0.17103 0.0128 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
4 5 1 1 1 0 0.01335 0.04211 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
4 7 1 1 1 0 0.0 0.20912 0.0 0 0 0 0 0 0.978 0.0 0.0 0.0 0.0 0.0 0.0
4 9 1 1 1 0 0.0 0.55618 0.0 0 0 0 0 0 0.969 0.0 0.0 0.0 0.0 0.0 0.0
5 6 1 1 1 0 0.0 0.25202 0.0 0 0 0 0 0 0.932 0.0 0.0 0.0 0.0 0.0 0.0
6 11 1 1 1 0 0.09498 0.19890 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
6 12 1 1 1 0 0.12291 0.25581 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
6 13 1 1 1 0 0.06615 0.13027 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
7 8 1 1 1 0 0.0 0.17615 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
7 9 1 1 1 0 0.0 0.11001 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
9 10 1 1 1 0 0.03181 0.08450 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
9 14 1 1 1 0 0.12711 0.27038 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
10 11 1 1 1 0 0.08205 0.19207 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
12 13 1 1 1 0 0.22092 0.19988 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
13 14 1 1 1 0 0.17093 0.34802 0.0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
'''
line = pd.compat.StringIO(ldata)
df = pd.read_csv(line, sep='\s+', header=None)
linedata = df.values
nbus = 14
fb = (linedata[:,0]) #from bus
fb = fb.astype(int)
tb = (linedata[:,1]) #to bus
tb = tb.astype(int)
r = linedata[:,6] #resistence
x = linedata[:,7] #reactance
b = (linedata[:,8])/2 #b/2
b = 1j*b
a = linedata[:,14] #transformer ratio
#change the 0s to 1s for the transformer tap changing ratio
for i in range(len(a)):
if a[i] == 0:
a[i] = 1
z = r + 1j*x
y = 1/z
nb = int(max(max(fb),max(tb)))
#print(max(tb))
nl = len(fb)
#print(nl)
Y = np.zeros((nb,nb))
Y = Y + 1j*Y
#forming the off-diagonal elements
for k in range(nl):
Y[fb[k]-1,tb[k]-1] = (Y[fb[k]-1,tb[k]-1]) - (y[k]/a[k])
Y[tb[k]-1,fb[k]-1] = Y[fb[k]-1,tb[k]-1]
#forming the diagonal elements
for n in range(nl):
Y[fb[n]-1,fb[n]-1] = Y[fb[n]-1,fb[n]-1]+(y[n]/(a[n]*a[n]))+b[n]
#elif tb[n]==m:
Y[tb[n]-1,tb[n]-1]=Y[tb[n]-1,tb[n]-1]+((y[n])+b[n])
#print(Y)
bdata = '''\
1 Bus 1 HV 1 1 3 1.060 0.0 0.0 0.0 0 0 0.0 1.060 0.0 0.0 0.0 0.0 0
2 Bus 2 HV 1 1 2 1.045 -4.98 21.7 12.7 40.0 42.4 0.0 1.045 50.0 -40.0 0.0 0.0 0
3 Bus 3 HV 1 1 2 1.010 -12.72 94.2 19.0 0.0 23.4 0.0 1.010 40.0 0.0 0.0 0.0 0
4 Bus 4 HV 1 1 0 1.019 -10.33 47.8 -3.9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0
5 Bus 5 HV 1 1 0 1.020 -8.78 7.6 1.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0
6 Bus 6 LV 1 1 2 1.070 -14.22 11.2 7.5 0.0 12.2 0.0 1.070 24.0 -6.0 0.0 0.0 0
7 Bus 7 ZV 1 1 0 1.062 -13.37 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0
8 Bus 8 TV 1 1 2 1.090 -13.36 0.0 0.0 0.0 17.4 0.0 1.090 24.0 -6.0 0.0 0.0 0
9 Bus 9 LV 1 1 0 1.056 -14.94 29.5 16.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.19 0
10 Bus 10 LV 1 1 0 1.051 -15.10 9.0 5.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0
11 Bus 11 LV 1 1 0 1.057 -14.79 3.5 1.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0
12 Bus 12 LV 1 1 0 1.055 -15.07 6.1 1.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0
13 Bus 13 LV 1 1 0 1.050 -15.16 13.5 5.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0
14 Bus 14 LV 1 1 0 1.036 -16.04 14.9 5.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0
'''
busd = pd.compat.StringIO(bdata)
dfbus = pd.read_csv(busd, sep='\s+', header=None)
dfbus = dfbus.loc[:,dfbus.dtypes != 'object']
busdata = dfbus.values
#print(busdata)
BMva = 100
bus = busdata[:,1]
type = busdata[:,4]
for i in range(len(type)):
if type[i] == 3:
type[i] = 1
elif type[i] == 0:
type[i] = 3
V = busdata[:,12]
for i in range(len(V)):
if V[i] == 0:
V[i] = 1
delta = busdata[:,15]
Pg = busdata[:,9]/BMva
Qg = busdata[:,10]/BMva
Pl = busdata[:,7]/BMva
Ql = busdata[:,8]/BMva
Qmin = busdata[:,14]/BMva
Qmax = busdata[:,13]/BMva
P = Pg - Pl
Q = Qg - Ql
Psp = P
Qsp = Q
#print(Qsp)
G = Y.real
B = Y.imag
pq = np.where(type == 3)
npq = len(pq[0])
pv = np.where(type < 3)
npv = len(pv[0])
Tol = 1
Iter = 1
while Tol > 1e-5:
P = np.zeros(nbus)
Q = np.zeros(nbus)
#Calculate P and Q
for i in range(nbus):
for k in range(nbus):
P[i] = P[i] + V[i]*V[k]*(G[i,k]*np.cos(delta[i]-delta[k]) + B[i,k]*np.sin(delta[i]-delta[k]))
Q[i] = Q[i] + V[i]*V[k]*(G[i,k]*np.sin(delta[i]-delta[k]) - B[i,k]*np.cos(delta[i]-delta[k]))
#Calculate the change from the specified value
dPa = Psp - P
dQa = Qsp - Q
k = 0
dQ = np.zeros(npq)
for i in range(nbus):
if type[i] == 3:
dQ[k] = dQa[i]
k = k + 1
dP = dPa[1:nbus]
#print(dP)
M = np.concatenate((dP,dQ))
#Jacobian Matrix formation
#J1-Derivative of real power injections with angles
J1 = np.zeros((nbus-1,nbus-1))
for i in range(nbus-1):
m = i + 1
for k in range(nbus-1):
n = k + 1
if n == m:
for n in range(nbus):
J1[i,k] = J1[i,k] + V[m]*V[n]*(-G[m,n]*np.sin(delta[m]-delta[n])+B[m,n]*np.cos(delta[m]-delta[n]))
J1[i,k] = J1[i,k] - (V[m]**2)*B[m,m]
else:
J1[i,k] = V[m]*V[n]*(G[m,n]*np.sin(delta[m]-delta[n])-B[m,n]*np.cos(delta[m]-delta[n]))
#J2 = Derivative of real power injections with V
J2 = np.zeros((nbus-1,npq))
for i in range(nbus-1):
m = i +1
for k in range(npq):
n = pq[0][k]
#print(n)
if n == m:
for n in range(nbus):
J2[i,k] = J2[i,k] + V[n]*(G[m,n]*np.cos(delta[m]-delta[n])+B[m,n]*np.sin(delta[m]-delta[n]))
J2[i,k] = J2[i,k] - V[m]*G[m,m]
else:
J2[i,k] = V[m]*(G[m,n]*np.cos(delta[m]-delta[n])+B[m,n]*np.sin(delta[m]-delta[n]))
#J3 = derivative of Reactive Power Injections with Angles
J3 = np.zeros((npq,nbus-1))
for i in range(npq):
m = pq[0][i]
for k in range(nbus-1):
n = k + 1
if n == m:
for n in range(nbus):
J3[i,k] = J3[i,k] + V[m]*V[n]*(G[m,n]*np.cos(delta[m]-delta[n])+B[m,n]*np.sin(delta[m]-delta[n]))
J3[i,k] = J3[i,k] - (V[m]**2)*G[m,m]
else:
J3[i,k] = V[m]*V[n]*(-G[m,n]*np.cos(delta[m]-delta[n])-B[m,n]*np.sin(delta[m]-delta[n]))
#J4 = Derivative of Reactive Power Injections with V
J4 = np.zeros((npq,npq))
for i in range(npq):
m = pq[0][i]
for k in range(npq):
n = pq[0][k]
if n == m:
for n in range(nbus):
J4[i,k] = J4[i,k] + V[n]*(G[m,n]*np.sin(delta[m]-delta[n])-B[m,n]*np.cos(delta[m]-delta[n]))
J4[i,k] = J4[i,k] - V[m]*B[m,m]
else:
J4[i,k] = V[m]*(G[m,n]*np.sin(delta[m]-delta[n])-B[m,n]*np.cos(delta[m]-delta[n]))
Jtemp1 = np.concatenate((J1,J3))
Jtemp2 = np.concatenate((J2,J4))
J = np.concatenate((Jtemp1,Jtemp2),axis=1)
#Jinvc = getMatrixInverse(J)
X = np.linalg.solve(J,M)
print(X)
dTh = X[0,0:nbus-1]
#print(dTh)
dV = X[0,nbus-1:2*npq+npv-1]
print(dV[0,8])
delta[1:nbus] = dTh + delta[1:nbus]
k = 0
for i in range(1,nbus):
if type[i] == 3:
V[i] = dV[0,k] + V[i]
k = k + 1
Iter = Iter + 1
Tol = max(abs(M))
print(V)
</code></pre>
<p>这是MATLAB代码:</p>
^{pr2}$
<p>ybusppg函数:</p>
<pre class="lang-matlab prettyprint-override"><code>function Y = ybusppg(num) % Returns Y
linedata = linedatas(num); % Calling Linedatas...
fb = linedata(:,1); % From bus number...
tb = linedata(:,2); % To bus number...
r = linedata(:,3); % Resistance, R...
x = linedata(:,4); % Reactance, X...
b = linedata(:,5); % Ground Admittance, B/2...
a = linedata(:,6); % Tap setting value..
z = r + i*x; % z matrix...
y = 1./z; % To get inverse of each element...
b = i*b; % Make B imaginary...
nb = max(max(fb),max(tb)); % No. of buses...
nl = length(fb); % No. of branches...
Y = zeros(nb,nb); % Initialise YBus...
% Formation of the Off Diagonal Elements...
for k = 1:nl
Y(fb(k),tb(k)) = Y(fb(k),tb(k)) - y(k)/a(k);
Y(tb(k),fb(k)) = Y(fb(k),tb(k));
end
% Formation of Diagonal Elements....
for m = 1:nb
for n = 1:nl
if fb(n) == m
Y(m,m) = Y(m,m) + y(n)/(a(n)^2) + b(n);
elseif tb(n) == m
Y(m,m) = Y(m,m) + y(n) + b(n);
end
end
end
%Y; % Bus Admittance Matrix
%Z = inv(Y);
</code></pre>
<p>linedata函数:</p>
<pre class="lang-matlab prettyprint-override"><code>function linedt = linedatas(num)
% | From | To | R | X | B/2 | X'mer |
% | Bus | Bus | pu | pu | pu | TAP (a) |
linedat14 = [1 2 0.01938 0.05917 0.0264 1
1 5 0.05403 0.22304 0.0246 1
2 3 0.04699 0.19797 0.0219 1
2 4 0.05811 0.17632 0.0170 1
2 5 0.05695 0.17388 0.0173 1
3 4 0.06701 0.17103 0.0064 1
4 5 0.01335 0.04211 0.0 1
4 7 0.0 0.20912 0.0 0.978
4 9 0.0 0.55618 0.0 0.969
5 6 0.0 0.25202 0.0 0.932
6 11 0.09498 0.19890 0.0 1
6 12 0.12291 0.25581 0.0 1
6 13 0.06615 0.13027 0.0 1
7 8 0.0 0.17615 0.0 1
7 9 0.0 0.11001 0.0 1
9 10 0.03181 0.08450 0.0 1
9 14 0.12711 0.27038 0.0 1
10 11 0.08205 0.19207 0.0 1
12 13 0.22092 0.19988 0.0 1
13 14 0.17093 0.34802 0.0 1 ];
switch num
case 14
linedt = linedat14;
end
</code></pre>
<p>总线数据功能:</p>
<pre class="lang-matlab prettyprint-override"><code>function busdt = busdatas(num)
% Type....
% 1 - Slack Bus..
% 2 - PV Bus..
% 3 - PQ Bus..
% |Bus | Type | Vsp | theta | PGi | QGi | PLi | QLi | Qmin | Qmax |
busdat14 = [1 1 1.060 0 0 0 0 0 0 0;
2 2 1.045 0 40 42.4 21.7 12.7 -40 50;
3 2 1.010 0 0 23.4 94.2 19.0 0 40;
4 3 1.0 0 0 0 47.8 -3.9 0 0;
5 3 1.0 0 0 0 7.6 1.6 0 0;
6 2 1.070 0 0 12.2 11.2 7.5 -6 24;
7 3 1.0 0 0 0 0.0 0.0 0 0;
8 2 1.090 0 0 17.4 0.0 0.0 -6 24;
9 3 1.0 0 0 0 29.5 16.6 0 0;
10 3 1.0 0 0 0 9.0 5.8 0 0;
11 3 1.0 0 0 0 3.5 1.8 0 0;
12 3 1.0 0 0 0 6.1 1.6 0 0;
13 3 1.0 0 0 0 13.5 5.8 0 0;
14 3 1.0 0 0 0 14.9 5.0 0 0;];
switch num
case 14
busdt = busdat14;
end
</code></pre>