MATLAB: Index in position 2 exceeds array bounds (must not exceed 13). , code to operate an IEEE 14/30 bus im looking into using WLS method. added a few details in order to plot the estimated vs real data. cannot work. please help!

Question

% Power System State Estimation using Weighted Least Square Method..

num = 30; % IEEE – 14 or IEEE – 30 bus system..(for IEEE-14 bus system replace 30 by 14)…

ybus = ybusppg(num); % Get YBus..

zdata = zdatas(num); % Get Measurement data..

bpq = bbusppg(num); % Get B data..

nbus = max(max(zdata(:,4)),max(zdata(:,5))); % Get number of buses..

type = zdata(:,2); % Type of measurement, Vi – 1, Pi – 2, Qi – 3, Pij – 4, Qij – 5, Iij – 6..

z = zdata(:,3); % Measuement values..

fbus = zdata(:,4); % From bus..

tbus = zdata(:,5); % To bus..

Ri = diag(zdata(:,6)); % Measurement Error..

V = ones(nbus,1); % Initialize the bus voltages..

del = zeros(nbus,1); % Initialize the bus angles..

E = [del(2:end); V]; % State Vector..

R = real(ybus);

UR = imag(ybus);

vi = find(type == 1); % Index of voltage magnitude measurements..

ppi = find(type == 2); % Index of real power injection measurements..

qi = find(type == 3); % Index of reactive power injection measurements..

pf = find(type == 4); % Index of real powerflow measurements..

qf = find(type == 5); % Index of reactive powerflow measurements..

nvi = length(vi); % Number of Voltage measurements..

npi = length(ppi); % Number of Real Power Injection measurements..

nqi = length(qi); % Number of Reactive Power Injection measurements..

npf = length(pf); % Number of Real Power Flow measurements..

nqf = length(qf); % Number of Reactive Power Flow measurements..

iter = 1;

tol = 5;

while(tol > 1e-4)

%Measurement Function, h

h1 = V(fbus(vi),1);

h2 = zeros(npi,1);

h3 = zeros(nqi,1);

h4 = zeros(npf,1);

h5 = zeros(nqf,1);

for i = 1:npi

m = fbus(ppi(i));

for k = 1:nbus

h2(i) = h2(i) + V(m)*V(k)*(R(m,k)*cos(del(m)-del(k)) + UR(m,k)*sin(del(m)-del(k)));

end

end

for i = 1:nqi

m = fbus(qi(i));

for k = 1:nbus

h3(i) = h3(i) + V(m)*V(k)*(R(m,k)*sin(del(m)-del(k)) – UR(m,k)*cos(del(m)-del(k)));

end

end

for i = 1:npf

m = fbus(pf(i));

n = tbus(pf(i));

h4(i) = -V(m)^2*R(m,n) – V(m)*V(n)*(-R(m,n)*cos(del(m)-del(n)) – UR(m,n)*sin(del(m)-del(n)));

end

for i = 1:nqf

m = fbus(qf(i));

n = tbus(qf(i));

h5(i) = -V(m)^2*(-UR(m,n)+bpq(m,n)) – V(m)*V(n)*(-R(m,n)*sin(del(m)-del(n)) + UR(m,n)*cos(del(m)-del(n)));

end

h = [h1; h2; h3; h4; h5];

% Residue..

r = z – h;

% Jacobian..

% H11 – Derivative of V with respect to angles.. All Zeros

H11 = zeros(nvi,nbus-1);

% H12 – Derivative of V with respect to V..

H12 = zeros(nvi,nbus);

for k = 1:nvi

for n = 1:nbus

if n == k

H12(k,n) = 1;

end

end

end

% H21 – Derivative of Real Power Injections with Angles..

H21 = zeros(npi,nbus-1);

for i = 1:npi

m = fbus(ppi(i));

for k = 1:(nbus-1)

if k+1 == m

for n = 1:nbus

H21(i,k) = H21(i,k) + V(m)* V(n)*(-R(m,n)*sin(del(m)-del(n)) + UR(m,n)*cos(del(m)-del(n)));

end

H21(i,k) = H21(i,k) – V(m)^2*UR(m,m);

else

H21(i,k) = V(m)* V(k+1)*(R(m,k+1)*sin(del(m)-del(k+1)) – UR(m,k+1)*cos(del(m)-del(k+1)));

end

end

end

% H22 – Derivative of Real Power Injections with V..

H22 = zeros(npi,nbus);

for i = 1:npi

m = fbus(ppi(i));

for k = 1:(nbus)

if k == m

for n = 1:nbus

H22(i,k) = H22(i,k) + V(n)*(R(m,n)*cos(del(m)-del(n)) + UR(m,n)*sin(del(m)-del(n)));

end

H22(i,k) = H22(i,k) + V(m)*R(m,m);

else

H22(i,k) = V(m)*(R(m,k)*cos(del(m)-del(k)) + UR(m,k)*sin(del(m)-del(k)));

end

end

end

% H31 – Derivative of Reactive Power Injections with Angles..

H31 = zeros(nqi,nbus-1);

for i = 1:nqi

m = fbus(qi(i));

for k = 1:(nbus-1)

if k+1 == m

for n = 1:nbus

H31(i,k) = H31(i,k) + V(m)* V(n)*(R(m,n)*cos(del(m)-del(n)) + UR(m,n)*sin(del(m)-del(n)));

end

