电力系统短路电流计算_matlab程序

电力系统短路电流计算_matlab程序 %电力系统极坐标下的牛顿-拉夫逊法潮流计算 disp('电力系统极坐标下的牛顿-拉夫逊法潮流计算:'); clear n=input('请输入结点数:n='); n1=input('请输入PV结点数:n1='); n2=input('请输入PQ结点数:n2='); isb=input('请输入平衡结点:isb='); pr=input('请输入精确度:pr='); K=input('请输入变比矩阵看:K='); C=input('请输入支路阻抗矩阵:C='); y=input('请输入支路导纳矩阵:y='); U=input('请输入结点电压矩阵:U='); S=input('请输入各结点的功率:S='); Z=zeros(1,n);N=zeros(n1+n2,n2);L=zeros(n2,n2);QT1=zeros(1,n1+n2); for m=1:n for R=1:n C(m,m)=C(m,m)+y(m,R); if K(m,R)~=0 C(m,m)=C(m,m)+1/((K(m,R)*C(m,R))/(K(m,R)-1)); C(R,R)=C(R,R)+1/((K(m,R)^2*C(m,R))/(1-K(m,R))); C(m,R)=C(m,R)*K(m,R); C(R,m)=C(m,R); end end end for m=1:n for R=1:n if m~=R Z(m)=Z(m)+1/C(m,R); end end end for m=1:n for R=1:n if m==R Y(m,m)=C(m,m)+Z(m); else Y(m,R)=-1/C(m,R); end end end disp('结点导纳矩阵:'); disp(Y); disp('迭代中的雅克比矩阵:'); G=real(Y); B=imag(Y); O=angle(U); U1=abs(U); k=0; PR=1; P=real(S); Q=imag(S); while PR>pr for m=1:n2 UD(m)=U1(m); end for m=1:n1+n2 for R=1:n PT(R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R))); end PT1(m)=sum(PT); PP(m)=P(m)-PT1(m); PP1(k+1,m)=PP(m); end for m=1:n2 for R=1:n QT(R)=U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R))); end QT1(m)=sum(QT); QQ(m)=Q(m)-QT1(m); QQ1(k+1,m)=QQ(m); end PR1=max(abs(PP)); PR2=max(abs(QQ)); PR=max(PR1,PR2); for m=1:n1+n2 for R=1:n1+n2 if m==R H(m,m)=U1(m)^2*B(m,m)+QT1(m); else H(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R))); end end end for m=1:n1+n2 for R=1:n2 if m==R N(m,m)=-U1(m)^2*G(m,m)-PT1(m); else N(m,R)=-U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R))); end end end for m=1:n2 for R=1:n1+n2 if m==R J(m,m)=U1(m)^2*G(m,m)-PT1(m); else J(m,R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R))); end end end for m=1:n2 for R=1:n2 if m==R L(m,m)=U1(m)^2*B(m,m)-QT1(m); else L(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R))); end end end JJ=[H N;J L]; disp(JJ); PQ=[PP';QQ']; DA=-inv(JJ)*PQ; DA1=DA'; for m=1:n1+n2 OO(m)=DA1(m); end for m=n:n1+n2+n2 UU1(m-n1-n2)=DA1(m); end UD2=diag(UD); UU=UU1*UD2; for m=1:n1+n2 O(m)=O(m)+OO(m); end for m=1:n2 U1(m)=U1(m)+UU(m); end for m=1:n1+n2 o(k+1,m)=180/pi*O(m); end for m=1:n2 u(k+1,m)=U1(m); end k=k+1; end for m=1:n b(m)=U1(m)*cos(O(m)); c(m)=U1(m)*sin(O(m)); end U=b+i*c; for R=1:n PH1(R)=U(isb)*conj(Y(isb,R))*conj(U(R)); end PH=sum(PH1); for m=1:n for R=1:n if m~=R C1(m,R)=1/C(m,R); else C1(m,m)=C(m,m); end end end for m=1:n for R=1:n if (C(m,R)~=inf)&(m~=R) SS(m,R)=U1(m)^2*conj(C1(m,m))+U(m)*(conj(U(m))-conj(U(R)))*conj(C1(m,R)); end end end disp('迭代中的?P:');disp(PP1); disp('迭代中的?Q:');disp(QQ1); disp('迭代中相角:');disp(o); disp('迭代中电压的模:');disp(u); disp('平衡结点的功率:');disp(PH); disp('全部线路功率分布:');disp(SS);