mypq30(未用无功越限函数)(57可用).m

电力系统潮流计算

clc;
clear;
format long;
n=118;%系统节点数
nbranch=0;%支路数据的行数
%****输入网络的数据***********************************正确改路径即可
%+++++++++++++输入 a5****************************
load a118.dat;
data0=a118;
%***************输入b5************************
load b118.dat;
data1=b118;
%***********************************************
fr_ij=data0(:,2);
to_ij=data0(:,3);
r_ij=data0(:,4);
x_ij=data0(:,5);
b_ij=data0(:,6)/2;
k_ij=data0(:,7);
[nbranch,wuyong]=size(fr_ij);%支路数据维数
pg_i=data1(:,2);%节点注入有功功率
pl_i=data1(:,3);%节点负荷消耗有功功率
qg_i=data1(:,4);%节点注入无功功率
ql_i=data1(:,5);%节点负荷消耗无功功率
UAngle=data1(:,6);%节点电压相角
U=data1(:,7);%电压幅值
type_i=data1(:,8);%节点类型
Capac=data1(:,9);%并联电容
Qmax=data1(:,10)%发电机无功上限,换成标幺值
Qmin=data1(:,11)%发电机无功下限,换成标幺值
%%%%%%%%%%%%%%%%%%%输入完成%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=0;%用于记PQ数
r=0;%用于记PV数
NoSlack=1;
NoPQ=zeros(n,1);
NoPV=zeros(n,1);
for a=1:n  %计算m  PV节点个数为(n-m)=r
    if type_i(a)==2%PQ节点typei==2
        m=m+1;
        NoPQ(m)=a;%将PQ节点的编号存到数组NoPQ(m)中
    end
    if type_i(a)==3
        NoSlack=a;%将平衡节点编号存入NoSlack中
    end
    if type_i(a)==1;%PV节点typei==1
        r=r+1;
        NoPV(r)=a;%将PV节点的编号存到数组NoPV(r)中
    end
end
m=m+1;
FlagPV=zeros(r,1);%记上次无功越限的节点值为1
FlagPVc=zeros(r,1);%处理无功越限时用
%ConQ;%PV转PQ时用来比较
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%(1)***************形成节点导纳阵********************
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%n=size(type_i)-1;
P=zeros(n,1);
Q=zeros(n,1);
y=zeros(n,n);
B1=zeros(n,n);%B'
B2=10^8*eye(n,n);%B''
for a=1:nbranch
       if k_ij(a)==1           %非变压器支路
         y(fr_ij(a),to_ij(a))=y(fr_ij(a),to_ij(a))-1/(r_ij(a)+x_ij(a)*i);%非对角元
         y(to_ij(a),fr_ij(a))=y(fr_ij(a),to_ij(a));
         %B1(fr_ij(a),to_ij(a))=-imag(y(fr_ij(a),to_ij(a)));
         B1(fr_ij(a),to_ij(a))=B1(fr_ij(a),to_ij(a))-x_ij(a)/(r_ij(a)^2+x_ij(a)^2);
         B1(to_ij(a),fr_ij(a))=B1(fr_ij(a),to_ij(a));
         B2(fr_ij(a),to_ij(a))=B2(fr_ij(a),to_ij(a))-1/x_ij(a);
         B2(to_ij(a),fr_ij(a))=B2(fr_ij(a),to_ij(a));
         
         y(fr_ij(a),fr_ij(a))=y(fr_ij(a),fr_ij(a))+1/(r_ij(a)+x_ij(a)*i)+b_ij(a)*i;%对角元
         y(to_ij(a),to_ij(a))=y(to_ij(a),to_ij(a))+1/(r_ij(a)+x_ij(a)*i)+b_ij(a)*i;
         B1(fr_ij(a),fr_ij(a))=B1(fr_ij(a),fr_ij(a))-B1(fr_ij(a),to_ij(a));
         B1(to_ij(a),to_ij(a))=B1(to_ij(a),to_ij(a))-B1(fr_ij(a),to_ij(a));
         B2(fr_ij(a),fr_ij(a))=B2(fr_ij(a),fr_ij(a))-B2(fr_ij(a),to_ij(a))-b_ij(a);%此处“-”
         B2(to_ij(a),to_ij(a))=B2(to_ij(a),to_ij(a))-B2(fr_ij(a),to_ij(a))-b_ij(a);
       end;
       if k_ij(a)~=1                         %变压器支路
         y(fr_ij(a),to_ij(a))=y(fr_ij(a),to_ij(a))-1/((r_ij(a)+x_ij(a)*i)*k_ij(a));
         y(to_ij(a),fr_ij(a))=y(fr_ij(a),to_ij(a));%非对角元
         B1(fr_ij(a),to_ij(a))=B1(fr_ij(a),to_ij(a))-x_ij(a)/((r_ij(a)^2+x_ij(a)^2)*k_ij(a));
         B1(to_ij(a),fr_ij(a))=B1(fr_ij(a),to_ij(a));
         B2(fr_ij(a),to_ij(a))=B2(fr_ij(a),to_ij(a))-1/(x_ij(a)*k_ij(a));
         B2(to_ij(a),fr_ij(a))=B2(fr_ij(a),to_ij(a));
         
         y(fr_ij(a),fr_ij(a))=y(fr_ij(a),fr_ij(a))+1/((r_ij(a)+x_ij(a)*i)*k_ij(a)*k_ij(a))+b_ij(a)*i;
         y(to_ij(a),to_ij(a))=y(to_ij(a),to_ij(a))+1/(r_ij(a)+x_ij(a)*i)+b_ij(a)*i;  %对角元
         B1(fr_ij(a),fr_ij(a))=B1(fr_ij(a),fr_ij(a))-B1(fr_ij(a),to_ij(a))/k_ij(a);
         B1(to_ij(a),to_ij(a))=B1(to_ij(a),to_ij(a))-B1(fr_ij(a),to_ij(a))*k_ij(a);
         B2(fr_ij(a),fr_ij(a))=B2(fr_ij(a),fr_ij(a))+1/(x_ij(a)*k_ij(a)*k_ij(a))-b_ij(a);%此处“+”
         B2(to_ij(a),to_ij(a))=B2(to_ij(a),to_ij(a))+1/x_ij(a)-b_ij(a);
       end; 
  end;
  for a=1:n%加上并联电容器
      if Capac(a)~=0
        y(a,a)=y(a,a)+Capac(a)*i;
        B1(a,a)=B1(a,a)-Capac(a);
        B2(a,a)=B2(a,a)-Capac(a);
      end
  end
