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



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

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