潮流样本.m

潮流样本。能运用于电力系统潮流计算

%============================潮流计算程序(P-Q分解法)及特征值计算程序==============================
%==注意事项== 发电机端带负荷  平衡节点的电压取值为1.00(第130行)
%已考虑到电容器,电抗器等元件,其数值写在支路数据文件的第五列,但其值要化为阻抗值
%--------------------------------------------------------------------------------
%支路原始数据需要：起始节点，终止节点，支路电阻，支路电抗，充电电容(相连的发电机的变比)
%如果相连的支路上有发电机，则在发电机的非标准变比侧的节点号前加上负号，以有所区别；
%第五列就改为发电机的变比。充电电容是起始节点的对地电容（当然也是终止节点的）。
%--------------------------------------------------------------------------------
%负荷数据文件包括:负荷所在节点号，负荷有功，负荷无功
%--------------------------------------------------------------------------------
%发电机数据文件包括：发电机所在节点号，发电机的P，发电机的Q（V）
%对于平衡节点，其P，Q均取为0.其节点的编号可不是最大,
%无穷大母线作为平衡节点时其数据应在发电机数据文件dfdj.mat的第一位
%PV节点的节点号要取为负，Q列取电压幅值



%程序功能：求解多机电力系统特征根，左右特征向量，参与矩阵，机电模式判别，励磁模式
% 判别，特征根对于阻尼系数、励磁系统增益及励磁系统时间常数的灵敏度。
%使用方法：运行潮流程序lxlpl，利用系统参数和潮流计算结果，在MATLAB编辑环境下，
%形成bzm,lpqv,gpqv,qpm四个数据文件；在MATLAB命令窗口输入发电机模型类型矩阵mt
%后，运行lxlmsda2 即可。
%例如对7机系统可取 mt=[0,1,2,3,4,3,5]，（0-无穷大母线；1--D+Q;2--Q;3--D;4-无DQ;5-D+Q+g），
%如有无穷大母线,其标识0必须是mt的第一个元素; mt元素的顺序必须与发电机参数阵dfcan
%中发电机的顺序一致
%无穷大母线的数据只写在发电机的数据文件gpqv
%此程序可以对单机无穷大系统计算线性化系数和特征根以及进行相关分析。 
%=========================================
%调入数据文件(电科院6机系统)
clear
load zhilu1.mat -ascii;   %支路数据文件                   
load fuhe1.mat -ascii;   %负荷数据文件                   
load fadian1.mat -ascii;  %发电机数据文件
load gpm1.mat -ascii;  %发电机参数文件
%load fhtx.mat -ascii; %负荷静态特性数组文件                   
%-----------------------------------------
nm=max(zhilu1(:,2:3)); nn=max(nm');%网络节点总数
bnn=size(zhilu1);      bn=bnn(1,1);%支路数
lnn=size(fuhe1);      ln=lnn(1,1);%负荷节点数
gnn=size(fadian1);     gn=gnn(1,1);%发电机节点数
%-----------------------------------------
%发电机的类型参数
mt1=[4 4 4]; 
%==================BX型PQ分解法=====================
%形成bp,不考虑充电电容及变压器非标准变比对导纳矩阵的影响
%但要考虑支路电阻的影响
ynb=zeros(nn,nn);
for l=1:bn;
   n=abs(zhilu1(l,1));
   m=abs(zhilu1(l,2));
   if n==m
      y1=0;
   else
    y1=1./(zhilu1(l,3)+j*zhilu1(l,4));
    ynb(n,n)=ynb(n,n)+y1;
    ynb(m,m)=ynb(m,m)+y1;
    ynb(n,m)=-y1;
    ynb(m,n)=-y1;
   end
end,
bp=imag(ynb);
%-----------------------------------------
%形成bpp,不考虑支路电阻的影响,
%但要考虑支路充电电容及变压器的非标准变比的影响
ynp=zeros(nn,nn);
for l=1:bn;
   n=abs(zhilu1(l,1));
   m=abs(zhilu1(l,2));
   if n==m
      ynp(n,n)=ynp(n,n)+1./(j*zhilu1(l,5));
   else
      y1=1./(j*zhilu1(l,4));%不考虑支路电阻
      if (zhilu1(l,1)>0&zhilu1(l,2)>0)       %不是变压器支路,是输电线支路 
         ynp(n,n)=ynp(n,n)+y1+j*zhilu1(l,5);
         ynp(m,m)=ynp(m,m)+y1+j*zhilu1(l,5);
         ynp(n,m)=-y1;
         ynp(m,n)=-y1;
      elseif zhilu1(l,1)<0   %非标准变比侧          
         kt=zhilu1(l,5);
         ynp(n,n)=ynp(n,n)+y1/(kt^2);
         ynp(m,m)=ynp(m,m)+y1;
         ynp(n,m)=-(ynp(n,m)+y1/kt);
         ynp(m,n)=ynp(n,m);
      elseif zhilu1(l,2)<0   %非标准变比侧
         kt=zhilu1(l,5);
         ynp(n,n)=ynp(n,n)+y1;
         ynp(m,m)=ynp(m,m)+y1/(kt^2);
         ynp(n,m)=-(ynp(n,m)+y1/kt);
         ynp(m,n)=ynp(n,m);
      end
   end
end,
bpp=imag(ynp);
%---------------------------------------
%求取导纳矩阵
yy=zeros(nn,nn);
for l=1:bn;
   n=abs(zhilu1(l,1));
   m=abs(zhilu1(l,2));
   if n==m
      yy(n,n)=yy(n,n)+1./(j*zhilu1(l,5));
   else
   y1=1./(zhilu1(l,3)+j*zhilu1(l,4));
   if (zhilu1(l,1)>0&zhilu1(l,2)>0)       %不是变压器支路,是输电线支路 
         yy(n,n)=yy(n,n)+y1+j*zhilu1(l,5);
         yy(m,m)=yy(m,m)+y1+j*zhilu1(l,5);
         yy(n,m)=-y1;
         yy(m,n)=-y1;
   elseif zhilu1(l,1)<0   %非标准变比侧          
         kt=zhilu1(l,5);
         yy(n,n)=yy(n,n)+y1/(kt^2);
         yy(m,m)=yy(m,m)+y1;
         yy(n,m)=-(yy(n,m)+y1/kt);
         yy(m,n)=yy(n,m);
    elseif zhilu1(l,2)<0   %非标准变比侧
         kt=zhilu1(l,5);
         yy(n,n)=yy(n,n)+y1;
         yy(m,m)=yy(m,m)+y1/(kt^2);
         yy(n,m)=-(yy(n,m)+y1/kt);
         yy(m,n)=yy(n,m);
      end
   end
end,
%-----------------------------------
%定义变量名
pv=0;          %PV节点的个数
v=ones(nn,1);  %电压
pg=zeros(nn,1);%发电机有功
qg=zeros(nn,1);%发电机无功
pl=zeros(nn,1);%负荷有功
ql=zeros(nn,1);%负荷无功
o=zeros(nn,1); %电压相角
%----------------------------------
%根据节点类型对bp，bpp进行修正原则是:将bp，bpp通取为n阶
%PO迭代:bp中对平衡节点所在的行和列,除对角线元素取作1外,其他元素均取为0.
%其相应的右端元素opv(i)也取为0(见后面计算);
%QV迭代:bpp中对平衡节点和PV节点所在的行和列,除对角线元素取作1外,其他元素均取为0.
%其相应的右端元素oqv(i)也取为0(见后面计算).
%取PQ节点的有功和无功功率和PV节点的有功,电压
for l=1:gnn;%发电机
   if fadian1(l,2)==0&fadian1(l,3)==0 %平衡节点
      f=fadian1(l,1);
      lxl=f;
      v(f)=1.02;
      bp(f,:)=0;bp(:,f)=0;bp(f,f)=1;
      bpp(f,:)=0;bpp(:,f)=0;bpp(f,f)=1;
   elseif abs(fadian1(l,1))>fadian1(l,1)%PV节点
      f=abs(fadian1(l,1));
      pg(f)=fadian1(l,2);
      v(f)=fadian1(l,3);
      bpp(f,:)=0;bpp(:,f)=0;bpp(f,f)=1;
      pv=pv+1;
      pvj(pv)=f;%记录PV节点的节点号
   else
      f=fadian1(l,1);
      pg(f)=fadian1(l,2);
      qg(f)=fadian1(l,3);
   end
end
for l=1:lnn %负荷
   f=fuhe1(l,1);
   pl(f)=fuhe1(l,2);
   ql(f)=fuhe1(l,3);
end;
%===============================================================================
%功率误差
oo=ones(nn,1);
ov=ones(nn,1);
opv=ones(nn,1);
oqv=ones(nn,1);
kop=0;
%--------------------------------------
while(max(abs(opv))>0.000001|max(abs(oqv))>0.000001)
   %---------PO迭代----------------
   kop=kop+1;
   if kop>100
      break;
   end
   
   if max(abs(opv))>0.000001
     av=0;%平均电压
     for i=1:nn
       su=0;
       for l=1:nn
         if yy(i,l)~=0
           do=o(i)-o(l); 
           su=su+v(l)*(real(yy(i,l))*cos(do)+imag(yy(i,l))*sin(do));
         end
       end
       opv(i)=(pg(i)-pl(i)-v(i)*su)/v(i);%计算常数项
       av=av+v(i);
     end
     opv(lxl)=0;%位于后面计算方便,令其维数为n
     av=(av+1.05)/nn;%平均电压
     oo=inv(bp)*opv/av;%计算误差量oo,
     o=o-oo;%修正相角
  end 
  %----------QV迭代---------------
  if max(abs(oqv))>0.000001
     for i=1:nn
       su=0;
       for l=1:nn
          if yy(i,l)~=0
             do=o(i)-o(l);
            su=su+v(l)*(real(yy(i,l))*sin(do)-imag(yy(i,l))*cos(do));
          end
       end
      oqv(i)=(qg(i)-ql(i)-v(i)*su)/v(i);%计算常数项
    end
    for i=1:pv
     l=pvj(i);
     oqv(l)=0;%对PV节点令其右端项为0
    end
   oqv(lxl)=0;
   ov=inv(bpp)*oqv;
   v=v-ov;%修正幅值
 end 
 %----------下一次迭代-------- 
end
%=================潮流已经收敛======================2002.3.28
%计算平衡节点的有功,无功
for i=1:nn
   if yy(lxl,i)~=0
      do=o(lxl)-o(i);
      pg(lxl)=pg(lxl)+v(lxl)*v(i)*(real(yy(lxl,i))*cos(do)+imag(yy(lxl,i))*sin(do));
      qg(lxl)=qg(lxl)+v(lxl)*v(i)*(real(yy(lxl,i))*sin(do)-imag(yy(lxl,i))*cos(do));
   end
end
%计算PV节点的无功功率
for i=1:pv
   k=pvj(i);
   for l=1:nn
      if yy(k,l)~=0
         do=o(k)-o(l);
         qg(k)=qg(k)+v(k)*v(l)*(real(yy(k,l))*sin(do)-imag(yy(k,l))*cos(do));
      end
   end
end
%================潮流计算完毕====================================
%与特征值计算程序衔接2002.3.28
%支路数据文件的衔接
bzm=zeros(bn,8);
for i=1:bn
   bzm(i,1)=1;bzm(i,2)=zhilu1(i,1);bzm(i,3)=zhilu1(i,2);
   bzm(i,4)=0; bzm(i,5)=zhilu1(i,3);bzm(i,6)=zhilu1(i,4);
   if zhilu1(i,1)<0|zhilu1(i,2)<0
      bzm(i,2)=abs(zhilu1(i,1));bzm(i,3)=abs(zhilu1(i,2));
      bzm(i,8)=zhilu1(i,5);
   else
      bzm(i,7)=zhilu1(i,5);
   end
end
%负荷数据的衔接
lpqv=zeros(ln,10);
for i=1:ln
   lpqv(i,1)=fuhe1(i,1);lpqv(i,2)=1;lpqv(i,5)=fuhe1(i,2);lpqv(i,6)=fuhe1(i,3);
   lpqv(i,7)=v(fuhe1(i,1));lpqv(i,8)=o(fuhe1(i,1))*180/pi;
   lpqv(i,9)=lpqv(i,7)*cos(o(fuhe1(i,1)));lpqv(i,10)=lpqv(i,7)*cos(o(fuhe1(i,1)));
end
%发电机节点数据的衔接
gpqv=zeros(gn,10);
for i=1:gn
   gpqv(i,1)=abs(fadian1(i,1));gpqv(i,2)=1;
   gpqv(i,3)=pg(gpqv(i,1));gpqv(i,4)=qg(gpqv(i,1));
   gpqv(i,7)=v(gpqv(i,1));gpqv(i,8)=o(gpqv(i,1))*180/pi;
   gpqv(i,9)=gpqv(i,7)*cos(o(gpqv(i,1)));gpqv(i,10)=gpqv(i,7)*sin(o(gpqv(i,1)));
   if i==lxl %平衡节点
      gpqv(i,2)=1;
   end
   for l=1:pv
      if pvj(l)==abs(fadian1(i,1));%PV节点
         gpqv(i,2)=-1;
      end
   end
end
%=========================衔接部分完成===========================
%--------------------------------------------------------------------------------------------------
ef=[1,1,1];   %各发电机励磁系统的阶数
%统计发电机的各个参数的个数及各发电机的阶数fact
eqpp=zeros(1,2);
ed=zeros(1,2);%edp
edpp=zeros(1,2);
fact=zeros(1,gn);
for i=1:gn           
   if mt1(i)==1       %D+Q
      eqpp(1,1)=eqpp(1,1)+1;
      edpp(1,1)=edpp(1,1)+1;
      fact(i)=5+ef(i);
   elseif mt1(i)==2   %Q
      edpp(1,1)=edpp(1,1)+1;
      fact(i)=4+ef(i);
   elseif mt1(i)==3   %D
      eqpp(1,1)=eqpp(1,1)+1;
      fact(i)=4+ef(i);
   elseif mt1(i)==4
      fact(i)=3+ef(i);
   elseif mt1(1,1)==5 %D+Q+g
      eqpp(1,1)=eqpp(1,1)+1;
      edpp(1,1)=edpp(1,1)+1;
      ed(1,1)=ed(1,1)+1;
      fact(i)=6+ef(i);
   end
end
%---------------------------------------------------------------------------------------------------
%====================================
%原始网络导纳矩阵yy
%====================================
%恒定阻抗负荷模型下的网络导纳矩阵ynl
yl=zeros(nn);
for l=1:ln,
    n=lpqv(l,1);
    yl(n,n)=(lpqv(l,5)-j*lpqv(l,6))./(lpqv(l,7)^2);
end,yl;
ynl=yy+yl;
%--------------------------------
%消去中间节点后的发电机节点间输电导纳矩阵yt
yff=ynl(1:gn,1:gn);
yfl=ynl(1:gn,gn+1:nn);
ylf=ynl(gn+1:nn,1:gn);
yll=ynl(gn+1:nn,gn+1:nn);
yt1=yff-yfl*inv(yll)*ylf; %发电机节点（含无穷大母线）间的导纳矩阵
%--------------------------------

if mt1(1,1)==0            %有无穷大母线时
   yt=yt1(2:gn,2:gn);     %发电机节点间导纳矩阵
   yt0=yt1(2:gn,1);       %无穷大母线与发电机节点间的互导纳向量
   u0=gpqv(1,7);          %无穷大母线电压幅值
   gpqv(1,:)=[];
   mt=mt1(1,2:gn);