y;%节点导纳矩阵

%%%%%%%%%%%%%%%形成完毕%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% 形成B'%%%%%%%%%%%%%%%%%%
    B1(NoSlack,NoSlack)=B1(NoSlack,NoSlack)+10^8;%划去B'平衡节点
for a=1:(m-1)
    B2(NoPQ(a),NoPQ(a))=B2(NoPQ(a),NoPQ(a))-10^8;%划去B''PV节点
end 

P=pg_i-pl_i;
Q=qg_i-ql_i;

%(2)%%%%%%%%%%%%平启动%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%U=ones(n,1) a5.txt中已经给出了
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%迭代循环部分%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DetP=zeros(n,1);%用于P不平衡量
DetQ=zeros(n,1);%用于Q不平衡量
k=0;
Uset=data1(:,7);%可以挑出电压设定值
Alpha=0.5;%减小数值振荡
while k<30
  %%%%%%%%解P方程%%%%%%%%%%%%%%%%%% 
  DetP=dP(y,U,UAngle,NoPQ,NoPV,m,P);
    for a=1:n
     DetP(a)=DetP(a)/U(a);
    end
  DetUAngle=B1\DetP;
  for a=1:(n-1)
     UAngle(a)=UAngle(a)+DetUAngle(a)/U(a);
  end   
 %%%%%%%%%解Q方程%%%%%%%%%%%%%%%%%%  
  DetQ=dQ(y,U,UAngle,NoPQ,m,Q);
    for a=1:n
      DetQ(a)=DetQ(a)/U(a);
    end
  DetU=B2\DetQ;
  for a=1:(n-1)    
    U(a)=U(a)+DetU(a);
  end
  % for a=1:(m-1)
   %  U(NoPQ(a))=U(NoPQ(a))+DetU(NoPQ(a)) ;   
   %end
       
  
   if (max(abs(DetP))<0.00001)&(max(abs(DetQ))<0.00001)%满足精度停止迭代
     break
   end
   k=k+1%计迭代次数
  DetP;
  DetQ;   
  U;
  UAngle;
  %%%%%%%%%%%%无功越限处理%%%%%%%%%%%%%%%%
  QPV=inQPV(y,U,UAngle,NoPV,m,Q);%发电机的注入功率
  ConQ=QPV-ql_i;
  %[B2,FlagPV,U,Q]=ExchPVPQ(ConQ,Q,data1,B2,NoPV,FlagPV,U);
 for a=1:r
   if FlagPV(a)==0
            if ConQ(NoPV(a))<Qmin(NoPV(a))%越下限
                Q(NoPV(a))=Qmin(NoPV(a)); 
                B2(NoPV(a),NoPV(a))=B2(NoPV(a),NoPV(a))-10^8;
                FlagPV(a)=2;
            end
             if ConQ(NoPV(a))>Qmax(NoPV(a))%越下限
                Q(NoPV(a))=Qmax(NoPV(a));
                B2(NoPV(a),NoPV(a))=B2(NoPV(a),NoPV(a))-10^8;
                FlagPV(a)=1;
            end
            break
      end      
     if FlagPV(a)==1
             if U(NoPV(a))>Uset(NoPV(a))
                 if ConQ(NoPV(a))<Qmin(NoPV(a))
                    Q(NoPV(a))=Qmin(NoPV(a));
                    FlagPV(a)=2;
                end
                 if (ConQ(NoPV(a))>=Qmin(NoPV(a)))&(ConQ(NoPV(a))<=Qmax(NoPV(a)))
                    U(NoPV(a))=U(NoPV(a))+Alpha*(Uset(NoPV(a))-U(NoPV(a)));
                    FlagPV(a)=0;
                    B2(NoPV(a),NoPV(a))=B2(NoPV(a),NoPV(a))+10^8;
                end
             end
             break
       end         
      if FlagPV(a)==2
             if U(NoPV(a))<Uset(NoPV(a))
                 if ConQ(NoPV(a))>Qmax(NoPV(a))
                    Q(NoPV(a))=Qmax(NoPV(a));
                    FlagPV(a)=1;
                end
                 if (ConQ(NoPV(a))>=Qmin(NoPV(a)))&(ConQ(NoPV(a))<=Qmax(NoPV(a)))
                    U(NoPV(a))=U(NoPV(a))+Alpha*(Uset(NoPV(a))-U(NoPV(a)));
                    FlagPV(a)=0;
                    B2(NoPV(a),NoPV(a))=B2(NoPV(a),NoPV(a))+10^8;
                end
             end  
             break
       end     
           
end  
  FlagPV
end