%============================潮流計(jì)算程序(P-Q分解法)及特征值計(jì)算程序==============================
%==注意事項(xiàng)== 發(fā)電機(jī)端帶負(fù)荷  平衡節(jié)點(diǎn)的電壓取值為1.00(第130行)
%已考慮到電容器,電抗器等元件,其數(shù)值寫在支路數(shù)據(jù)文件的第五列,但其值要化為阻抗值
%--------------------------------------------------------------------------------
%支路原始數(shù)據(jù)需要：起始節(jié)點(diǎn)，終止節(jié)點(diǎn)，支路電阻，支路電抗，充電電容(相連的發(fā)電機(jī)的變比)
%如果相連的支路上有發(fā)電機(jī)，則在發(fā)電機(jī)的非標(biāo)準(zhǔn)變比側(cè)的節(jié)點(diǎn)號(hào)前加上負(fù)號(hào)，以有所區(qū)別；
%第五列就改為發(fā)電機(jī)的變比。充電電容是起始節(jié)點(diǎn)的對(duì)地電容（當(dāng)然也是終止節(jié)點(diǎn)的）。
%--------------------------------------------------------------------------------
%負(fù)荷數(shù)據(jù)文件包括:負(fù)荷所在節(jié)點(diǎn)號(hào)，負(fù)荷有功，負(fù)荷無功
%--------------------------------------------------------------------------------
%發(fā)電機(jī)數(shù)據(jù)文件包括：發(fā)電機(jī)所在節(jié)點(diǎn)號(hào)，發(fā)電機(jī)的P，發(fā)電機(jī)的Q（V）
%對(duì)于平衡節(jié)點(diǎn)，其P，Q均取為0.其節(jié)點(diǎn)的編號(hào)可不是最大,
%無窮大母線作為平衡節(jié)點(diǎn)時(shí)其數(shù)據(jù)應(yīng)在發(fā)電機(jī)數(shù)據(jù)文件dfdj.mat的第一位
%PV節(jié)點(diǎn)的節(jié)點(diǎn)號(hào)要取為負(fù)，Q列取電壓幅值



%程序功能：求解多機(jī)電力系統(tǒng)特征根，左右特征向量，參與矩陣，機(jī)電模式判別，勵(lì)磁模式
% 判別，特征根對(duì)于阻尼系數(shù)、勵(lì)磁系統(tǒng)增益及勵(lì)磁系統(tǒng)時(shí)間常數(shù)的靈敏度。
%使用方法：運(yùn)行潮流程序lxlpl，利用系統(tǒng)參數(shù)和潮流計(jì)算結(jié)果，在MATLAB編輯環(huán)境下，
%形成bzm,lpqv,gpqv,qpm四個(gè)數(shù)據(jù)文件；在MATLAB命令窗口輸入發(fā)電機(jī)模型類型矩陣mt
%后，運(yùn)行l(wèi)xlmsda2 即可。
%例如對(duì)7機(jī)系統(tǒng)可取 mt=[0,1,2,3,4,3,5]，（0-無窮大母線；1--D+Q;2--Q;3--D;4-無DQ;5-D+Q+g），
%如有無窮大母線,其標(biāo)識(shí)0必須是mt的第一個(gè)元素; mt元素的順序必須與發(fā)電機(jī)參數(shù)陣dfcan
%中發(fā)電機(jī)的順序一致
%無窮大母線的數(shù)據(jù)只寫在發(fā)電機(jī)的數(shù)據(jù)文件gpqv
%此程序可以對(duì)單機(jī)無窮大系統(tǒng)計(jì)算線性化系數(shù)和特征根以及進(jìn)行相關(guān)分析。 
%=========================================
%調(diào)入數(shù)據(jù)文件(電科院6機(jī)系統(tǒng))
clear
load zhilu1.mat -ascii;   %支路數(shù)據(jù)文件                   
load fuhe1.mat -ascii;   %負(fù)荷數(shù)據(jù)文件                   
load fadian1.mat -ascii;  %發(fā)電機(jī)數(shù)據(jù)文件
load gpm1.mat -ascii;  %發(fā)電機(jī)參數(shù)文件
%load fhtx.mat -ascii; %負(fù)荷靜態(tài)特性數(shù)組文件                   
%-----------------------------------------
nm=max(zhilu1(:,2:3)); nn=max(nm');%網(wǎng)絡(luò)節(jié)點(diǎn)總數(shù)
bnn=size(zhilu1);      bn=bnn(1,1);%支路數(shù)
lnn=size(fuhe1);      ln=lnn(1,1);%負(fù)荷節(jié)點(diǎn)數(shù)
gnn=size(fadian1);     gn=gnn(1,1);%發(fā)電機(jī)節(jié)點(diǎn)數(shù)
%-----------------------------------------
%發(fā)電機(jī)的類型參數(shù)
mt1=[4 4 4]; 
%==================BX型PQ分解法=====================
%形成bp,不考慮充電電容及變壓器非標(biāo)準(zhǔn)變比對(duì)導(dǎo)納矩陣的影響
%但要考慮支路電阻的影響
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,不考慮支路電阻的影響,
%但要考慮支路充電電容及變壓器的非標(biāo)準(zhǔn)變比的影響
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   %非標(biāo)準(zhǔn)變比側(cè)          
         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   %非標(biāo)準(zhǔn)變比側(cè)
         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);
