capsmith1985
Electrical
- Nov 12, 2008
- 3
Hi everyone.
I'm trying to modify a code to double-check Newton-Raphson, and I'm having trouble with matrix dimensions. The error reads the following:
Error in ==> nrlfppg2 at 80
X = inv(J)*M; % Correction Vector
I'm not sure how to remedy this, as everything seems alright to me and the original code worked (w/ 30 buses). Any help on this would be greatly appreciated. Here is the code:
Y = ybusppg(); % Calling ybusppg.m to get Bus Admittance Matrix..
busdata = bus5(); % Calling busdata30.m to get busdatas..
bus = busdata
,1); % Bus
type = busdata
,2); % Type of Bus 1-Slack, 2-PV, 3-PQ
V = busdata
,3); % Specified Voltage.
Th = busdata
,4); % Voltage Angle.
P = busdata
,5); % PG
Q = busdata
,6); % QG
n = max(bus); % To get no. of buses
P = P; % Pi = PGi - PLi
Q = Q; % Qi = QGi - QLi
Tol = 10; % Tolerence kept at high value.
Iter = 1; % iteration starting
Psp = P;
Qsp = Q;
G = real(Y); % Conductance..
B = imag(Y); % Susceptance..
while (Tol > 0.00001)
P = zeros(n,1);
Q = zeros(n,1);
% Calculate P and Q
for i = 1:n
for k = 1:n
P(i) = P(i) + V(i)*V(k)*abs(Y(i,k))*cos(Th(i)-Th(k)-angle(Y(i,k)));
Q(i) = Q(i) + V(i)*V(k)*abs(Y(i,k))*sin(Th(i)-Th(k)-angle(Y(i,k)));
end
end
% Calculate change from specified value
dPa = Psp-P;
dQa = Qsp-Q;
k = 1;
for i = 1:n
if type(i) == 3
dQ(k,1) = dQa(i);
k = k+1;
end
end
dP = dPa(2:n);
M = [dP; dQ]; % Mismatch Vector
% Jacobian
for i = 2:n
for k = 2:n
ii = i-1;
kk = k-1;
if k ~= i
J1(ii,kk) = V(i) * V(k) * abs(Y(i,k)) * sin(Th(i)-Th(k)-angle(Y(i,k)));
J2(ii,kk) = V(i) * abs(Y(i,k)) * cos(Th(i)-Th(k)-angle(Y(i,k)));
J3(ii,kk) = -V(i) * V(k) * abs(Y(i,k)) * cos(Th(i)-Th(k)-angle(Y(i,k)));
J4(ii,kk) = V(i) * abs(Y(i,k)) * sin(Th(i)-Th(k)-angle(Y(i,k)));
else
J1(ii,ii) = -Q(i) - V(i)^2*B(i,i);
J2(ii,ii) = P(i)/V(i) + V(i)*G(i,i);
J3(ii,ii) = P(i) - V(i)^2*G(i,i);
J4(ii,ii) = -Q(i)/V(i) - V(i)*B(i,i);
end
end
end
% To remove rows & columns corresponding to PV buses from Jacobian.
r = find(type == 2);
l = length(r);
for i = l:-1:1
J2
,(r(i)-1)) = [];
J4
,(r(i)-1)) = [];
end
for i = l:-1:1
J3((r(i)-1),
= [];
J4((r(i)-1),
= [];
end
J = [J1 J2; J3 J4]; % Jacobian
X = inv(J)*M; % Correction Vector
dTh = X(1:n-1);
dV = X(n:length(X));
Th(2:n) = dTh + Th(2:n);
k = 1;
for i = 2:n
if type(i) == 3
V(i) = dV(k) + V(i);
k = k+1;
end
end
Iter = Iter + 1;
Tol = max(abs(X));
end
Iter % Number of Iterations took..
V;
del = 180/pi*Th;
[V del] % Bus Voltages and angles..
The bus5.m file follows:
function busdata = bus5() % Returns busdata.
busdata = [ 1 1 1.02 0 0.000 0.000;
2 3 1.00 0 -0.60 -0.30;
3 2 1.04 0 1.000 0.000;
4 3 1.00 0 -0.40 -0.10;
5 3 1.00 0 -0.60 -0.20 ];
lines.m follows:
function linedata = lines() % Returns linedata.
linedata = [1 2 0.00 0.40;
1 4 0.00 0.60;
1 5 0.00 0.20;
2 3 0.00 0.20;
2 4 0.00 0.40;
3 5 0.00 0.20;];
ybusppg.m follows:
function ybus = ybusppg(); % Returns ybus
linedata = lines(); % Calling "lines.m" for Line Data...
fb = linedata
,1); % From bus number...
tb = linedata
,2); % To bus number...
r = linedata
,3); % Resistance, R...
x = linedata
,4); % Reactance, X...
z = r + i*x; % Z matrix...
y = 1./z; % To get inverse of each element...
nbus = max(max(fb),max(tb)); % no. of buses...
nbranch = length(fb); % no. of branches...
ybus = zeros(nbus,nbus); % Initialise YBus...
% Formation of the Off Diagonal Elements...
for k=1:nbranch
ybus(fb(k),tb(k)) = -y(k);
ybus(tb(k),fb(k)) = ybus(fb(k),tb(k));
end
% Formation of Diagonal Elements....
for m=1:nbus
for n=1:nbranch
if fb
== m | tb
== m
ybus(m,m) = ybus(m,m) + y
;
end
end
end
%ybus; % Bus Admittance Matrix
%zbus = inv(ybus); % Bus Impedance Matrix
I'm trying to modify a code to double-check Newton-Raphson, and I'm having trouble with matrix dimensions. The error reads the following:
Error in ==> nrlfppg2 at 80
X = inv(J)*M; % Correction Vector
I'm not sure how to remedy this, as everything seems alright to me and the original code worked (w/ 30 buses). Any help on this would be greatly appreciated. Here is the code:
Y = ybusppg(); % Calling ybusppg.m to get Bus Admittance Matrix..
busdata = bus5(); % Calling busdata30.m to get busdatas..
bus = busdata
type = busdata
V = busdata
Th = busdata
P = busdata
Q = busdata
n = max(bus); % To get no. of buses
P = P; % Pi = PGi - PLi
Q = Q; % Qi = QGi - QLi
Tol = 10; % Tolerence kept at high value.
Iter = 1; % iteration starting
Psp = P;
Qsp = Q;
G = real(Y); % Conductance..
B = imag(Y); % Susceptance..
while (Tol > 0.00001)
P = zeros(n,1);
Q = zeros(n,1);
% Calculate P and Q
for i = 1:n
for k = 1:n
P(i) = P(i) + V(i)*V(k)*abs(Y(i,k))*cos(Th(i)-Th(k)-angle(Y(i,k)));
Q(i) = Q(i) + V(i)*V(k)*abs(Y(i,k))*sin(Th(i)-Th(k)-angle(Y(i,k)));
end
end
% Calculate change from specified value
dPa = Psp-P;
dQa = Qsp-Q;
k = 1;
for i = 1:n
if type(i) == 3
dQ(k,1) = dQa(i);
k = k+1;
end
end
dP = dPa(2:n);
M = [dP; dQ]; % Mismatch Vector
% Jacobian
for i = 2:n
for k = 2:n
ii = i-1;
kk = k-1;
if k ~= i
J1(ii,kk) = V(i) * V(k) * abs(Y(i,k)) * sin(Th(i)-Th(k)-angle(Y(i,k)));
J2(ii,kk) = V(i) * abs(Y(i,k)) * cos(Th(i)-Th(k)-angle(Y(i,k)));
J3(ii,kk) = -V(i) * V(k) * abs(Y(i,k)) * cos(Th(i)-Th(k)-angle(Y(i,k)));
J4(ii,kk) = V(i) * abs(Y(i,k)) * sin(Th(i)-Th(k)-angle(Y(i,k)));
else
J1(ii,ii) = -Q(i) - V(i)^2*B(i,i);
J2(ii,ii) = P(i)/V(i) + V(i)*G(i,i);
J3(ii,ii) = P(i) - V(i)^2*G(i,i);
J4(ii,ii) = -Q(i)/V(i) - V(i)*B(i,i);
end
end
end
% To remove rows & columns corresponding to PV buses from Jacobian.
r = find(type == 2);
l = length(r);
for i = l:-1:1
J2
J4
end
for i = l:-1:1
J3((r(i)-1),
J4((r(i)-1),
end
J = [J1 J2; J3 J4]; % Jacobian
X = inv(J)*M; % Correction Vector
dTh = X(1:n-1);
dV = X(n:length(X));
Th(2:n) = dTh + Th(2:n);
k = 1;
for i = 2:n
if type(i) == 3
V(i) = dV(k) + V(i);
k = k+1;
end
end
Iter = Iter + 1;
Tol = max(abs(X));
end
Iter % Number of Iterations took..
V;
del = 180/pi*Th;
[V del] % Bus Voltages and angles..
The bus5.m file follows:
function busdata = bus5() % Returns busdata.
busdata = [ 1 1 1.02 0 0.000 0.000;
2 3 1.00 0 -0.60 -0.30;
3 2 1.04 0 1.000 0.000;
4 3 1.00 0 -0.40 -0.10;
5 3 1.00 0 -0.60 -0.20 ];
lines.m follows:
function linedata = lines() % Returns linedata.
linedata = [1 2 0.00 0.40;
1 4 0.00 0.60;
1 5 0.00 0.20;
2 3 0.00 0.20;
2 4 0.00 0.40;
3 5 0.00 0.20;];
ybusppg.m follows:
function ybus = ybusppg(); % Returns ybus
linedata = lines(); % Calling "lines.m" for Line Data...
fb = linedata
tb = linedata
r = linedata
x = linedata
z = r + i*x; % Z matrix...
y = 1./z; % To get inverse of each element...
nbus = max(max(fb),max(tb)); % no. of buses...
nbranch = length(fb); % no. of branches...
ybus = zeros(nbus,nbus); % Initialise YBus...
% Formation of the Off Diagonal Elements...
for k=1:nbranch
ybus(fb(k),tb(k)) = -y(k);
ybus(tb(k),fb(k)) = ybus(fb(k),tb(k));
end
% Formation of Diagonal Elements....
for m=1:nbus
for n=1:nbranch
if fb
ybus(m,m) = ybus(m,m) + y
end
end
end
%ybus; % Bus Admittance Matrix
%zbus = inv(ybus); % Bus Impedance Matrix