programengine.m

状态估计原始程序,本程序可以进行状态估计的计算利用最小加权二乘法计算WLS

        Hqij(i,j)=0;
        if j==temp(1)
            Hqij(i,j)=v_bus(temp(1))*v_bus(temp(2))*(G(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2)))+B(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2))));
        end
        if j==temp(2)
            Hqij(i,j)=v_bus(temp(1))*v_bus(temp(2))*(-G(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2)))-B(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2))));
        end
        if j==temp(1)+nbus
            Hqij(i,j)=2*v_bus(temp(1))*B(temp(1),temp(2))+v_bus(temp(2))*(G(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2)))-B(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2))));
        end
        if j==temp(2)+nbus
            Hqij(i,j)=v_bus(temp(1))*(G(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2)))-B(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2))));
        end
    end
end
%Hqij

%formation of Hpji matrix corresponding to real power flow between bus j
%and i mismatches
%for i=1:length(mmpji)
 %   for j=1:nstate
  %      temp=line_tracker(:,i);
   %     trash=temp(1);
    %    temp(1)=temp(2);
     %   temp(2)=trash;
      %  Hpji(i,j)=0;
       % if j==temp(1)
        %    Hpji(i,j)=v_bus(temp(1))*v_bus(temp(2))*(-G(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2)))+B(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2))));
        %end
       % if j==temp(2)
           % Hpji(i,j)=v_bus(temp(1))*v_bus(temp(2))*(G(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2)))-B(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2))));
           %end
        %if j==temp(1)+nbus
         %   Hpji(i,j)=-2*v_bus(temp(1))*G(temp(1),temp(2))+v_bus(temp(2))*(G(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2)))+B(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2))));
         %end
        %if j==temp(2)+nbus
        %    Hpji(i,j)=v_bus(temp(1))*(G(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2)))+B(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2))));
        %end
        %end
        %end
%Hpji

%formation of Hqji matrix corresponding to reactive power flow between bus
%j and i mismatches
%for i=1:length(mmqji)
 %   for j=1:nstate
  %      temp=line_tracker(:,i);
   %     trash=temp(1);
   %     temp(1)=temp(2);
   %     temp(2)=trash;
   %     Hqji(i,j)=0;
   %     if j==temp(1)
   %         Hqji(i,j)=v_bus(temp(1))*v_bus(temp(2))*(G(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2)))+B(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2))));
   %     end
  %      if j==temp(2)
  %          Hqji(i,j)=v_bus(temp(1))*v_bus(temp(2))*(-G(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2)))-B(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2))));
  %      end
 %       if j==temp(1)+nbus
 %           Hqji(i,j)=2*v_bus(temp(1))*B(temp(1),temp(2))+v_bus(temp(2))*(G(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2)))-B(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2))));
 %       end
%        if j==temp(2)+nbus
%            Hqji(i,j)=v_bus(temp(1))*(G(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2)))-B(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2))));
%        end
%    end
%end
%Hqji

%birth of Jacobian
%H=[Hv; Hp; Hq; Hpij; Hqij; Hpji; Hqji]
H=[Hv; Hp; Hq; Hpij; Hqij];
hsize=size(H);
if count==0
   fid=fopen('output.txt','w');
   fprintf(fid,'Iteration number %d\n\n',count);
   fprintf(fid,'The Jacobian\n\n');
    for i=1:hsize(1)
    fprintf(fid,'%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',H(i,:));
end
fclose(fid);
end
if count>0
    fid=fopen('output.txt','a');
    fprintf(fid,'Iteration number %d\n\n',count);
    fprintf(fid,'The Jacobian\n\n');
    for i=1:hsize(1)
        fprintf(fid,'%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',H(i,:));
    end
    fclose(fid);
end
%determination of weight matrix
%v_wts
%p_wts
%q_wts
%pij_wts
%qij_wts
%pji_wts
%qji_wts
%W=eye(nmeas,nmeas);
W=zeros(nmeas,nmeas);
%W
for i=1:nbus
    %i
    for j=1:nmeas
     %   j
        W(i,j)=0;
        if j==i
            W(i,j)=v_wts(i);
        end
    end
    %W
end
i_temp=i;

for i=i_temp+1:i_temp+nbus
    %i
    for j=1:nmeas
        W(i,j)=0;
        if j==i
            W(i,j)=p_wts(i-i_temp);
        end
    end
    %W
end
i_temp=i;
for i=i_temp+1:i_temp+nbus
    %i
    for j=1:nmeas
        W(i,j)=0;
        if j==i
            W(i,j)=q_wts(i-i_temp);
        end
    end
    %W
end
i_temp=i;
for i=i_temp+1:i_temp+nlines
    for j=1:nmeas
        W(i,j)=0;
        if j==i
            W(i,j)=pij_wts(i-i_temp);
        end
    end
    %W
end
i_temp=i;
for i=i_temp+1:i_temp+nlines
    for j=1:nmeas
        W(i,j)=0;
        if j==i
            W(i,j)=qij_wts(i-i_temp);
        end
    end
    %W
