1、程序代碼如下：111111.%讀入數(shù)據(jù)clcclearfilename=123.txt;a=textread(filename)n=a(1,1);pinghengjd=a(1,2);phjddianya=a(1,3);jingdu=a(1,4);b=zeros(1,9);j1=0;m1,n1=size(a);for i1=1:m1if a(i1,1)=0j1=j1+1;b(j1)=i1;endendb;%矩陣分塊a1=a(b(1)+1:b(2)-b(1)+1,1:n1);a2=a(b(2)+1:b(3)-1,1:n1);a3=a(b(3)+1:b(4)-1,1:n1);a4=a(b(4)+1:

2、b(5)-1,1:n1);a5=a(b(5)+1:b(6)-1,1:n1);%設(shè)置初值vcz=1;dcz=0;kmax=20;k1=0;%求節(jié)點(diǎn)導(dǎo)納矩陣a11=zeros(4,6);for i0=1:3for j0=1:6a11(i0,j0)=a1(i0,j0);a11(4,j0)=a2(1,j0);endendall;linei=a11(1:4,2);linej=a11(1:4,3);liner=a11(1:4,4);linex=a11(1:4,5);lineb=a11(1:4,6);branchi=0;branchj=0;branchb=0;G=zeros(4,4);B=zeros(4,4

3、);for k=1:4i2=linei(k,1);j2=linej(k,1);r=liner(k,1);x=linex(k,1);b=0;GIJ=r/(r*r+x*x);BIJ=-x/(r*r+x*x);if k=4 & lineb(k)=0k0=lineb(k);G(i2,j2)=-GIJ/k0;G(j2,i2)=G(i2,j2);B(i2,j2)=-BIJ/k0;B(j2,i2)=B(i2,j2);G(i2,i2)=G(i2,i2)+GIJ/k0/k0;B(i2,i2)=B(i2,i2)+BIJ/k0/k0; elseG(j2,i2)=-GIJ;G(i2,j2)=G(j2,i2);B(j2

4、,i2)=-BIJ;B(i2,j2)=B(j2,i2);G(i2,i2)=G(i2,i2)+GIJ;b=lineb(k);B(i2,i2)=B(i2,i2)+BIJ+b;endG(j2,j2)=G(j2,j2)+GIJ;B(j2,j2)=B(j2,j2)+BIJ+b;endG;B;B=B.*i;Yf=G+BY=abs(Yf);alf=angle(Yf);%賦Jacobian矩陣參數(shù)P=zeros(n,1);Q=zeros(n,1);Pd=zeros(1,n);Qd=zeros(1,n);dP=zeros(1,n);dQ=zeros(1,n);PG=a4(:,3);PD=a4(:,5);QG=a

5、4(:,4);QD=a4(:,6);i8=a4(:,2);for j8=1:length(i8)P(i8(j8)=PG(i8(j8)-PD(i8(j8);Q(i8(j8)=QG(i8(j8)-QD(i8(j8);enddelt=zeros(n,1);V=ones(n,1);V(3)=1.10;V(4)=1.05;ddelt=zeros(n,1);dV=zeros(n,1);A=zeros(2*n,2*n);B=zeros(2*n,1);Jacobian=Jaco(V,delt,n,Y,alf)%求取矩陣功率for j5=1:kmaxdisp(第int2str(j5)次計(jì)算結(jié)果)if k=kma

6、xbreakend for i10=1:4 Pd(i10)=0; Qd(i10)=0;for j10=1:nPd(i10)=Pd(i10)+V(i10)*Y(i10,j10)*V(j10)*cos(d elt(i10)-delt(j10)-alf(i10,j10);Qd(i10)=Qd(i10)+V(i10)*Y(i10,j10)*V(j10)*sin(d elt(i10)-delt(j10)-alf(i10,j10);endendfor i4=1:3 dP(i4)=P(i4)-Pd(i4);end for j4=1:2 dQ(j4)=Q(j4)-Qd(j4);endA=Jaco(V,delt

7、,n,Y,alf) for i14=1:n B(i14*2-1)=-dP(i14);B(i14*2)=-dQ(i14);end if max(abs(B)jingdu X=AB;for i16=1:n ddelt(i16)=X(2*i16-1);dV(i16)=X(2*i16)*V(i16);endV=V+dVdelt=delt+ddeltelse break enddisp()end%流氓算法%for ii=1:2%V(ii)=V(ii)+dV(ii);%end%V222222.function A=Jaco(V,delt,n,Y,alf)%計(jì)算Jacobian矩陣for i7=1:nHd1

8、(i7)=0;Jd1(i7)=0;for j7=1:nHd1(i7)=Hd1(i7)+V(i7)*Y(i7,j7)*V(j7)*sin(delt(i7)-delt (j7)-alf(i7,j7);Jd1(i7)=Jd1(i7)+V(i7)*Y(i7,j7)*V(j7)*cos(delt(i7)-delt (j7)-alf(i7,j7);end endfor i6=1:nfor j6=1:n if i6=j6H(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*sin(delt(i6)-de lt(j6)-alf(i6,j6);N(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)

9、*cos(delt(i6)-de lt(j6)-alf(i6,j6);J(i6,j6)=-N(i6,j6);L(i6,j6)=H(i6,j6);elseH(i6,i6)=Hd1(i6)-V(i6)*Y(i6,i6)*V(i6)*sin(delt (i6)-delt(j6)-alf(i6,j6);J(i6,j6)=-Jd1(i6)+V(i6)*Y(i6,j6)*V(j6)*cos(del t(i6)-delt(j6)-alf(i6,j6);N(i6,j6)=-Jd1(i6)-V(i6)*Y(i6,i6)*V(i6)*cos(alf (i6,i6);L(i6,i6)=-Hd1(i6)+V(i6)*

10、Y(i6,i6)*V(i6)*sin(alf (i6,i6);endendend%修正Jacobian矩陣for j9=3for i9=1:nN(i9,j9)=0;L(i9,j9)=0;J(j9,i9)=0;L(j9,i9)=0;endendL(j9,j9)=1;for j9=4for i9=1:nH(i9,j9)=0;N(i9,j9)=0;J(i9,j9)=0;L(i9,j9)=0;H(j9,i9)=0;N(j9,i9)=0;J(j9,i9)=0;L(j9,i9)=0;endendH(j9,j9)=1;L(j9,j9)=1;%Jaco=H N;J L;%Jaco=zeros(2*n,2*n);for i11=1:nfor j11=1:nJaco(2*i11-1,2*j11-1)=H(i11,j11);Jaco(2*i11-1,2*j11)=N(i11,j11);Jaco(2*i11,2*j11-1)=J(i11,j11);Jaco(2*i11,2*