配电网基于Matlab实现辐射状8节点配电网的状态估计 上传.zip

%Copyright(c) North Electrical Univ. All Rights Reserved %Description :Power System State Estimation Main Function %Author: S.D.H 2012.8 YY=Calculate_YY; zdata=Measurement; BusNode=max(max(zdata(:,4)),max(zdata(:,5))); type=zdata(:,2); z=zdata(:,3); FromNode=zdata(:,4); ToNode=zdata(:,5); Ri=diag(zdata(:,6)); V=ones(BusNode,1); Ang=zeros(BusNode,1); E=[Ang(2:end);V]; G=real(YY); B=imag(YY); vi=find(type==1); ppi=find(type==2); qi=find(type==3); pf=find(type==4); qf=find(type==5); nvi=length(vi); npi=length(ppi); nqi=length(qi); npf=length(pf); nqf=length(qf); iter=1;%迭代次数 tol=5;%收敛判据 while(tol>1e-6) % 计算h(x) h1=V(FromNode(vi),1); h2=zeros(npi,1); h3=zeros(nqi,1); h4=zeros(npf,1); h5=zeros(nqf,1); for i=1:npi m=FromNode(ppi(i)); for k=1:BusNode h2(i)=h2(i)+V(m)*V(k)*G(m,k)*cos(Ang(m)-Ang(k))+B(m,k)*sin(Ang(m)-Ang(k)); end end for i=1:nqi m=FromNode(qi(i)); for k=1:BusNode h3(i)=h3(i)+V(m)*V(k)*G(m,k)*sin(Ang(m)-Ang(k))-B(m,k)*cos(Ang(m)-Ang(k)); end end for i=1:npf m=FromNode(pf(i)); n=ToNode(pf(i)); h4(i)=-V(m)^2*G(m,n)+V(m)*V(n)*G(m,n)*cos(Ang(m)-Ang(n))+V(m)*V(n)*B(m,n)*sin(Ang(m)-Ang(n)); end for i=1:nqf m=FromNode(qf(i)); n=ToNode(qf(i)); h5(i)=V(m)^2*B(m,n)+V(m)*V(n)*G(m,n)*sin(Ang(m)-Ang(n))-V(m)*V(n)*B(m,n)*cos(Ang(m)-Ang(n)); end h=[h1;h2;h3;h4;h5]; %% %残差 r=z-h; %计算雅可比矩阵 %H11 电压量测对相角求导 H11=zeros(nvi,BusNode-1); %H12 电压量测对幅值求导 H12=zeros(nvi,BusNode); for k=1:nvi for n=1:BusNode if n==FromNode(vi(k)) H12(k,n)=1; end end end H21=zeros(npi,BusNode-1); for i=1:npi m=FromNode(ppi(i)); for k=1:(BusNode-1) if k+1==m for n=1:BusNode H21(i,k)=H21(i,k)+V(m)*V(n)*(-G(m,n)*sin(Ang(m)-Ang(n))+B(m,n)*cos(Ang(m)-Ang(n))); end H21(i,k)=H21(i,k)-V(m)^2*B(m,n); else H21(i,k)=V(m)*V(k+1)*(G(m,k+1)*sin(Ang(m)-Ang(k+1))-B(m,k+1)*cos(Ang(m)-Ang(k+1))); end end end H22=zeros(npi,BusNode); for i=1:npi m=FromNode(ppi(i)); for k=1:BusNode if k==m for n=1:BusNode H22(i,k)=H22(i,k)+V(n)*(G(m,n)*cos(Ang(m)-Ang(n))+B(m,n)*sin(Ang(m)-Ang(n))); end H22(i,k)=H22(i,k)+V(m)*G(m,m); else H22(i,k)=V(m)*(G(m,k)*cos(Ang(m)-Ang(k))+B(m,k)*sin(Ang(m)-Ang(k))); end end end H31 = zeros(nqi,BusNode-1); for i = 1:nqi m = FromNode(qi(i)); for k = 1:(BusNode-1) if k+1 == m for n = 1:BusNode H31(i,k) = H31(i,k) + V(m)* V(n)*(G(m,n)*cos(Ang(m)-Ang(n)) + B(m,n)*sin(Ang(m)-Ang(n))); end H31(i,k) = H31(i,k) - V(m)^2*G(m,m); else H31(i,k) = V(m)* V(k+1)*(-G(m,k+1)*cos(Ang(m)-Ang(k+1)) - B(m,k+1)*sin(Ang(m)-Ang(k+1))); end end end H32 = zeros(nqi,BusNode); for i = 1:nqi m = FromNode(qi(i)); for k = 1:(BusNode) if k == m for n = 1:BusNode H32(i,k) = H32(i,k) + V(n)*(G(m,n)*sin(Ang(m)-Ang(n)) - B(m,n)*cos(Ang(m)-Ang(n))); end H32(i,k) = H32(i,k) - V(m)*B(m,m); else H32(i,k) = V(m)*(G(m,k)*sin(Ang(m)-Ang(k)) - B(m,k)*cos(Ang(m)-Ang(k))); end end end H41=zeros(npf,BusNode-1); for i=1:npf m=FromNode(pf(i)); n=ToNode(pf(i)); for k=1:(BusNode-1) if k+1==m H41(i,k)=V(m)*V(n)*(-G(m,n)*sin(Ang(m)-Ang(n))+B(m,n)*cos(Ang(m)-Ang(n))); else if k+1==n H41(i,k)=-V(m)*V(n)*(-G(m,n)*sin(Ang(m)-Ang(n))+B(m,n)*cos(Ang(m)-Ang(n))); else H41(i,k)=0; end end end end H42=zeros(npf,BusNode); for i = 1:npf m = FromNode(pf(i)); n = ToNode(pf(i)); for k = 1:BusNode if k == m H42(i,k) = -V(n)*(-G(m,n)*cos(Ang(m)-Ang(n)) - B(m,n)*sin(Ang(m)-Ang(n))) - 2*G(m,n)*V(m); else if k == n H42(i,k) = -V(m)*(-G(m,n)*cos(Ang(m)-Ang(n)) - B(m,n)*sin(Ang(m)-Ang(n))); else H42(i,k) = 0; end end end end H51 = zeros(nqf,BusNode-1); for i = 1:nqf m = FromNode(qf(i)); n = ToNode(qf(i)); for k = 1:(BusNode-1) if k+1 == m H51(i,k) = -V(m)* V(n)*(-G(m,n)*cos(Ang(m)-Ang(n)) - B(m,n)*sin(Ang(m)-Ang(n))); else if k+1 == n H51(i,k) = V(m)* V(n)*(-G(m,n)*cos(Ang(m)-Ang(n)) - B(m,n)*sin(Ang(m)-Ang(n))); else H51(i,k) = 0; end end end end H52 = zeros(nqf,BusNode); for i = 1:nqf m = FromNode(qf(i)); n = ToNode(qf(i)); for k = 1:BusNode if k == m H52(i,k) = -V(n)*(-G(m,n)*sin(Ang(m)-Ang(n)) + B(m,n)*cos(Ang(m)-Ang(n))) - 2*V(m)*(-B(m,n)); else if k == n H52(i,k) = -V(m)*(-G(m,n)*sin(Ang(m)-Ang(n)) + B(m,n)*cos(Ang(m)-Ang(n))); else H52(i,k) = 0; end end end end %雅可比矩阵 H=[H11 H12;H21 H22;H31 H32;H41 H42;H51 H52]; Gm=H'*inv(Ri)*H; J=sum(inv(Ri)*r.^2); dE=inv(Gm)*(H'*inv(Ri)*r); E=E+dE; Ang(2:end)=E(1:BusNode-1); V=E(BusNode:end); iter=iter+1; tol=max(abs(dE)); end CvE=diag(inv(H'*inv(Ri)*H)); Angel=180/pi*Ang; E2=[V Angel]; disp('-------------state estimation--------------'); disp('-------------------------------------------'); disp('|Bus | V | Angel |'); disp('|No | pu | Degree |'); for m=1:BusNode fprintf('%4g',m);fprintf(' %8.4f ',V(m));fprintf(' %8.4f ',Angel(m));fprintf('\n'); end disp('------------------------------------------'); %%