I am trying to convert my program from MATLAB to Python. Everything seems to work fine until I try to invert one of my matrices (22x22). I get different values for the inverse matrix. In MATLAB I have used the function inv(J) and also A\b to invert the matrix, while in python I have used J.I and also np.linalg.inv(J). While both functions in python return the same result, they are different from from what I get in MATLAB. Even the determinant values are different for the two cases. While MATLAB returns a value in the order of 10^23, python returns a value of 10^20. Is it possible for a matrix to have more than one inverse? Or does it have anything to do with the precision of the calculation?
There is a similar discussion in this thread: Discrepancies between Matlab and Numpy matrix inverse. However, I did not find any conclusive answer in that discussion and hence I am posting it again.
Edit: The matrix is not very well conditioned.
Here is the python program:
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)
This is the MATLAB code:
nbus = 14; % IEEE-14, IEEE-30, IEEE-57..
Y = ybusppg(nbus); % Calling ybusppg.m to get Y-Bus Matrix..
busd = busdatas(nbus); % Calling busdatas..
BMva = 100; % Base MVA..
bus = busd(:,1); % Bus Number..
type = busd(:,2); % Type of Bus 1-Slack, 2-PV, 3-PQ..
V = busd(:,3); % Specified Voltage..
del = busd(:,4); % Voltage Angle..
Pg = busd(:,5)/BMva; % PGi..
Qg = busd(:,6)/BMva; % QGi..
Pl = busd(:,7)/BMva; % PLi..
Ql = busd(:,8)/BMva; % QLi..
Qmin = busd(:,9)/BMva; % Minimum Reactive Power Limit..
Qmax = busd(:,10)/BMva; % Maximum Reactive Power Limit..
P = Pg - Pl; % Pi = PGi - PLi..
Q = Qg - Ql; % Qi = QGi - QLi..
Psp = P; % P Specified..
Qsp = Q; % Q Specified..
G = real(Y); % Conductance matrix..
B = imag(Y); % Susceptance matrix..
pv = find(type == 2 | type == 1); % PV Buses..
pq = find(type == 3); % PQ Buses..
npv = length(pv); % No. of PV buses..
npq = length(pq); % No. of PQ buses..
Tol = 1;
Iter = 1;
while (Tol > 1e-5) % Iteration starting..
P = zeros(nbus,1);
Q = zeros(nbus,1);
% Calculate P and Q
for i = 1:nbus
for k = 1:nbus
P(i) = P(i) + V(i)* V(k)*(G(i,k)*cos(del(i)-del(k)) + B(i,k)*sin(del(i)-del(k)));
Q(i) = Q(i) + V(i)* V(k)*(G(i,k)*sin(del(i)-del(k)) - B(i,k)*cos(del(i)-del(k)));
end
end
% Checking Q-limit violations..
if Iter <= 7 && Iter > 2 % Only checked up to 7th iterations..
for n = 2:nbus
if type(n) == 2
QG = Q(n)+Ql(n);
if QG < Qmin(n)
V(n) = V(n) + 0.01;
elseif QG > Qmax(n)
V(n) = V(n) - 0.01;
end
end
end
end
% Calculate change from specified value
dPa = Psp-P;
dQa = Qsp-Q;
k = 1;
dQ = zeros(npq,1);
for i = 1:nbus
if type(i) == 3
dQ(k,1) = dQa(i);
k = k+1;
end
end
dP = dPa(2:nbus);
M = [dP; dQ]; % Mismatch Vector
% Jacobian
% J1 - Derivative of Real Power Injections with Angles..
J1 = zeros(nbus-1,nbus-1);
for i = 1:(nbus-1)
m = i+1;
for k = 1:(nbus-1)
n = k+1;
if n == m
for n = 1:nbus
J1(i,k) = J1(i,k) + V(m)* V(n)*(-G(m,n)*sin(del(m)-del(n)) + B(m,n)*cos(del(m)-del(n)));
end
J1(i,k) = J1(i,k) - V(m)^2*B(m,m);
else
J1(i,k) = V(m)* V(n)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)-del(n)));
end
end
end
% J2 - Derivative of Real Power Injections with V..
J2 = zeros(nbus-1,npq);
for i = 1:(nbus-1)
m = i+1;
for k = 1:npq
n = pq(k);
if n == m
for n = 1:nbus
J2(i,k) = J2(i,k) + V(n)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n)));
end
J2(i,k) = J2(i,k) + V(m)*G(m,m);
else
J2(i,k) = V(m)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n)));
end
end
end
% J3 - Derivative of Reactive Power Injections with Angles..
J3 = zeros(npq,nbus-1);
for i = 1:npq
m = pq(i);
for k = 1:(nbus-1)
n = k+1;
if n == m
for n = 1:nbus
J3(i,k) = J3(i,k) + V(m)* V(n)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n)));
end
J3(i,k) = J3(i,k) - V(m)^2*G(m,m);
else
J3(i,k) = V(m)* V(n)*(-G(m,n)*cos(del(m)-del(n)) - B(m,n)*sin(del(m)-del(n)));
end
end
end
% J4 - Derivative of Reactive Power Injections with V..
J4 = zeros(npq,npq);
for i = 1:npq
m = pq(i);
for k = 1:npq
n = pq(k);
if n == m
for n = 1:nbus
J4(i,k) = J4(i,k) + V(n)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)-del(n)));
end
J4(i,k) = J4(i,k) - V(m)*B(m,m);
else
J4(i,k) = V(m)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)-del(n)));
end
end
end
J = [J1 J2; J3 J4]; % Jacobian Matrix..
condJ = cond(J)
detJ = det(J);
X = J\M;% Correction Vector
dTh = X(1:nbus-1); % Change in Voltage Angle..
dV = X(nbus:end); % Change in Voltage Magnitude..
% Updating State Vectors..
del(2:nbus) = dTh + del(2:nbus); % Voltage Angle..
k = 1;
for i = 2:nbus
if type(i) == 3
V(i) = dV(k) + V(i); % Voltage Magnitude..
k = k+1;
end
end
Iter = Iter + 1;
Tol = max(abs(M)); % Tolerance..
end
The ybusppg function:
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);
The linedata function:
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
The bus data function:
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
invis a bad choice. Just useA\binstead ofinv(A)*b. This happens because you have a badly conditioned matrix. - Ander Biguriinvin a computer, I did not mean in a specific programming language. Usenumpy.linalg.solvefor numpy. - Ander Biguri113.3328, while in python it is3948.1338. - gaurav