潮流计算P-Q分解法MATLAB语言程序（希望对大家有用哦）

Roczhang · 发表于 2011-7-15 17:17:17

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

您需要登录才可以下载或查看，没有账号？立即加入

×

[x]=xlsread('pqinput.xls','A2:A2');
[y]=xlsread('pqinput.xls','B2:B2');
e=xlsread('pqinput.xls','B4:B4');
[point]=xlsread('pqinput.xls','D3:H100');
[zhilu]=xlsread('pqinput.xls','J3:R100');
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=zeros(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);
W=zhilu(:,9);
%------------------------------求节点导纳矩阵Y
for 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));
end
end
for 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))+(KT(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
end
end
end
G=real(Y);
%-----------------------求B'矩阵及其逆矩阵B1
B=imag(Y);
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''及其逆矩阵B2
phpv=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(:,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;
end
disp('有功功率不平衡量')
disp(k)
disp(deltaPi)
max1=max(abs(deltaPi));
if max1<=e
   kp=0;
   if kq==0
   break
else
   for m=1:pqnum
      sum2=0;
      for n=1:x
sum2=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)=[];
   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);
   end
disp('电压相角修正')
disp(k)
disp(deltaa)
disp(a)
   kq=1;
   for m=1:pqnum
      sum2=0;
      for n=1:x
sum2=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);
end
disp('电压修正')
disp(k)
disp(deltaU)
disp(U)
kp=1;
k=k+1;
end
end
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))*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))*cos(a(J(m)))-i*U(J(m))*sin(a(J(m))))*((KT(m)-1)/KT(m)*conj(1/ZT(m)))+(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    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))))*((KT(m)-1)/KT(m)*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))*cos(a(J(m)))-i*U(J(m))*sin(a(J(m))))*((1-KT(m))/(KT(m))^2*conj(1/ZT(m)))+(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))));
      end
   end
end
deltaSij=Sij+Sji;
S=zeros(x,1)+i*zeros(x,1);
for b=1:x
for m=1:y
      if I(m)==b
         S(b)=S(b)+Sij(m);
      else
         if J(m)==b
         S(b)=S(b)+Sji(m);
      else
         S(b)=S(b);
      end
  end
end
end
P=real(S);
Q=imag(S);
sumdeltaS=sum(S);
fid=fopen('shiyan
fprintf(fid,'          ******************************电力系统06 白云霄 0603010213********************************\n');jieguo.txt','wt');
fprintf(fid,'                                     *****************潮流上机输出结果*************\n');
fprintf(fid,'                                     迭代次数k为: %d \n',k);
fprintf(fid,'                                     ================================================\n');
fprintf(fid,'                                     平衡节点%d的复功率Sph为: %f+j*(%f) \n',ph,real(Sph),imag(Sph));
fprintf(fid,'                                     ================================================\n');
fprintf(fid,'                                     节点电压U为:  \n');
for m=1:x
fprintf(fid,'                                  第%d个节点电压: %f\n',m,U(m));
end
fprintf(fid,'                                     ================================================\n');
fprintf(fid,'                                     节点相角a为:  \n');
for m=1:x
fprintf(fid,'                                  第%d个节点相角: %f\n',m,a(m));
end
fprintf(fid,'                                     ================================================\n');
fprintf(fid,'                                     节点复功率S为:  \n');
for m=1:x
fprintf(fid,'                                  第%d个节点复功率: %f+i*(%f)\n',m,real(S(m)),imag(S(m)));
end
fprintf(fid,'                                     ================================================\n');
fprintf(fid,'                                     节点有功功率P为:  \n');
for m=1:x
fprintf(fid,'                                  第%d个节点有功功率: %f\n',m,P(m));
end
fprintf(fid,'                                     ================================================\n');
fprintf(fid,'                                     节点无功功率Q为:  \n');
for m=1:x
fprintf(fid,'                               第%d个节点无功功率Q: %f\n',m,Q(m));
end
fprintf(fid,'                                     ================================================\n');
fprintf(fid,'                                     线路功率Sij和Sji为:  \n');
for m=1:y
fprintf(fid,'                               节点%d到节点%d的功率为:  %f+i*(%f) \n',I(m),J(m),real(Sij(m)),imag(Sij(m)));
fprintf(fid,'                                  节点%d到节点%d的功率为:  %f+i*(%f) \n',J(m),I(m),real(Sji(m)),imag(Sji(m)));
end
fprintf(fid,'                                     ================================================\n');
fprintf(fid,'                                     网络总损耗sumdeltaS为: %f+i*(%f) \n',real(sumdeltaS),imag(sumdeltaS));
fprintf(fid,'                                     ================================================\n');
fprintf(fid,'                                     线路功率损耗deltaSij为:  \n');
for m=1:y
fprintf(fid,'                               %d--%d线路的功率损耗为:  %f+i*(%f) \n',I(m),J(m),real(deltaSij(m)),imag(deltaSij(m)));
end

davidmacintosh · 发表于 2011-7-15 17:59:27

这是楼主自己写的？？华硕PQ法至今还是搞不太懂啊……

Roczhang · 发表于 2011-7-17 11:18:47

回复 2# davidmacintosh
其实是前一段时间做课程设计时候在网上搜集的，觉得挺好用的就分享了，我们这门课下学期才开。华电研究生？大有前途啊。

davidmacintosh · 发表于 2011-7-17 19:47:19

回复 3# Roczhang

偶不是华电的哦，上电本科的。前途渺茫，明年毕业了。

facat · 发表于 2011-7-22 10:26:39

PQ法现在已经很少用了

Roczhang · 发表于 2011-7-22 12:35:03

回复 5# facat

为什么呢？因为收敛性比较差吗？

facat · 发表于 2011-7-22 13:17:24

回复  facat
7 I* M1 d( ]! u, Y+ X# z
- p* X' b1 i/ ~; Z. q) O9 [- r& t: o# c  ^5 Z- ?7 h7 @9 M) ~
为什么呢？因为收敛性比较差吗？
* |% A  x1 j9 J1 x. SRoczhang 发表于 2011-7-22 12:35
1 L9 g/ O' r5 q

登录/注册后可看大图

对
当时研究出PQ法是因为计算机的性能无法满足大规模电网计算的要求，现在计算机的性能已经非常富余，没必要用PQ了。

账号		自动登录	找回密码
密码			立即加入

[讨论] 潮流计算P-Q分解法MATLAB语言程序（希望对大家有用哦）

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

相关帖子

浏览过的版块