%---------------------------------------
%求取導(dǎo)納矩陣
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   %非標(biāo)準(zhǔn)變比側(cè)          
         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   %非標(biāo)準(zhǔn)變比側(cè)
         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節(jié)點(diǎn)的個(gè)數(shù)
v=ones(nn,1);  %電壓
pg=zeros(nn,1);%發(fā)電機(jī)有功
qg=zeros(nn,1);%發(fā)電機(jī)無功
pl=zeros(nn,1);%負(fù)荷有功
ql=zeros(nn,1);%負(fù)荷無功
o=zeros(nn,1); %電壓相角
%----------------------------------
%根據(jù)節(jié)點(diǎn)類型對(duì)bp，bpp進(jìn)行修正原則是:將bp，bpp通取為n階
%PO迭代:bp中對(duì)平衡節(jié)點(diǎn)所在的行和列,除對(duì)角線元素取作1外,其他元素均取為0.
%其相應(yīng)的右端元素opv(i)也取為0(見后面計(jì)算);
%QV迭代:bpp中對(duì)平衡節(jié)點(diǎn)和PV節(jié)點(diǎn)所在的行和列,除對(duì)角線元素取作1外,其他元素均取為0.
%其相應(yīng)的右端元素oqv(i)也取為0(見后面計(jì)算).
%取PQ節(jié)點(diǎn)的有功和無功功率和PV節(jié)點(diǎn)的有功,電壓
for l=1:gnn;%發(fā)電機(jī)
   if fadian1(l,2)==0&fadian1(l,3)==0 %平衡節(jié)點(diǎn)
      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節(jié)點(diǎn)
      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節(jié)點(diǎn)的節(jié)點(diǎn)號(hào)
   else
      f=fadian1(l,1);
      pg(f)=fadian1(l,2);
      qg(f)=fadian1(l,3);
   end
end
for l=1:lnn %負(fù)荷
   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);%計(jì)算常數(shù)項(xiàng)
       av=av+v(i);
     end
     opv(lxl)=0;%位于后面計(jì)算方便,令其維數(shù)為n
     av=(av+1.05)/nn;%平均電壓
     oo=inv(bp)*opv/av;%計(jì)算誤差量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);%計(jì)算常數(shù)項(xiàng)
    end
    for i=1:pv
     l=pvj(i);
     oqv(l)=0;%對(duì)PV節(jié)點(diǎn)令其右端項(xiàng)為0
    end
   oqv(lxl)=0;
   ov=inv(bpp)*oqv;
   v=v-ov;%修正幅值
 end 
 %----------下一次迭代-------- 
end
%=================潮流已經(jīng)收斂======================2002.3.28
%計(jì)算平衡節(jié)點(diǎn)的有功,無功
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
%計(jì)算PV節(jié)點(diǎn)的無功功率
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
%================潮流計(jì)算完畢====================================
%與特征值計(jì)算程序銜接2002.3.28
%支路數(shù)據(jù)文件的銜接
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
%負(fù)荷數(shù)據(jù)的銜接
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
%發(fā)電機(jī)節(jié)點(diǎn)數(shù)據(jù)的銜接
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 %平衡節(jié)點(diǎn)
      gpqv(i,2)=1;
   end
   for l=1:pv
      if pvj(l)==abs(fadian1(i,1));%PV節(jié)點(diǎn)
         gpqv(i,2)=-1;
      end
   end
end
%=========================銜接部分完成===========================
%--------------------------------------------------------------------------------------------------
ef=[1,1,1];   %各發(fā)電機(jī)勵(lì)磁系統(tǒng)的階數(shù)
%統(tǒng)計(jì)發(fā)電機(jī)的各個(gè)參數(shù)的個(gè)數(shù)及各發(fā)電機(jī)的階數(shù)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
%---------------------------------------------------------------------------------------------------
%====================================
%原始網(wǎng)絡(luò)導(dǎo)納矩陣yy
%====================================
%恒定阻抗負(fù)荷模型下的網(wǎng)絡(luò)導(dǎo)納矩陣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;
%--------------------------------
%消去中間節(jié)點(diǎn)后的發(fā)電機(jī)節(jié)點(diǎn)間輸電導(dǎo)納矩陣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; %發(fā)電機(jī)節(jié)點(diǎn)（含無窮大母線）間的導(dǎo)納矩陣
%--------------------------------

if mt1(1,1)==0            %有無窮大母線時(shí)
   yt=yt1(2:gn,2:gn);     %發(fā)電機(jī)節(jié)點(diǎn)間導(dǎo)納矩陣
   yt0=yt1(2:gn,1);       %無窮大母線與發(fā)電機(jī)節(jié)點(diǎn)間的互導(dǎo)納向量
   u0=gpqv(1,7);          %無窮大母線電壓幅值
   gpqv(1,:)=[];
   mt=mt1(1,2:gn);