end
%i_temp=i;
%for i=i_temp+1:i_temp+nlines
 %   for j=1:nmeas
  %      W(i,j)=0;
   %     if j==i
    %        W(i,j)=pji_wts(i-i_temp);
    %   end
    %end
   % W
   %end
%i_temp=i;
%for i=i_temp+1:i_temp+nlines
 %   for j=1:nmeas
  %      W(i,j)=0;
   %     if j==i
    %        W(i,j)=qji_wts(i-i_temp);
    %   end
    %end
    %W
    %end
 %   size(H)
  %  size(H')
   % size(W)

%determination of state mismatch vector
%W;
Gain=H'*W*H;
%Gain
Gain(1,1)=10000000;
gsize=size(Gain);
fid=fopen('output.txt','a');
fprintf(fid,'\n\nThe following is the Gain matrix.\n\n');
for i=1:gsize(1);
    fprintf(fid,'%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',Gain(i,:));
end
fclose(fid);
niag=inv(Gain);
smm=inv(Gain)*H'*W*mm;
%test=smm*1000
count=count+1;
%v_bus
%ang_bus
smm(1)=0; %to ensure that bus angle 1 remains as reference. no change in it
%smm
smma=smm(1:nbus,:);
smmv=smm(nbus+1:2*nbus,:);
%determination of tolerance 
plz_giveta=[];
for i=1:length(smm)
    if smm(i)>tolerance
        plz_giveta(i)=1;
    end
end
if length(plz_giveta)>0
   flag=1;
    v_bus=v_bus+smmv;
   ang_bus=ang_bus+smma;
end
%v_bus
%ang_bus
if length(plz_giveta)==0
    flag=0;
   % clc;
   fid=fopen('output.txt','a');
    fprintf(fid,'\n\nConvergence has occured in %d iterations.\n',count);
    fprintf(fid,'Estimated bus voltages.\n');
    fprintf(fid,'Bus number  Voltage\n');
    for i=1:nbus
        fprintf(fid,'%f %f\n',bus_number(i),v_bus(i));
    end
    fprintf(fid,'Estimated bus angles\n.');
    fprintf(fid,'Bus number  Angle\n');
    for i=1:nbus
        fprintf(fid,'%f %f\n',bus_number(i),ang_bus(i));
    end
    %Real power injection calculation
p_calc=[];
for i=1:nbus
    p_calc(i)=0;
    for j=1:nbus
        p_calc(i)=p_calc(i)+v_bus(i)*v_bus(j)*(G(i,j)*cos(ang_bus(i)-ang_bus(j))+B(i,j)*sin(ang_bus(i)-ang_bus(j)));
    end
end
fprintf(fid,'Estimated real power bus injections.\n');
    fprintf(fid,'Bus number  P\n');
    for i=1:nbus
        fprintf(fid,'%f %f\n',bus_number(i),p_calc(i));
    end

%reactive power injection calculations
q_calc=[];
for i=1:nbus
    q_calc(i)=0;
    for j=1:nbus
        q_calc(i)=q_calc(i)+v_bus(i)*v_bus(j)*(G(i,j)*sin(ang_bus(i)-ang_bus(j))-B(i,j)*cos(ang_bus(i)-ang_bus(j)));
    end
end
fprintf(fid,'Estimated reactive power bus injections.\n');
    fprintf(fid,'Bus number  Q\n');
    for i=1:nbus
        fprintf(fid,'%f %f\n',bus_number(i),q_calc(i));
    end

start_bus=x;
end_bus=y;

%i to j real flows calculation
pij_calc=[];
length(x);
for i=1:length(x)
    q=start_bus(i);
    r=end_bus(i);
    pij_calc(i)=-((v_bus(q))^2)*G(q,r)+v_bus(q)*v_bus(r)*(cos(ang_bus(q)-ang_bus(r))*G(q,r)+sin(ang_bus(q)-ang_bus(r))*B(q,r));
end
fprintf(fid,'Real line flows\n.');
fprintf(fid,'From_bus   To_bus  Real power flow\n');
for i=1:length(start_bus)
fprintf(fid,'%f %f  %f\n',start_bus(i),end_bus(i),pij_calc(i));

end

%i to j reactive flows calculation
qij_calc=[];
for i=1:length(x)
    q=start_bus(i);
    r=end_bus(i);
    qij_calc(i)=((v_bus(q))^2)*B(q,r)+v_bus(q)*v_bus(r)*(sin(ang_bus(q)-ang_bus(r))*G(q,r)-cos(ang_bus(q)-ang_bus(r))*B(q,r));
    
end
fprintf(fid,'Reactive line flows\n.');
fprintf(fid,'From_bus   To_bus  Reactive power flow\n');
for i=1:length(start_bus)
fprintf(fid,'%f %f  %f\n',start_bus(i),end_bus(i),qij_calc(i));

end
fclose(fid);
  break;
end
%if count==1
 %   fid=fopen('output.txt','w');
  %  fprintf(fid,'%6.2f %6.ff %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',smm');
   %fclose(fid);
   %end
%if count>1
 %  fid=fopen('output.txt','a');
  %fprintf(fid,'%6.2f %6.ff %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',smm');
 %fclose(fid);
 %end
%count
end