电力系统潮流计算-仿真结果自动保存到excel文件中-源码

clear %先求节点导纳矩阵 N1=30; %节点数 L1=43; %支路数 %N1=5;L1=7; B1=xlsread('XIANLU2.xlsx','B2:F44'); %B1=xlsread('XIANLU.xlsx','B2:F44'); G=zeros(N1); B=zeros(N1); %每行存储一条支路 %第一列存储支路的一个端点I,变压器所在端加个负号，励磁支路为0 %第二列存储另一个端点J，变压器所在端加个负号，励磁支路为0 %第三列存储支路的电阻R %第四列存储支路电抗X %第五列存储支路对地电纳或变压器变比k,k在首端为正 for m1=1:L1 I=B1(m1,1);J=B1(m1,2);R=B1(m1,3);X=B1(m1,4);k=B1(m1,5); if I*J==0 %励磁支路？ if I==0 G(J,J)=G(J,J)+R; B(J,J)=B(J,J)+X; end if J==0 G(I,I)=G(I,I)+R; B(I,I)=B(I,I)+X; end else if I*J>0 %线路支路？ B(I,I)= B(I,I)+ k; B(J,J)= B(J,J)+ k; k=1; else if I<0 t=I; I=J; J=t; end J=abs(J); if k<0 k=-1/k; end end G(I,J)=-(1/k)*R/(R^2+X^2); G(J,I)=G(I,J); B(I,J)=(1/k)*X/(R^2+X^2); B(J,I)=B(I,J); G(I,I)=G(I,I)+R/(R^2+X^2); B(I,I)=B(I,I)-X/(R^2+X^2); G(J,J)=G(J,J)+1/(k^2)*R/(R^2+X^2); B(J,J)=B(J,J)-1/(k^2)*X/(R^2+X^2); end end Y=G+B*1i; %求差. b=xlsread('JIEDIAN2.xlsx','C2:F31'); %b=xlsread('JIEDIAN.xlsx','C2:F31'); precision=1; %误差 t=0; %迭代次数 pq=24; %pq节点数 %pq=4; U=b(pq+2:N1,1); %开始牛拉法 while precision>0.00001 P1=zeros(pq,1); Q=zeros(pq,1); P2=zeros(N1-pq-1,1); U2=zeros(N1-pq-1,1); deltpqu=zeros(2*(N1-1),1); deltu2=zeros(N1-pq-1,1); for i=1:pq for j=1:N1 P1(i,1)=P1(i,1)+b(i+1,1)*(G(i+1,j)*b(j,1)-B(i+1,j)*b(j,2))+b(i+1,2)*(G(i+1,j)*b(j,2)+B(i+1,j)*b(j,1)); Q(i,1) =Q(i,1) +b(i+1,2)*(G(i+1,j)*b(j,1)-B(i+1,j)*b(j,2))-b(i+1,1)*(G(i+1,j)*b(j,2)+B(i+1,j)*b(j,1)); end end for i=1:N1-pq-1 for j=1:N1 P2(i,1)=P2(i,1)+b(pq+i+1,1)*(G(pq+i+1,j)*b(j,1)-B(pq+i+1,j)*b(j,2))+b(pq+i+1,2)*(G(pq+i+1,j)*b(j,2)+B(pq+i+1,j)*b(j,1)); end U2(i,1)=b(pq+i+1,1)^2+b(pq+i+1,2)^2; end deltp1=b(2:pq+1,3)-P1(1:pq,1); deltq =b(2:pq+1,4)-Q(1:pq,1); deltp2=b(pq+2:N1,3)-P2(1:N1-pq-1,1); for i=1:N1-pq-1 deltu2(i,1)=U(i,1)^2-U2(i,1); end deltu2; for i=1:pq deltpqu(2*i-1,1)=deltp1(i); deltpqu(2*i,1)=deltq(i); end for i=1:N1-pq-1 deltpqu(2*pq+2*i-1,1)=deltp2(i); deltpqu(2*pq+2*i,1)=deltu2(i); end deltpqu; %delta %求雅克比 aii=zeros(N1);bii=zeros(N1); H1=zeros(N1);N11=zeros(N1); for i=1:N1 for j=1:N1 aii(i,1)=aii(i,1)+G(i,j)*b(j,1)-B(i,j)*b(j,2); bii(i,1)=bii(i,1)+G(i,j)*b(j,2)+B(i,j)*b(j,1); end end %H、N for i=1:N1 for j=1:N1 if i==j H1(i,i)=-B(i,i)*b(i,1)+G(i,i)*b(i,2)+bii(i,1); N11(i,i)= G(i,i)*b(i,1)+B(i,i)*b(i,2)+aii(i,1); else H1(i,j)=-B(i,j)*b(i,1)+G(i,j)*b(i,2); N11(i,j)= G(i,j)*b(i,1)+B(i,j)*b(i,2); end end end H=H1(2:N1,2:N1);N=N11(2:N1,2:N1); %J,L,R,S J1=zeros(N1);L1=zeros(N1); R1=zeros(N1);S1=zeros(N1); for i=1:N1 for j=1:N1 if i==j J1(i,i)=-G(i,i)*b(i,1)-B(i,i)*b(i,2)+aii(i); L1(i,i)=-B(i,i)*b(i,1)+G(i,i)*b(i,2)-bii(i); else J1(i,j)=-N11(i,j);L1(i,j)=H1(i,j); end end end J=J1(2:pq+1,2:N1);L=L1(2:pq+1,2:N1); for i=1:N1 for j=1:N1 if i==j R1(i,i)=2*b(i,2);S1(i,i)=2*b(i,1); else R1(i,j)=0;S1(i,j)=0; end end end R=R1((pq+2):N1,2:N1); S=S1((pq+2):N1,2:N1); %雅克比 Jacobi=zeros(2*(N1-1)); for i=1:N1-1 for j=1:N1-1 Jacobi(2*i-1,2*j-1)=H(i,j); Jacobi(2*i-1,2*j) =N(i,j); end end for i=1:pq for j=1:N1-1 Jacobi(2*i,2*j-1)=J(i,j); Jacobi(2*i,2*j) =L(i,j); end end for i=1:N1-pq-1 for j=1:N1-1 Jacobi(2*pq+2*i,2*j-1)=R(i,j); Jacobi(2*pq+2*i,2*j) =S(i,j); end end Jacobi; %修正 Correction=Jacobi^-1*deltpqu; for i=1:N1-1 b(i+1,1)=b(i+1,1)+Correction(2*i); b(i+1,2)=b(i+1,2)+Correction(2*i-1); end precision=max(abs(deltpqu)); %precision=max(abs(Correction)); t=t+1; b; end %平衡节点功率 pp=0; for j=1:N1 pp=pp+conj(Y(1,j))*(b(j,1)-b(j,2)*1i); end b(1,3)=real(b(1,1)*pp); b(1,4)=imag(b(1,1)*pp); %各线路功率 S=zeros(N1); y=zeros(N1); y0=zeros(N1); for i=1:N1 for j=1:N1 if i~=j y(i,j)=-Y(i,j); end end end y; for m1=1:L1 I=B1(m1,1);J=B1(m1,2);R=B1(m1,3);X=B1(m1,4);k=B1(m1,5); if I*J>0 y0(I,J)=k*1i; y0(J,I)=k*1i; end if I*J<0 if I<0 t=I; I=J; J=t; end J=abs(J); if k<0 k=-1/k; end y0(I,J)=(R/(R^2+X^2)-X/(R^2+X^2)*1i)*(k-1)/k; y0(J,I)=(R/(R^2+X^2)-X/(R^2+X^2)*1i)*(1-k)/k^2; end end y0; SS=0; for m=1:N1 for n=1:N1 if (m~=n)&&y(m,n)~=0 S(m,n)=(b(m,1)+b(m,2)*1i)*((b(m,1)-b(m,2)*1i)*conj(y0(m,n))+((b(m,1)-b(m,2)*1i)-(b(n,1)-b(n,2)*1i))*conj(y(m,n))); SS=SS+S(m,n); else S(m,n)=0; end end end SB=zeros(N1,1); for m=1:N1 for n=1:N1 SB(m)=SB(m)+S(m,n); if m>pq b(m,3)=real(SB(m)); b(m,4)=imag(SB(m)); end end end %线路上损耗的功率 DS=zeros(N1); for i=1:N1 for j=1:N1 DS(i,j)=S(i,j)+S(j,i); end end siu=zeros(43,1); for m1=1:L1 I=B1(m1,1);J=B1(m1,2); siu(m1,1)=DS(I,J); end t=t-1 precision b; S; siu; SS; xlswrite('eifi.xlsx',b,'D2:G31'); xlswrite('XIANLU2.xlsx',real(siu(1:43)),'G2:G44'); xlswrite('XIANLU2.xlsx',imag(siu(1:43)),'H2:H44'); winopen('XIANLU2.xlsx'); winopen('eifi.xlsx');