matlab isb,牛顿-拉夫逊法进行潮流计算matlab源程序

%本程序的功能是用牛顿-拉夫逊法进行潮流计算 n=input('请输入节点数:n='); nl=input('请输入支路数:nl='); isb=input('请输入平衡母线节电号:isb='); pr=input('请输入误差精度:pr='); B1=input('请输入由支路参数形成的矩阵:B1=');%变压器侧为1，否则为0 B2=input('请输入各节点参数形成的矩阵:B2='); X=input('请输入由节点号及其对地阻抗形成的矩阵:X='); Y=zeros(n);U=zeros(1,n);cta=zeros(1,n);V=zeros(1,n);O=zeros(1,n);S1=zeros(nl); for i=1:n     if X(i,2)~=0;         p=X(i,1);         Y(p,p)=X(i,2);     end end for i=1:nl     if B1(i,6)==0         p=B1(i,1);q=B1(i,2);     else p=B1(i,2);q=B1(i,1);     end     Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));     Y(q,p)=Y(p,q);     Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;     Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; end %求导纳矩阵 G=real(Y);B=imag(Y); for i=1:n     cta(i)=angle(B2(i,3));     U(i)=abs(B2(i,3));     %V(i)=B2(i,4); end for i=1:n     S(i)=B2(i,1)-B2(i,2);     B(i,i)=B(i,i)+B2(i,5); end P=real(S);Q=imag(S); ICT1=0;IT2=1; while IT2~=0     IT2=0;t1=1;t2=1;     for i=1:n         if i~=isb             C(i)=0;             D(i)=0;             for j1=1:n                 C(i)=C(i)+U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));                 D(i)=D(i)+U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));             end             DP(t1)=P(i)-C(i);             t1=t1+1;             if B2(i,6)==2                 DQ(t2)=Q(i)-D(i);                 t2=t2+1;             end         end     end     t1=t1-1;t2=t2-1;     DPQ=[DP';DQ']; %求DP,DQ     for i=1:t1+t2         if abs(DPQ(i))>pr             IT2=IT2+1;         end     end     H=zeros(t1,t1);N=zeros(t1,t2);K=zeros(t2,t1);L=zeros(t2,t2);     for i=1:t1         for j1=1:t1             if j1~=isb&j1~=i                H(i,j1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));            elseif j1~=isb&j1==i                H(i,j1)=U(i)^2*B(i,j1)+D(i);            end         end     end     for i=1:t1         for j1=1:t2             if j1~=isb&j1~=i                 N(i,j1)=0-U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));             elseif j1~=isb&j1==i                 N(i,j1)=0-U(i)^2*G(i,j1)-C(i);             end         end     end     for  i=1:t2         for j1=1:t1             if j1~=isb&j1~=i                 K(i,j1)= U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));             elseif j1~=isb&j1==i                 K(i,j1)=U(i)^2*G(i,j1)-C(i);             end         end     end     for i=1:t2         for j1=1:t2             if j1~=isb&j1~=i                 L(i,j1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));             elseif j1~=isb&j1==i                 L(i,j1)=U(i)^2*B(i,j1)-D(i);             end         end     end     J=[H,N;K,L];%求雅可比矩阵     modify=-J\DPQ;     Dcta=modify([1:t1],:);     t3=U(:,[1:t2]);     DU=diag(t3,0)*modify([t1+1:t1+t2],:);     t4=1;     for i=1:t1         if B2(i,6)~=1         cta(1,i)=cta(1,i)+Dcta(t4,1);         t4=t4+1;         end     end     t5=1;     for i=1:t2         if B2(i,6)==2         U(1,i)=U(1,i)+DU(t5,1);         t5=t5+1;         end     end     ICT1=ICT1+1; end   %修正原值 for i=1:n     UU(i)=U(i)*cos(cta(i))+1i*U(i)*sin(cta(i)); end for p=1:n     c(p)=0;     for q=1:n     c(p)=c(p)+conj(Y(p,q))*conj(UU(q));     end     s(p)=UU(p)*c(p); end disp('--------------------------------------------------------------------------------'); disp('各节点电压U为(节点从小到大排列):'); disp(UU); disp('--------------------------------------------------------------------------------'); disp('各节点电压相角为(节点从小到大排列):'); disp(180*angle(UU)/pi); disp('--------------------------------------------------------------------------------'); disp('按公式计算全部线路功率，结果如下:'); for i=1:nl     if B1(i,6)==0         p=B1(i,1);q=B1(i,2);     else p=B1(i,2);q=B1(i,1);     end     Si(p,q)=UU(p)*(conj(UU(p))*conj(B1(i,4)./2)+(conj(UU(p)*B1(i,5))-conj(UU(q)))*conj(1./(B1(i,3)*B1(i,5))));%各条支路首端功率Si     f=[p,q,Si(p,q)];     disp(f); end for i=1:nl     if B1(i,6)==0         p=B1(i,1);q=B1(i,2);     else p=B1(i,2);q=B1(i,1);     end     Sj(q,p)=UU(q)*(conj(UU(q))*conj(B1(i,4)./2)+(conj(UU(q)./B1(i,5))-conj(UU(p)))*conj(1./(B1(i,3)*B1(i,5))));%各条支路末端功率Sj     f=[q,p,Sj(q,p)];     disp(f); end disp('--------------------------------------------------------------------------------'); disp('各条支路的功率损耗DS为(顺序同您输入B1时一样):'); for i=1:nl     if B1(i,6)==0         p=B1(i,1);q=B1(i,2);     else p=B1(i,2);q=B1(i,1);         end     DS(i)=Si(p,q)+Sj(q,p);%各条支路功率损耗DS     disp(DS(i)); end Sp=0; for i=1:n     Sp=Sp+UU(isb)*conj(Y(isb,i))*conj(UU(i)); end disp('平衡节点的功率:'); disp(Sp);