PQ.rar_PQ分解法_PQ分解法matlab_PQ分解法潮流_matlab pq_电力系统PQ

function[tab2,tab4,tab6,cycle]=PQ(filename) base=100; b=textread('case30.txt'); %读文件 [m1,n1]=size(b); %求节点数，支路数 n=0; %节点数 m=0; %支路数 for k=1:200 %最大节点数200 if b(k,1)~=-999 n=n+1; else break; end end c=k+1; %指向第一个-999的后一行数 for k=c:1000 %最大支路数 if b(k,1)~=-999 m=m+1; else break; end end p=zeros(n,n1+1);%节点矩阵 z=zeros(m,n1+2);%支路矩阵 p(1:n,1:n1)=b(1:n,1:n1); z(1:m,1:n1)=b(n+2:m1-1,1:n1); %筛选pq，pv，平衡节点，并排序 x1=0; for i=1:n if p(i,5)==1 x1=x1+1; p(i,n1+1)=x1; end end x2=x1; for i=1:n if p(i,5)==2 x2=x2+1; p(i,n1+1)=x2; end end x3=x2; for i=1:n if p(i,5)==3 x3=x3+1; p(i,n1+1)=x3; end end for i=1:m for j=1:n if z(i,1)==p(j,1) z(i,22)=p(j,22);break; end end for j=1:n if z(i,2)==p(j,1) z(i,23)=p(j,22);break; end end end for k=1:(n-1) for j=1:(n-k) if p(j+1,22)<p(j,22) bj=p(j+1,:); p(j+1,:)=p(j,:); p(j,:)=bj; end end end for i=1:m for j=1:n if z(i,1)==p(j,1) z(i,22)=p(j,22);break; end end for j=1:n if z(i,2)==p(j,1) z(i,23)=p(j,22);break; end end end %%%%%%%%%%%%%形成节点导纳矩阵B q1=0;%pq节点个数 v1=0;%pv节点个数 for k=1:n if p(k,5)==2 v1=v1+1; %这样只要求pv节点类型有值就行 end end q1=n-v1-1; B=zeros(n,n); B1=zeros(n-1,n-1); B2=zeros(q1,q1); for i=1:m z(i,6)=z(i,6);%电阻 z(i,7)=3*z(i,7);%电抗 z(i,8)=z(i,8);%电纳 end for i=1:m j=z(i,22); k=z(i,23); if z(i,5)==1 %%%%%支路为线路 B(j,j)=B(j,j)+1i*z(i,8)/2+1/(z(i,6)+1i*z(i,7)); B(k,k)=B(k,k)+1i*z(i,8)/2+1/(z(i,6)+1i*z(i,7)); B(j,k)=-1/(z(i,6)+1i*z(i,7)); B(k,j)=-1/(z(i,6)+1i*z(i,7)); else if z(i,5)==2 %%%%%%%%%%支路为变压器 变压器的形式：i k：1 Zt j 的形式 B(k,k)=B(k,k)+1i*z(i,8); B(k,k)=B(k,k)+(z(i,14)-1)/(z(i,14)*(z(i,6)+1i*z(i,7)))+1/(z(i,14)*(z(i,6)+1i*z(i,7))); B(j,j)=B(j,j)+(1-z(i,14))/((z(i,14))^2*(z(i,6)+1i*z(i,7)))+1/(z(i,14)*(z(i,6)+1i*z(i,7))); B(j,k)=-1/(z(i,14)*(z(i,6)+1i*z(i,7))); B(k,j)=-1/(z(i,14)*(z(i,6)+1i*z(i,7))); end end end for k=1:n B(k,k)=B(k,k)+1i*p(k,17); %加上对地电容 end Y=B; G=real(B); B=imag(B); B2=B(1:q1,1:q1); B1=B(1:n-1,1:n-1); %%求初始差值 P2=zeros(n-1,1);%真实的P Q2=zeros(q1,1); Q3=zeros(n-1,1); P21=zeros(n-1,1);%P差值 Q21=zeros(q1,1); U=zeros(n,1); U2=zeros(q1,1); A=zeros(n,1); U21=zeros(q1,1); A21=zeros(n-1,1); U=p(1:n,13); U2=U(1:q1,1); A=p(1:n,7)./(180/pi); P2=(p(1:n-1,10)-p(1:n-1,8))./base; Q3=(p(:,11)-p(:,9))./base; Q2=Q3(1:q1,1); maxError=0; Node=0; for cycle=1:2000 P21=P2(1:n-1,1); Q21=Q2; for i=1:n-1 for j=1:n P21(i,1)=P21(i,1)-U(i,1)*U(j,1)*(G(i,j)*cos(A(i,1)-A(j,1))+B(i,j)*sin(A(i,1)-A(j,1))); end if P21(i,1)>maxError maxError=P21(i,1); Node=i; end end P21=P21./U(1:n-1,1); A21=-1*inv(B1)*P21; A21=A21./U(1:n-1,1); A(1:n-1,1)=A(1:n-1,1)+A21; for i=1:n-1 for j=1:n if(i<=q1) Q21(i,1)=Q21(i,1)-U(i,1)*U(j,1)*(G(i,j)*sin(A(i,1)-A(j,1))-B(i,j)*cos(A(i,1)-A(j,1))); end end if (i<=q1)&Q21(i,1)>maxError maxError=Q21(i,1); Node=i; end end Q21=Q21./U2; U21=-1*inv(B2)*Q21; U2=U2+U21; for k=1:q1 U(k,1)=U2(k,1); end if(abs(all(P21(:))))<1.0e-5&(abs(all(Q21(:))))<1.0e-5 disp('success'); break; end end P=zeros(n,1); P(1:n-1,1)=P2.*base; Q=zeros(n,1); Q(1:q1,1)=Q2.*base; for i=q1+1:n for j=1:n Q(i,1)=Q(i,1)+base*U(i,1)*U(j,1)*(G(i,j)*sin( (A(i,1)-A(j,1)))-B(i,j)*cos( (A(i,1)-A(j,1)))); if(i==n) P(i,1)=P(i,1)+base*U(i,1)*U(j,1)*(G(i,j)*cos( (A(i,1)-A(j,1)))+B(i,j)*sin( (A(i,1)-A(j,1)))); end end end for k=1:(n-1) for j=1:(n-k) if p(j+1,1)<p(j,1) bj=p(j+1,:); p(j+1,:)=p(j,:); p(j,:)=bj; bj=U(j+1,:); U(j+1,:)=U(j,:); U(j,:)=bj; bj=A(j+1,:); A(j+1,:)=A(j,:); A(j,:)=bj; bj=P(j+1,:); P(j+1,:)=P(j,:); P(j,:)=bj; bj=Q(j+1,:); Q(j+1,:)=Q(j,:); Q(j,:)=bj; end end end U0=U.*(cos(A)+sin(A).*1i); A=A*180/pi; S1=zeros(m,1); S2=zeros(m,1); S21=zeros(m,1); for k=1:m i=z(k,1); j=z(k,2); if z(k,5)==1 S1(k,1)=U(i,1)^2*(-1i*z(k,8)/2)+U0(i,1)*(conj(U0(i,1))-conj(U0(j,1)))*conj(1/(z(k,6)+1i*z(k,7))); S2(k,1)=U(j,1)^2*(-1i*z(k,8)/2)+U0(j,1)*(conj(U0(j,1))-conj(U0(i,1)))*conj(1/(z(k,6)+1i*z(k,7))); S21(k,1)=S1(k,1)+S2(k,1); else if z(k,5)==2 S1(k,1)=U(i,1)^2*conj(1i*z(k,8)/2+(1-z(k,14))/((z(k,14))^2*(z(k,6)+1i*z(k,7))))+U0(i,1)*(conj(U0(i,1))-conj(U0(j,1)))*conj(1/(z(k,14)*(z(k,6)+1i*z(k,7)))); S2(k,1)=U(j,1)^2*conj(1i*z(k,8)/2+(z(k,14)-1)/(z(k,14)*(z(k,6)+1i*z(k,7))))+U0(j,1)*(conj(U0(j,1))-conj(U0(i,1)))*conj(1/(z(k,14)*((z(k,6)+1i*z(k,7))))); S21(k,1)=S1(k,1)+S2(k,1); end end end S21=S21*100; Ps=real(S21); Qs=imag(S21); tab1=[{'节点'},{'节点电压'},{'节点相角'}]; tab2=[ p(:,1) U A ]; tab3=[{'最大误差'},{'最大误差节点'}]; tab4=[ maxError Node ]; tab5=[{'始节点'},{'终结点'},{'有功损耗'},{'无功损耗'}]; tab6=[ z(:,1) z(:,2) Ps Qs]; xlswrite('PQ.xls',tab1,1,'A1'); xlswrite('PQ.xls',tab2,1,'A2'); xlswrite('PQ.xls',tab3,1,'D1'); xlswrite('PQ.xls',tab4,1,'D2'); xlswrite('PQ.xls',tab5,1,'F1'); xlswrite('PQ.xls',tab6,1,'F2');