电力系统潮流计算代码

附录程序的主要代码：n=input('请输入节点数n=');na=input('请输入支路数na=');isb=input('请输入平衡节点母线号isb=');jd=input('请输入误差精度jd=');B1=input('请输入由支路参数形成的矩阵B1=');B2=input('请输入由节点参数形成的矩阵B2=');L=input('请输入由节点号及其对地阻抗形成的矩阵L=');nb=input('请输入P-Q节点数nb=');Y=zeros(n);Z=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);O=zeros(1,n);fori=1:naifB1(i,6)==0a=B1(i,1);b=B1(i,2);elsea=B1(i,2);b=B1(i,1);endY(a,b)=Y(a,b)-1./(B1(i,3)*B1(i,5));Z(a,b)=Z(a,b)-1./(B1(i,3));Y(b,a)=Y(a,b);Z(b,a)=Z(a,b);Y(b,b)=Y(b,b)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;Z(b,b)=Z(b,b)+1./(B1(i,3));Y(a,a)=Y(a,a)+1./(B1(i,3))+B1(i,4)./2;Z(a,a)=Z(a,a)+1./(B1(i,3));endG=real(Y);B=imag(Z);CI=imag(Y);fori=1:nS(i)=B2(i,1)-B2(i,2);CI(i,i)=CI(i,i)+B2(i,5);endP=real(S);Q=imag(S);fori=1:ne(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4);endfori=1:nifB2(i,6)==2V(i)=sqrt(e(i)^2+f(i)^2);O(i)=atan(f(i)./e(i));endendfori=2:nifi==nB(i,i)=1./B(i,i);elseIT1=i+1;forj1=IT1:nB(i,j1)=B(i,j1)./B(i,i);endB(i,i)=1./B(i,i);fork=i+1:nforj1=i+1:nB(k,j1)=B(k,j1)-B(k,i)*B(i,j1);endendendenda=0;b=0;fori=1:nifB2(i,6)==2a=a+1;k=0;forj1=1:nifB2(j1,6)==2k=k+1;A(a,k)=CI(i,j1);endendendendfori=1:nbifi==naA(i,i)=1./A(i,i);elsek=i+1;forj1=k:nbA(i,j1)=A(i,j1)./A(i,i);endA(i,i)=1./A(i,i);fork=i+1:nbforj1=i+1:nbA(k,j1)=A(k,j1)-A(k,i)*A(i,j1);endendendendNT2=1;NT1=0;kp=1;kq=1;K=1;NCT=0;NT3=1;whileNT2~=0|NT3~=0NT2=0;NT3=0;fori=1:nifi~=isbC(i)=0;fork=1:nC(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+CI(i,k)*sin(O(i)-O(k)));endCP1(i)=P(i)-V(i)*C(i);CP(i)=CP1(i)./V(i);NCT=abs(CP1(i));ifNCT=jdNT2=NT2+1;endendendNp(k)=NT2;ifNT2~=0fori=2:nCP(i)=B(i,i)*CP(i);ifi~=nIT1=i+1;fork=IT1:nCP(k)=CP(k)-B(k,i)*CP(i);endelseforLZ=3:iL=i+3-LZ;NC4=L-1;forMZ=2:NC4I=NC4+2-MZ;CP(I)=CP(I)-B(I,L)*CP(L);endendendendfori=2:nO(i)=O(i)-CP(i);endkq=1;L=0;fori=1:nifB2(i,6)==2C(i)=0;L=L+1;fork=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-CI(i,k)*cos(O(i)-O(k)));endDQ1(i)=Q(i)-V(i)*C(i);DQ(L)=DQ1(i)./V(i);NCT=abs(DQ1(i));ifNCT=jdNT3=NT3+1;endendendelsekp=0;ifkq~=0;L=0;fori=1:nifB2(i,6)==2C(i)=0;L=L+1;fork=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-CI(i,k)*cos(O(i)-O(k)));endDQ1(i)=Q(i)-V(i)*C(i);DQ(L)=DQ1(i)./V(i);NCT=abs(DQ1(i));endendendendNq(K)=NT3;ifNT3~=0L=0;fori=1:nbDQ(i)=A(i,i)*DQ(i);ifi==nbforLZ=2:iL=i+2-LZ;NC4=L-1;forMZ=1:NC4I=NC4+1-MZ;DQ(I)=DQ(I)-A(I,L)*DQ(L);endendelseIT1=i+1;fork=IT1:nbDQ(k)=DQ(k)-A(k,i)*DQ(i);endendendL=0;fori=1:nifB2(i,6)==2L=L+1;V(i)=V(i)-DQ(L);endendkp=1;K=K+1;elsekq=0;ifkp~=0K=K+1;endendfori=1:nDp(K-1,i)=V(i);endenddisp('迭代次数');disp(K);disp('每次没有达到精度要求的有功功率个数为');disp(Np);disp('每次没有达到精度要求的无功功率个数为');disp(Nq);fork=1:nE(k)=V(k)*cos(O(k))+V(k)*sin(O(k))*j;O(k)=O(k)*180./pi;enddisp('各节点的实际电压标幺值E(节点号从小到大排列):');disp(E);disp('各节点的电压大小V为(节点号从小到大排列):');disp(V);disp('各节点的电压相角O为(节点号从小到大排列):');disp(O);fora=1:nC(a)=0;forb=1:nC(a)=C(a)+conj(Y(a,b))*conj(E(b));endS(a)=E(a)*C(a);enddisp('各节点的功率S为(节点号从小到大排列):');disp(S);disp('各条支路的首端功率Si为(顺序同您输入B1时一样):');fori=1:naifB1(i,6)==0a=B1(i,1);b=B1(i,2);elsea=B1(i,2);b=B1(i,1);endSi(a,b)=E(a)*(conj(E(a))*conj(B1(i,4)./2)+(conj(E(a)*B1(i,5))-conj(E(b)))*conj(1./(B1(i,3)*B1(i,5))));disp(Si(a,b));enddisp('各条支路的末端功率Sj为(顺序同您输入B1时一样):');fori=1:naifB1(i,6)==0a=B1(i,1);b=B1(i,2);elsea=B1(i,2);b=B1(i,1);endSj(b,a)=E(b)*(conj(E(b))*conj(B1(i,4)./2)+(conj(E(b)./B1(i,5))-conj(E(a)))*conj(1./(B1(i,3)*B1(i,5))));disp(Sj(b,a));enddisp('各条支路的功率损耗DS为(顺序同您输入B1时一样)::');fori=1:naifB1(i,6)==0a=B1(i,1);b=B1(i,2);elsea=B1(i,2);b=B1(i,1);endDS(i)=Si(a,b)+Sj(b,a);disp(DS(i));endfori=1:KCs(i)=i;forj=1:nDp(K,j)=Dp(K-1,j);endenddisp('以下是每次迭代后各节点的电压值（如图所示）');plot(Cs,Dp),xlabel('迭代次数'),ylabel('电压'),title('电压迭代次数曲线');