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

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

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

=xlsread('pqinput.xls','A2:A2');
=xlsread('pqinput.xls','B2:B2');
e=xlsread('pqinput.xls','B4:B4');
=xlsread('pqinput.xls','D3:H100');
=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));
          elseY(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


   为什么呢？因为收敛性比较差吗？
Roczhang 发表于 2011-7-22 12:35 https://tech.cepsc.com/images/common/back.gif

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

查看完整版本: 潮流计算P-Q分解法MATLAB语言程序（希望对大家有用哦）