《PQ分解法程序.doc》由会员分享，可在线阅读，更多相关《PQ分解法程序.doc（12页珍藏版）》请在金锄头文库上搜索。
1、PQ分解法程序x=xlsread(潮流参数.xls,A2:A2);y=xlsread(潮流参数.xls,B2:B2);e=xlsread(潮流参数.xls,B7:B7);point=xlsread(潮流参数.xls,D3:H6);zhilu=xlsread(潮流参数.xls,J3:R6);TYPE=zeros(x,1);U=zeros(x,1);a=zeros(x,1);P=zeros(x,1);Q=zeros(x,1);I=zeros(y,1);J=zeros(y,1);Rij=zeros(y,1);Xij=zeros(y,1);Zij=Rij+j*Xij; Y=zeros(x); G=ze
2、ros(x); B=zeros(x); B0=zeros(y,1);RT=zeros(y,1);XT=zeros(y,1);ZT=RT+j*XT;KT=zeros(y,1); %-矩阵赋初值：TYPE=point(:,1);U=point(:,2);a=point(:,3);P=point(:,4);Q=point(:,5);I=zhilu(:,1);J=zhilu(:,2);Rij=zhilu(:,3);Xij=zhilu(:,4);Zij=Rij+j*Xij;B0=zhilu(:,5);RT=zhilu(:,6);XT=zhilu(:,7);ZT=RT+j*XT;KT=zhilu(:,8)
3、;W=zhilu(:,9);%-求节点导纳矩阵Yfor m=1:y if KT(m)=0 Y(I(m),J(m)=-1/Zij(m); Y(J(m),I(m)=-1/Zij(m); else Y(I(m),J(m)=-1/(KT(m)*ZT(m); Y(J(m),I(m)=-1/(KT(m)*ZT(m); endendfor m=1:x for n=1:y if KT(n)=0 if(I(n)=m|J(n)=m) Y(m,m)=Y(m,m)-Y(I(n),J(n)+j*B0(n)/2; end else if W=0 if I(n)=m Y(m,m)=Y(m,m)-Y(I(n),J(n)+(K
4、T(n)-1)/KT(n)*(1/ZT(n); elseif J(n)=m Y(m,m)=Y(m,m)-Y(I(n),J(n)+(1-KT(n)/(KT(n)2)*(1/ZT(n); end else if I(n)=m Y(m,m)=Y(m,m)-Y(I(n),J(n)+(1-KT(n)/(KT(n)2)*(1/ZT(n); elseif J(n)=m Y(m,m)=Y(m,m)-Y(I(n),J(n)+(KT(n)-1)/KT(n)*(1/ZT(n); else Y(m,m)=Y(m,m); end end endendendG=real(Y);%-求B矩阵及其逆矩阵B1B=imag(Y)
5、; ph=find(TYPE(:,1)=3); BB=B; BB(:,ph)=; BB(ph,:)=; B1=BB; B1=inv(B1);disp(B矩阵)disp(BB)disp(B逆矩阵)disp(B1) %-%求B及其逆矩阵B2phpv=find(TYPE(:,1)1); BB=B; BB(:,phpv)=; BB(phpv,:)=; B2=BB; B2=inv(B2);disp(BB)disp(B2) %-计算各节点有功功率不平衡量deltaPi k=0; kp=1; kq=1;while(kp=0)|(kq=0)&k=20) kp=1; kq=1; notph=find(TYPE(
6、:,1)3); deltaPi=zeros(x-1,1); pq=find(TYPE(:,1)=1); pqnum=size(B2); pqnum=pqnum(1); deltaQi=zeros(pqnum,1); for m=1:(x-1) sum1=0; for n=1:x sum1=sum1+U(notph(m)*U(n)*(G(notph(m),n)*cos(a(notph(m)-a(n)+B(notph(m),n)*sin(a(notph(m)-a(n); end deltaPi(m)=P(notph(m)-sum1; enddisp(有功功率不平衡量)disp(k)disp(del
7、taPi) max1=max(abs(deltaPi); if max1=e kp=0; if kq=0 break else for m=1:pqnum sum2=0; for n=1:xsum2=sum2+U(pq(m)*U(n)*(G(pq(m),n)*sin(a(pq(m)-a(n)-B(pq(m),n)*cos(a(pq(m)-a(n); end deltaQi(m)=Q(pq(m)-sum2; end max2=max(abs(deltaQi); if max2=e kq=0; if kp=0 break else k=k+1; end else Uq=U; Uq(phpv)=;
8、Upq=Uq; deltaU=-B2*(deltaQi./Upq); for m=1:pqnum U(pq(m)=U(pq(m)+deltaU(m); end kp=1; k=k+1; end end else Up=U; Up(ph)=; Unotph=Up; deltaa=(-B1*(deltaPi./Unotph)./Unotph); for m=1:(x-1) a(notph(m)=a(notph(m)+deltaa(m); enddisp(电压相角修正)disp(k)disp(deltaa)disp(a) kq=1; for m=1:pqnum sum2=0; for n=1:xsu
9、m2=sum2+U(pq(m)*U(n)*(G(pq(m),n)*sin(a(pq(m)-a(n)-B(pq(m),n)*cos(a(pq(m)-a(n); end deltaQi(m)=Q(pq(m)-sum2; end disp(wu功功率不平衡量)disp(k)disp(deltaQi) Uq=U; Uq(phpv)=; Upq=Uq; deltaU=-B2*(deltaQi./Upq); for m=1:pqnum U(pq(m)=U(pq(m)+deltaU(m); enddisp(电压修正)disp(k)disp(deltaU)disp(U) kp=1; k=k+1; endend
10、 sum3=0+j*0; for m=1:x sum3=sum3+conj(Y(ph,m)*(U(m)*cos(a(m)-i*U(m)*sin(a(m) end Sph=(U(ph)*cos(a(ph)+j*U(ph)*sin(a(ph)*sum3; %-求线路功率Sij和Sji Sij=zeros(y,1); Sji=zeros(y,1); for m=1:y if KT(m)=0 Sij(m)=(U(I(m)*cos(a(I(m)+i*U(I(m)*sin(a(I(m)*(U(I(m)*cos(a(I(m)-i*U(I(m)*sin(a(I(m)*(-i*B0(m)/2)+(U(I(m)*
11、cos(a(I(m)-i*U(I(m)*sin(a(I(m)-U(J(m)*cos(a(J(m)+i*U(J(m)*sin(a(J(m)*conj(-Y(I(m),J(m); Sji(m)=(U(J(m)*cos(a(J(m)+i*U(J(m)*sin(a(J(m)*(U(J(m)*cos(a(J(m)-i*U(J(m)*sin(a(J(m)*(-i*B0(m)/2)+(U(J(m)*cos(a(J(m)-i*U(J(m)*sin(a(J(m)-U(I(m)*cos(a(I(m)+i*U(I(m)*sin(a(I(m)*conj(-Y(I(m),J(m); else if W=1 Sij(m)=(U(I(m)*cos(a(I(m)+i*U(I(m)*sin(a(I(m)*(U(I(m)*cos(a(I(m)-i*U(I(m)*sin(a(I(m)*(1-KT(m)/(KT(m)2*conj(1/ZT(m)+(U(I(m)*cos(a(I(m)-i*U(I(m)*sin(a(I(m)-U(J(m)*cos(a(J(m)+i*U(J(m)*sin(a(J(m)*conj(-Y(I(m),J(m); Sji(m)=(U(J(m)*cos(a(J(m)+i*U(J(m)*sin(a(J(m)*(U(J(m)*c
