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!

1 view (last 30 days)
Ritaj Alsarawi
Ritaj Alsarawi on 9 May 2020
Commented: Ritaj Alsarawi on 16 May 2020
% 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')

Sign in to answer this question.

Accepted Answer

Guru Mohanty
Guru Mohanty on 13 May 2020
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.

More Answers (0)

Categories

Find more on Matrix Operations 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!