H31(i,k) = H31(i,k) – V(m)^2*R(m,m);

else

H31(i,k) = V(m)* V(k+1)*(-R(m,k+1)*cos(del(m)-del(k+1)) – UR(m,k+1)*sin(del(m)-del(k+1)));

end

end

end

% H32 – Derivative of Reactive Power Injections with V..

H32 = zeros(nqi,nbus);

for i = 1:nqi

m = fbus(qi(i));

for k = 1:(nbus)

if k == m

for n = 1:nbus

H32(i,k) = H32(i,k) + V(n)*(R(m,n)*sin(del(m)-del(n)) – UR(m,n)*cos(del(m)-del(n)));

end

H32(i,k) = H32(i,k) – V(m)*UR(m,m);

else

H32(i,k) = V(m)*(R(m,k)*sin(del(m)-del(k)) – UR(m,k)*cos(del(m)-del(k)));

end

end

end

% H41 – Derivative of Real Power Flows with Angles..

H41 = zeros(npf,nbus-1);

for i = 1:npf

m = fbus(pf(i));

n = tbus(pf(i));

for k = 1:(nbus-1)

if k+1 == m

H41(i,k) = V(m)* V(n)*(-R(m,n)*sin(del(m)-del(n)) + UR(m,n)*cos(del(m)-del(n)));

else if k+1 == n

H41(i,k) = -V(m)* V(n)*(-R(m,n)*sin(del(m)-del(n)) + UR(m,n)*cos(del(m)-del(n)));

else

H41(i,k) = 0;

end

end

end

end

% H42 – Derivative of Real Power Flows with V..

H42 = zeros(npf,nbus);

for i = 1:npf

m = fbus(pf(i));

n = tbus(pf(i));

for k = 1:nbus

if k == m

H42(i,k) = -V(n)*(-R(m,n)*cos(del(m)-del(n)) – UR(m,n)*sin(del(m)-del(n))) – 2*R(m,n)*V(m);

else if k == n

H42(i,k) = -V(m)*(-R(m,n)*cos(del(m)-del(n)) – UR(m,n)*sin(del(m)-del(n)));

else

H42(i,k) = 0;

end

end

end

end

% H51 – Derivative of Reactive Power Flows with Angles..

H51 = zeros(nqf,nbus-1);

for i = 1:nqf

m = fbus(qf(i));

n = tbus(qf(i));

for k = 1:(nbus-1)

if k+1 == m

H51(i,k) = -V(m)* V(n)*(-R(m,n)*cos(del(m)-del(n)) – UR(m,n)*sin(del(m)-del(n)));

else if k+1 == n

H51(i,k) = V(m)* V(n)*(-R(m,n)*cos(del(m)-del(n)) – UR(m,n)*sin(del(m)-del(n)));

else

H51(i,k) = 0;

end

end

end

end

% H52 – Derivative of Reactive Power Flows with V..

H52 = zeros(nqf,nbus);

for i = 1:nqf

m = fbus(qf(i));

n = tbus(qf(i));

for k = 1:nbus

if k == m

H52(i,k) = -V(n)*(-R(m,n)*sin(del(m)-del(n)) + UR(m,n)*cos(del(m)-del(n))) – 2*V(m)*(-UR(m,n)+ bpq(m,n));

else if k == n

H52(i,k) = -V(m)*(-R(m,n)*sin(del(m)-del(n)) + UR(m,n)*cos(del(m)-del(n)));

else

H52(i,k) = 0;

end

end

end

end

% Measurement Jacobian, H..

H = [H11 H12; H21 H22; H31 H32; H41 H42; H51 H52];

% Gain Matrix, Gm..

Gm = H'*inv(Ri)*H;

%Objective Function..

J = sum(inv(Ri)*r.^2);

% State Vector..

dE = inv(Gm)*(H'*inv(Ri)*r);

E = E + dE;

del(2:end) = E(1:nbus-1);

V = E(nbus:end);

iter = iter + 1;

tol = max(abs(dE));

end

CvE = diag(inv(H'*inv(Ri)*H)); % Covariance matrix..

Del = 180/pi*del;

E2 = [V Del]; % Bus Voltages and angles..

disp('——– State Estimation ——————');

disp('————————–');

disp('| Bus | V | Angle | ');

disp('| No | pu | Degree | ');

disp('————————–');

for m = 1:n

fprintf('%4g', m); fprintf(' %8.4f', V(m)); fprintf(' %8.4f', Del(m)); fprintf('\n');

end

disp('———————————————');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Theta=mpc.bus(:,30); % actual voltage angles for comparison, to be estimated

z2=H2*Theta/180*pi; %this is measurement from bus meters, it shall be read from mpc=Gen-Load: z2=mpc.gen(:,2)-mpc.bus(:,3)

z=[z1;z2]; % overall measurement vector z

theta=inv(H'*H)*H'*z; % estimation method

theta=theta-theta(1); % making bus 1 angle zero, as a reference bus

theta=mod(theta/pi*180,180)-180; % convert to degree

theta(1)=theta(1)+180; % bus 1 to zero

plot(1:30,Theta,'b-',1:30,theta,'r–')

legend('Original', 'Estimated')

Answer 1

Hi, I tried to execute your code with some random input. The error in your code is due to the following line

tol = 5;
while(tol > 1e-4)

This condition always returns logical 1. So the loop executes continuously and the code after the while loop never executes. After changing the condition the code should work.