Continue to Site

Eng-Tips is the largest engineering community on the Internet

Intelligent Work Forums for Engineering Professionals

  • Congratulations Danlap on being selected by the Eng-Tips community for having the most helpful posts in the forums last week. Way to Go!

Inner Matrix Dimensions Must Agree Trouble 1

Status
Not open for further replies.

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(n) == m | tb(n) == m
ybus(m,m) = ybus(m,m) + y(n);
end
end
end
%ybus; % Bus Admittance Matrix
%zbus = inv(ybus); % Bus Impedance Matrix
 
Replies continue below

Recommended for you

I reckon M has dimensions (2n-1) x 1
As for J...
Well J1, J2, J3 and J4 are initially defined as (n-1) x (n-1)

As one of the entries in busdata is of type 2, 1 row is removed from J2, 1 column is removed from J3, and 1 row and 1 column is removed from J4.

J1: (n-1)x(n-1)
J2: (n-2)x(n-1)
J3: (n-1)x(n-2)
J4: (n-2)x(n-2)


By my reckconing, the program should fail at the line J = [J1 J2; J3 J4];

M



--
Dr Michael F Platten
 
Thanks, MikeyP...that's exactly what the problem was. Does anyone have any ideas how I could go about to remedy this?
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor