我正在尝试将我的程序从MATLAB转换成Python。在我尝试反转矩阵(22x22)之前,一切似乎都很正常。我得到逆矩阵的不同值。在MATLAB中,我使用了函数inv(J)
和{J.I
和{
本帖中也有类似的讨论: Discrepancies between Matlab and Numpy matrix inverse。 然而,我在讨论中没有找到任何结论性的答案,因此我再次发帖。在
编辑:矩阵条件不是很好。在
下面是python程序:
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)
这是MATLAB代码:
^{pr2}$ybusppg函数:
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);
linedata函数:
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
总线数据功能:
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
更详细的回答: 当矩阵条件不好时,这意味着它有一些非常小的特征值。让我们使用奇异值分解,如果矩阵是
X=L @ S @ R.T
,那么逆矩阵是Xi=R @ np.sqrt(S) @ L.T
。中间的部分是对特征值进行反演,如果有一些非常小的特征值,浮点误差在它们的逆函数中会有很大的差异,因此你可以看到它们之间的差异。在但这有关系吗?不,如果你的矩阵条件不好,那么你计算的逆矩阵是不可靠的,你可以用它乘以你原来的矩阵来测试,看看不是每个元素都接近单位矩阵。在
解决办法是什么? 1使用前置条件 2如果你看到原始矩阵的特征值发生了很大的变化,只要扔掉那些不重要的特征值和特征向量,重新构造矩阵,那么不管你使用的平台/软件包/操作系统是什么,你的逆矩阵都应该是一样的。 请注意,第一种方法与第二种方法基本相同,只是数学步骤不同。在
没有一个矩阵不能有一个以上的逆,但它可以没有逆,数学上。在
但是在计算机上进行的所有计算都是用有限的浮点表示法进行的,并且在计算的每一步都可能出现各种各样的误差。在
矩阵数值反演的结果取决于矩阵的“良好条件”。 即使一个矩阵可以倒过来,如果它没有很好的条件,那么输入的小变化也会导致输出的大变化。因为矩阵中的值是用有限的浮点表示法表示的,内部表示法稍有不同就可能造成很大的差异。另外,由于所有的计算都是由有限的浮点表示来执行的,所以所使用的算法的细节可能会导致与病态矩阵的差异。在
相关问题 更多 >
编程相关推荐