Matrix dimension must agree HELP

6 views (last 30 days)
Ahmad Muaz
Ahmad Muaz on 7 Sep 2021
Commented: Sunil on 15 Jan 2024
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')

Sign in to answer this question.

Answers (1)

Jan
Jan on 7 Sep 2021
Edited: Jan on 7 Sep 2021
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
Show 1 older comment
Jan
Jan on 7 Sep 2021
Edited: Jan 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
Sunil on 15 Jan 2024
When I am applying this to IEEE 39 bus system , then it is showing NAN values

Sign in to comment.

Categories

Find more on Semiconductors in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!