matlab潮流计算程序14节点

JD=14;ZL=20;WC=0.01;ZP=[1,2,0.01335,0.04211,0.0,0.0;1,3,0.0,0.20912,0.0,0.0;1,4,0.0,0.55618,0.0,0.0;1,10,0.05811,0.17632,0.0,0.0340;1,11,0.06701,0.17103,0.0,0.0128;2,10,0.05695,0.17388,0.0,0.0346;2,12,0.0,0.25202,0.0,0.0;2,14,0.05403,0.22304,0.0,0.0492;3,4,0.0,0.11001,0.0,0.0;3,13,0.0,0.17615,0.0,0.0;4,5,0.03181,0.08450,0.0,0.0;4,9,0.12711,0.27038,0.0,0.0;5,6,0.08205,0.19207,0.0,0.0;6,12,0.09498,0.19890,0.0,0.0;7,8,0.22092,0.19988,0.0,0.0;7,12,0.12291,0.25581,0.0,0.0;8,9,0.17093,0.34802,0.0,0.0;8,12,0.06615,0.13027,0.0,0.0;10,11,0.04699,0.19797,0.0,0.0438;10,14,0.01938,0.05917,0.0,0.0528];JP=[1,1,-0.478,0.039,1.0,0.0;2,1,-0.076,-0.016,1.0,0.0;3,1,0.0,0.0,1.0,0.0;4,1,-0.295,-0.166,1.0,0.0;5,1,-0.09,-0.058,1.0,0.0;6,1,-0.035,-0.018,1.0,0.0;7,1,-0.061,-0.016,1.0,0.0;8,1,-0.135,-0.058,1.0,0.0;9,1,-0.149,-0.05,1.0,0.0;10,2,0.183,0.0,1.045,0.0;11,2,-0.942,0.0,1.010,0.0;12,2,-0.112,0.047,1.70,0.0;13,2,0.0,0.174,1.90,0.0;14,0,0.0,0.0,1.060,0.0];[N1,N2]=size(ZP);p=max(ZP(:,1));q=max(ZP(:,2));max1=max(p,q);Y=zeros(max1);forn=1:N1i1=ZP(n,1);j1=ZP(n,2);r=ZP(n,3);x=ZP(n,4);g=ZP(n,5);b=ZP(n,6);m=r*r+x*x;ifi1==j1Y(i1,i1)=Y(i1,i1)+r/m-j*x/m+j*b;elseY(i1,j1)=Y(i1,j1)-r/m+j*x/m;Y(j1,i1)=Y(i1,j1);Y(i1,i1)=Y(i1,i1)+r/m-j*x/m+j*b;Y(j1,j1)=Y(j1,j1)+r/m-j*x/m+j*b;endend%导纳矩阵生成V=zeros(1,JD);XJ=zeros(1,JD);P=zeros(1,JD);Q=zeros(1,JD);V2=zeros(1,JD);G=real(Y);B=imag(Y);fori=1:JDV(i)=JP(i,5);XJ(i)=JP(i,6);P(i)=JP(i,3);Q(i)=JP(i,4);V2(i)=V(i)*V(i);end%赋值DD=0;PD=1;fori=1:JDe(i)=V(i)*cos(XJ(i));f(i)=V(i)*sin(XJ(i));endwhilePD~=0PD=0;DD=DD+1;n=JD-1;H=zeros(n);N=zeros(n);J=zeros(n);L=zeros(n);R=zeros(n);S=zeros(n);fori=1:nfork=1:nif(i~=k)H(i,k)=-B(i,k)*e(i)+G(i,k)*f(i);N(i,k)=G(i,k)*e(i)+B(i,k)*f(i);J(i,k)=-N(i,k);L(i,k)=H(i,k);R(i,k)=0;S(i,k)=0;elsea1=0;b1=0;fori1=1:na1=a1+G(i,i1)*e(i1)-B(i,i1)*f(i1);b1=1+G(i,i1)*f(i1)-B(i,i1)*e(i1);endH(i,k)=-B(i,k)*e(i)+G(i,k)*f(i)+b1;N(i,k)=G(i,k)*e(i)+B(i,k)*f(i)+a1;J(i,k)=-G(i,k)*e(i)-B(i,k)*f(i)+a1;L(i,k)=-B(i,k)*e(i)+G(i,k)*f(i)-b1;R(i,k)=2*f(i);S(i,k)=2*e(i);endendendJ2n=zeros(2*n);fori=1:nif(JP(i,2)==1)fork=1:nJ2n(2*i-1,2*k-1)=H(i,k);J2n(2*i-1,2*k)=N(i,k);J2n(2*i,2*k-1)=J(i,k);J2n(2*i,2*k)=L(i,k);endelsefork1=1:nJ2n(2*i-1,2*k1-1)=H(i,k1);J2n(2*i-1,2*k1)=N(i,k1);J2n(2*i,2*k1-1)=R(i,k1);J2n(2*i,2*k1)=S(i,k1);endendendfori=1:nXP(i)=0;XQ(i)=0;XV2(i)=e(i)*e(i)+f(i)*f(i);fork=1:nXP(i)=XP(i)+e(i)*(G(i,k)*e(k)-B(i,k)*f(k))+f(i)*(G(i,k)*f(k)+B(i,k)*e(k));XQ(i)=XQ(i)+f(i)*(G(i,k)*e(k)-B(i,k)*f(k))-e(i)*(G(i,k)*f(k)+B(i,k)*e(k));endDP(i)=P(i)-XP(i);DQ(i)=Q(i)-XQ(i);DV2(i)=V2(i)-XV2(i);endPQU=zeros(2*n,1);fori=1:nPQU(2*i-1,1)=DP(i);if(JP(i,2)==1)PQU(2*i,1)=DQ(i);elsePQU(2*i,1)=DV2(i);endendmax2=max(abs(PQU));disp(max2);ifmax2WCPD=1;endFE=inv(J2n)*PQU;fori=1:ne(i)=e(i)+FE(2*i);f(i)=f(i)+FE(2*i-1);endenddisp(PQU);disp(DD);