MATLAB.rar_潮流计算

程序代码如下: 111111. %读入数据 clc clear filename='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:m1 if a(i1,1)==0 j1=j1+1; b(j1)=i1; end end b; %矩阵分块 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:b(5)-1,1:n1); a5=a(b(5)+1:b(6)-1,1:n1); %设置初值 vcz=1; dcz=0; kmax=20; k1=0; %求节点导纳矩阵 a11=zeros(4,6); for i0=1:3 for j0=1:6 a11(i0,j0)=a1(i0,j0); a11(4,j0)=a2(1,j0); end end a11; 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); for k=1:4 i2=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)~=0 k0=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; else G(j2,i2)=-GIJ; G(i2,j2)=G(j2,i2); B(j2,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; end G(j2,j2)=G(j2,j2)+GIJ; B(j2,j2)=B(j2,j2)+BIJ+b; end G; B; B=B.*i; Yf=G+B Y=abs(Yf); alf=angle(Yf); %赋Jacobian矩阵参数 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=a4(:,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)); end delt=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:kmax disp(['第' int2str(j5) '次计算结果']) if k>=kmax break end for i10=1:4 Pd(i10)=0; Qd(i10)=0; for j10=1:n Pd(i10)=Pd(i10)+V(i10)*Y(i10,j10)*V(j10)*cos(delt(i10)-delt(j10)-alf(i10,j10)); Qd(i10)=Qd(i10)+V(i10)*Y(i10,j10)*V(j10)*sin(delt(i10)-delt(j10)-alf(i10,j10)); end end for i4=1:3 dP(i4)=P(i4)-Pd(i4); end for j4=1:2 dQ(j4)=Q(j4)-Qd(j4); end A=Jaco(V,delt,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=A\B; for i16=1:n ddelt(i16)=X(2*i16-1); dV(i16)=X(2*i16)*V(i16); end V=V+dV delt=delt+ddelt else break end disp('----------------') end %流氓算法 % for ii=1:2 % V(ii)=V(ii)+dV(ii); % end % V 222222. function A=Jaco(V,delt,n,Y,alf) %计算Jacobian矩阵 for i7=1:n Hd1(i7)=0; Jd1(i7)=0; for j7=1:n Hd1(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 end for i6=1:n for j6=1:n if i6~=j6 H(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*sin(delt(i6)-delt(j6)-alf(i6,j6)); N(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6)); J(i6,j6)=-N(i6,j6); L(i6,j6)=H(i6,j6); else H(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(delt(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)*Y(i6,i6)*V(i6)*sin(alf(i6,i6)); end end end %修正Jacobian矩阵 for j9=3 for i9=1:n N(i9,j9)=0; L(i9,j9)=0; J(j9,i9)=0; L(j9,i9)=0; end end L(j9,j9)=1; for j9=4 for i9=1:n H(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; end end H(j9,j9)=1; L(j9,j9)=1; %Jaco=[H N;J L]; %Jaco=zeros(2*n,2*n); for i11=1:n for j11=1:n Jaco(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*j11)=L(i11,j11); end end A=Jaco; 33333. 数据: 4 4 1.05 0.00001 0 1 1 2 0.1 0.40 0.01528 2 1 4 0.12 0.50 0.01920 3 2 4 0.08 0.40 0.01413 0 1 1 3 0 0.3 0.90909091 0 0 1 1 0 0 0.30 0.18 2 2 0 0 0.55 0.13 3 3 0.5 0 0 0 0 1 3 1.10 0 0 0