Matrix dimension must agree HELP
Show older comments
This is my code
this part said matrix dimension must agree (X = inv(J1).*M; % Correction Vector)
clc;
clear all;
Vbase = 230*10^3;
Sbase = 100*10^6;
Zbase = Vbase^2/Sbase;
Ybase = 1/Zbase;
% | From | To | R | X | B/2 | X'mer |
% | Bus | Bus | pu | pu | pu | TAP (a) |
linedata= [ 1 2 10.58*5/Zbase 21.16i*5/Zbase 2.11i*5/Ybase 1.0;
1 4 8.82*3/Zbase 35.27i*3/Zbase 3.52i*3/Ybase 1.0;
2 3 1.76*15/Zbase 8.82i*15/Zbase 1.06i*15/Ybase 1.0;
2 4 0.88*30/Zbase 1.76i*30/Zbase 0.17i*30/Ybase 1.0;
3 6 0.53*20/Zbase 2.65i*20/Zbase 0.26i*20/Ybase 1.0;
4 5 4.23*25/Zbase 8.46i*25/Zbase 0.84i*25/Ybase 1.0;
5 6 1.32*40/Zbase 3.96i*40/Zbase 0.40i*40/Ybase 1.0];
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 + x; % z matrix...
y = 1./z; % To get inverse of each element...
b = 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(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); % Bus Impedance Matrix
% |Bus | Type | Vsp | theta | PGi | QGi | Qmin | Qmax | PLi | QLi |
busd = [1 1 1.07 0 0.0 0 -300*10^6 500*10^6 0 0;
2 2 1.04 0 100*10^6 0 -200*10^6 300*10^6 0 0;
3 2 1.05 0 125*10^6 0 -100*10^6 200*10^6 0 0;
4 3 1.00 0 0.0 0 0 0 120*10^6 25*10^6;
5 3 1.00 0 0.0 0 0 0 80*10^6 15*10^6;
6 3 1.00 0 0.0 0 0 0 60*10^6 20*10^6];
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)/Sbase; % PGi..
Qg = busd(:,6)/Sbase; % QGi..
Pl = busd(:,9)/Sbase; % PLi..
Ql = busd(:,10)/Sbase; % QLi..
Qmin = busd(:,7)/Sbase; % Minimum Reactive Power Limit..
Qmax = busd(:,8)/Sbase; % Maximum Reactive Power Limit..
Psp = Pg - Pl; % Pspecified = PGi - PLi..
Qsp = Qg - Ql; % Qspecified = QGi - QLi..
G = real(Y); % Conductance matrix..
B = imag(Y); % Susceptance matrix..
pv = find(type == 2); % 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..
nbus = nb;
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
Iter
P
Q
% 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);
dQ = dQa(3:nbus);
M = [dP; dQ]; % Mismatch Vector
% Jacobian
% H - Derivative of Real Power Injections with Angles..
H = zeros(nbus-1,nbus-1);
for i = 1:(nbus-1)
m = i+1;
for k = 1:(nbus-1)
n = k+1;
if n == m
H(i,k) = Q(m) + V(m)^2*B(m,m);
else
H(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
% J - Derivative of Reactive Power Injections with Angles..
J = zeros(npq,nbus-1);
for i = 1:npq
m = pq(i);
for k = 1:(nbus-1)
n = k+1;
if n == m
J(i,k) = -P(m) + V(m)^2*G(m,m);
else
J(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
% N - Derivative of Real Power Injections with V..
N = zeros(nbus-1,npq);
for i = 1:(nbus-1)
m = i+1;
for k = 1:npq
n = pq(k);
if n == m
N(i,k) = -P(m) -V(m)^2*G(m,m);
else
N(i,k) = -J(k,i);
end
end
end
% L - Derivative of Reactive Power Injections with V..
L = zeros(npq,npq);
for i = 1:npq
m = pq(i);
for k = 1:npq
n = pq(k);
if n == m
L(i,k) = -Q(m) + V(m)^2*B(m,m);
else
L(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
J1 = [H N; J L]; % Jacobian Matrix..
size(J1)
size(M)
X = inv(J1).*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
V (3:nbus) = -dV.*V (3:nbus) + V (3: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
Iter = Iter -1;
Voltage = real(V);
Delta = real (del.*180/pi);
nb=nbus;
pol2rect = V.*cos(del) + 1i*V.*sin(del);
Y; % Calling Ybus function..
%lined = linedata(nb); % Get linedats..
%busd = bus(nb); % Get busdatas..
Vm = pol2rect; % Converting polar to rectangular..
Del = 180/pi*del; % Bus Voltage Angles in Degree...
fb = linedata(:,1); % From bus number...
tb = linedata(:,2); % To bus number...
nl = length(fb); % No. of Branches..
Pl = busd(:,9); % PLi..
Ql = busd(:,10); % QLi..
Iij = zeros(nb,nb);
Sij = zeros(nb,nb);
Si = zeros(nb,1);
% Bus Current Injections..
I = Y*Vm;
Im = abs(I);
Ia = angle(I);
%Line Current Flows..
for m = 1:nl
p = fb(m); q = tb(m);
Iij(p,q) = -(Vm(p) - Vm(q))*Y(p,q); % Y(m,n) = -y(m,n)..
Iij(q,p) = -Iij(p,q);
end
%Iij = sparse(Iij); % Commented out..
Iijm = abs(Iij);
Iija = angle(Iij);
% Line Power Flows..
for m = 1:nb
for n = 1:nb
if m ~= n
Sij(m,n) = Vm(m)*conj(Iij(m,n))*Sbase;
end
end
end
%Sij = sparse(Sij); % Commented out..
Pij = real(Sij);
Qij = imag(Sij);
% Line Losses..
Lij = zeros(nl,1);
for m = 1:nl
p = fb(m); q = tb(m);
Lij(m) = Sij(p,q) + Sij(q,p);
end
Lpij = real(Lij);
Lqij = imag(Lij);
% Bus Power Injections..
for i = 1:nb
for k = 1:nb
Si(i) = Si(i) + conj(Vm(i))* Vm(k)*Y(i,k)*Sbase;
end
end
Pi = real(Si);
Qi = -imag(Si);
Pg = Pi+Pl;
Qg = Qi+Ql;
disp('#########################################################################################');
disp('-----------------------------------------------------------------------------------------');
disp(' Newton Raphson Loadflow Analysis ');
disp('-----------------------------------------------------------------------------------------');
disp('| Bus | V | Angle | Injection | Generation | Load |');
disp('| No | pu | Degree | MW | MVar | MW | Mvar | MW | MVar | ');
for m = 1:nb
disp('-----------------------------------------------------------------------------------------');
fprintf('%3g', m); fprintf(' %8.4f', V(m)); fprintf(' %8.4f', Del(m));
fprintf(' %8.3f', Pi(m)/10^6); fprintf(' %8.3f', Qi(m)/10^6);
fprintf(' %8.3f', Pg(m)/10^6); fprintf(' %8.3f', Qg(m)/10^6);
fprintf(' %8.3f', Pl(m)/10^6); fprintf(' %8.3f', Ql(m)/10^6); fprintf('\n');
end
disp('-----------------------------------------------------------------------------------------');
fprintf(' Total ');fprintf(' %8.3f', sum(Pi)/10^6); fprintf(' %8.3f', sum(Qi)/10^6);
fprintf(' %8.3f', sum(Pi+Pl)/10^6); fprintf(' %8.3f', sum(Qi+Ql)/10^6);
fprintf(' %8.3f', sum(Pl)/10^6); fprintf(' %8.3f', sum(Ql)/10^6); fprintf('\n');
disp('-----------------------------------------------------------------------------------------');
disp('#########################################################################################');
disp('-------------------------------------------------------------------------------------');
disp(' Line FLow and Losses ');
disp('-------------------------------------------------------------------------------------');
disp('|From|To | P | Q | From| To | P | Q | Line Loss |');
disp('|Bus |Bus| MW | MVar | Bus | Bus| MW | MVar | MW | MVar |');
for m = 1:nl
p = fb(m); q = tb(m);
disp('-------------------------------------------------------------------------------------');
fprintf('%4g', p); fprintf('%4g', q); fprintf(' %8.3f', Pij(p,q)/10^6); fprintf(' %8.3f', Qij(p,q)/10^6);
fprintf(' %4g', q); fprintf('%4g', p); fprintf(' %8.3f', Pij(q,p)/10^6); fprintf(' %8.3f', Qij(q,p)/10^6);
fprintf(' %8.3f', Lpij(m)/10^6); fprintf(' %8.3f', Lqij(m)/10^6);
fprintf('\n');
end
disp('-------------------------------------------------------------------------------------');
fprintf(' Total Loss ');
fprintf(' %8.3f', sum(Lpij)/10^6); fprintf(' %8.3f', sum(Lqij)/10^6); fprintf('\n');
disp('-------------------------------------------------------------------------------------');
disp('#####################################################################################');
subplot(211)
bar(V);
title('bus voltage magnitude')
xlabel('bus number')
ylabel('voltage in pu')
subplot(212)
bar((180/pi)*del);
title('bus voltage phase')
xlabel('bus number')
ylabel('phase in degrees')
Answers (1)
Are you sure, that you mean the elementwise multiplication and not a matrix multiplication?
inv(J1) .* M % Elementwise multiplication
inv(J1) * M % Matrix multiplication
By the way, the multiplication by the inverse is slower and less stable than the matrix division:
J1 \ M % Instead of inv(J1) * M
3 Comments
Ahmad Muaz
on 7 Sep 2021
Okay. Then use the debugger to stop, when the error occurs:
dbstop if error
Run the code again. When Matlab stops, check the sizes:
size(J1)
size(M)
What do you get?
I get 8x8 and 9x1. This cannot work.
Due to the lack of useful comments, I do not have a chance to find out, where the problem is in the definition of M:
dP = dPa(2:nbus);
dQ = dQa(3:nbus);
M = [dP; dQ]; % Mismatch Vector
Sunil
on 15 Jan 2024
When I am applying this to IEEE 39 bus system , then it is showing NAN values
Categories
Find more on RF Network Construction in Help Center and File Exchange
on 7 Sep 2021
on 15 Jan 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
