gn=gn-1;               %實際發電機節點數
 else                     %無無窮大母線時
   yt=yt1;                %發電機節點間導納矩陣
   yt0=zeros(gn,1);       %無窮大母線與發電機節點間的互導納置0。
   u0=0;
   mt=mt1;
end
%====================================
%發電機參數對角矩陣
xd=diag(gpm1(:,2));    xdp=diag(gpm1(:,3));    xdpp=diag(gpm1(:,4));
xq=diag(gpm1(:,5));    xqpp=diag(gpm1(:,6));   ra=diag(gpm1(:,7));
tdp=diag(gpm1(:,8));   tdpp=diag(gpm1(:,9));   tqpp=diag(gpm1(:,10));
tj=diag(gpm1(:,11));   ke=diag(gpm1(:,12));    te=diag(gpm1(:,13));
D=diag(gpm1(:,14));
%采用六階模型,加入
xqp=diag(gpm1(:,15));  tqp=diag(gpm1(:,16));
%-------------------------------
%根據發電機模型類型修正發電機參數
xqq=zeros(gn,gn);   xqqp=zeros(gn,gn);
for n1=1:gn,
   if mt(1,n1)==1,         %五階模型考慮阻尼繞組D+Q
      D(n1,n1)=0;
   end
   if mt(1,n1)==2,         %只考慮交軸阻尼繞組Q 
      xdpp(n1,n1)=xdp(n1,n1);
      D(n1,n1)=0;
   elseif mt(1,n1)==3,     %只考慮直軸d阻尼繞組D
      xqpp(n1,n1)=xq(n1,n1);
      D(n1,n1)=0;
   elseif mt(1,n1)==4,     %三階模型考慮綜合阻尼繞組
      xdpp(n1,n1)=xdp(n1,n1);
      xqpp(n1,n1)=xq(n1,n1);
   elseif mt(1,n1)==5,     %六階模型考慮阻尼繞組D+g+Q
      D(n1,n1)=0;
      xqq(n1,n1)=xq(n1,n1)-xqp(n1,n1);
      xqqp(n1,n1)=xqp(n1,n1)-xqpp(n1,n1);
      xq(n1,n1)=xqpp(n1,n1);
   end
end
%--------------------------------
xdd=xd-xdp;           xddp=xdp-xdpp;         xqqpp=xq-xqpp+xqqp;
xqd=xqpp-xdpp;        z=ra+j*xdpp;           
%====================================
%D-Q坐標系與d-q坐標系間的變換矩陣T 
PQ=gpqv(:,3)+j*gpqv(:,4); %發電機復功率
U=gpqv(:,9)+j*gpqv(:,10); %發電機機端電壓 I=conj(PQ./U);                      %發電機定子電流
EQ=U+(ra+j*xq)*I;                   %虛構電勢
d0=angle(EQ)-pi/2;                  %d-q坐標系超前D-Q坐標系的角度
T=diag(exp(j*d0));                  % D-Q坐標系與d-q坐標系間的變換矩陣T 
%====================================
%變量初值計算
i=inv(T)*I;            id=diag(real(i));     iq=diag(imag(i));
u=inv(T)*U;            ud=diag(real(u));     uq=diag(imag(u)); 
ut=diag(abs(u));       epp=u+z*i-xqd*imag(i);
edp=diag(real(epp));   eqp=diag(imag(epp));
%====================================
%系數矩陣k1-k15 （gn*gn階）
y=inv(T)*inv(inv(yt)+z)*T;             yr=real(y);   yi=imag(y);
y0=diag(inv((inv(yt)+z))*inv(yt)*yt0); y0r=real(y0);y0i=imag(y0);
md=-yi*edp-yr*eqp-yi*xqd*iq;           mq=yr*edp-yi*eqp+yr*xqd*iq;
pd=md-diag(sum(md'))+y0i*u0;           pq=mq-diag(sum(mq'))-y0r*u0;
yq=inv(eye(gn)-yi*xqd)*yr;             yd=-yi+yr*xqd*yq;
fq=inv(eye(gn)-yi*xqd)*pq;             fd=pd+yr*xqd*fq;
kq=inv(eye(gn)-yi*xqd)*yi;             kd=yr+yr*xqd*kq;
bd=edp+xqd*iq;                         bq=eqp+xqd*id;
cd=ud*inv(ut)*ra+uq*inv(ut)*xdpp;      cq=ud*inv(ut)*xqpp-uq*inv(ut)*ra;
%----------------------------------
k1=bd*fd+bq*fq;  k2=bd*yd+bq*yq+iq;    k3=xdd*yd;         k4=xdd*fd;
k5=cq*fq-cd*fd;  k6=cq*yq-cd*yd+uq;    k7=cq*kd-cd*kq+ud; k8=xddp*fd;
k9=xddp*yd;      k12=bd*kd+bq*kq+id;
k13=xdd*kd;      k14=xddp*kd;         
%六階模型需計算---------------------
k10=-xqqpp*fq;    k11=-xqqpp*kq;         k15=-xqqpp*yq;       
k18=-xqq*fq;      k17=-xqq*kq;           k16=-xqq*yq;
%====================================
%六階模型需要加入狀態變量edp
%狀態向量x=[δ,ω,eqp,eqpp,edp,edpp,efd]'7*gn維
a=zeros(7*gn);

a(1:gn,gn+1:2*gn)=100*pi*eye(gn);

a(gn+1:2*gn,1:gn)=-inv(tj)*k1;
a(gn+1:2*gn,gn+1:2*gn)=-inv(tj)*D;
a(gn+1:2*gn,3*gn+1:4*gn)=-inv(tj)*k2;
a(gn+1:2*gn,5*gn+1:6*gn)=-inv(tj)*k12;

a(2*gn+1:3*gn,1:gn)=-inv(tdp)*k4;
a(2*gn+1:3*gn,2*gn+1:3*gn)=-inv(tdp);
a(2*gn+1:3*gn,3*gn+1:4*gn)=-inv(tdp)*k3;
a(2*gn+1:3*gn,5*gn+1:6*gn)=-inv(tdp)*k13;
a(2*gn+1:3*gn,6*gn+1:7*gn)=inv(tdp);

a(3*gn+1:4*gn,1:gn)=-inv(tdpp)*k8-inv(tdp)*k4;
a(3*gn+1:4*gn,2*gn+1:3*gn)=inv(tdpp)-inv(tdp);
a(3*gn+1:4*gn,3*gn+1:4*gn)=-inv(tdpp)-inv(tdpp)*k9-inv(tdp)*k3;
a(3*gn+1:4*gn,5*gn+1:6*gn)=-inv(tdpp)*k14-inv(tdp)*k13;
a(3*gn+1:4*gn,6*gn+1:7*gn)=inv(tdp);
%處理tqp的不可逆的情況                                         
for n1=1:gn,
   if mt(1,n1)<5,
      tqp(n1,n1)=1;
   end
end   
vtqp=inv(tqp);
for n1=1:gn,
   if mt(1,n1)<5,
      vtqp(n1,n1)=0;
   end
end   
a(4*gn+1:5*gn,1:gn)=-vtqp*k18;
a(4*gn+1:5*gn,3*gn+1:4*gn)=-vtqp*k16;
a(4*gn+1:5*gn,4*gn+1:5*gn)=-vtqp;
a(4*gn+1:5*gn,5*gn+1:6*gn)=-vtqp*k17;


a(5*gn+1:6*gn,1:gn)=-inv(tqpp)*k10-vtqp*k18;
a(5*gn+1:6*gn,3*gn+1:4*gn)=-inv(tqpp)*k15-vtqp*k16;
a(5*gn+1:6*gn,4*gn+1:5*gn)=inv(tqpp)-vtqp;
a(5*gn+1:6*gn,5*gn+1:6*gn)=-inv(tqpp)-inv(tqpp)*k11-vtqp*k17;


a(6*gn+1:7*gn,1:gn)=-inv(te)*ke*k5;
a(6*gn+1:7*gn,3*gn+1:4*gn)=-inv(te)*ke*k6;
a(6*gn+1:7*gn,5*gn+1:6*gn)=-inv(te)*ke*k7;
a(6*gn+1:7*gn,6*gn+1:7*gn)=-inv(te);
%--------------------------------
%根據發電機模型類型修正A 陣
for n=1:gn,
   if mt(1,n)==1,                         %無g繞組,
      a(5*gn+n,:)=a(5*gn+n,:)-a(4*gn+n,:);%edpp行-edp行
      a(4*gn+n,:)=0;                      %edp行置0
      a(:,4*gn+n)=0;                      %edp列置0
   elseif mt(1,n)==2,                     %無D+g繞組，
      a(3*gn+n,:)=0;                      %eqpp 行置0，
      a(:,2*gn+n)=a(:,2*gn+n)+a(:,3*gn+n);% eqpp 列+ eqp列
      a(:,3*gn+n)=0;                      %eqpp列置0
      
      a(5*gn+n,:)=a(5*gn+n,:)-a(4*gn+n,:);%edpp行-edp行
      a(4*gn+n,:)=0;                      %edp行置0
      a(:,4*gn+n)=0;                      %edp列置0

   elseif mt(1,n)==3                      %無g+Q繞組
      a(5*gn+n,:)=a(5*gn+n,:)-a(4*gn+n,:);%edpp行-edp行
      a(4*gn+n,:)=0;                      %edp行置0
      a(:,4*gn+n)=0;                      %edp列置0

      a(5*gn+n,:)=0;                      %edpp 行置0，
      a(:,5*gn+n)=0;                      %edpp列置0
   elseif mt(1,n)==4,                     %無D+g+Q繞組，
      a(3*gn+n,:)=0;                      %eqpp 行置0，
      a(:,2*gn+n)=a(:,2*gn+n)+a(:,3*gn+n);% eqpp 列+ eqp列
      a(:,3*gn+n)=0;                      %eqpp列置0
      
      a(5*gn+n,:)=a(5*gn+n,:)-a(4*gn+n,:);%edpp行-edp行
      a(4*gn+n,:)=0;                      %edp行置0
      a(:,4*gn+n)=0;                      %edp列置0

      a(5*gn+n,:)=0;                      %edpp 行置0，
      a(:,5*gn+n)=0;                      %edpp列置0
   end
end
%----------------------------------
%根據發電機模型類型壓縮A 陣
dp=0;qp=0;gp=0;
for n=1:gn,
   k=mt(1,n);
   if k==1,                                      % 無g繞組
      a(4*gn+n-gp-dp,:)=[];a(:,4*gn+n-gp-dp)=[]; % 刪edp 行和列
      a(7*gn,7*gn)=1;       gp=gp+1;             %末元素置1
   end   
   if k==2,                                      % 無D+g繞組
      a(3*gn+n-dp,:)=[];    a(:,3*gn+n-dp)=[];   % 刪eqpp 行和列
      a(7*gn,7*gn)=1;       dp=dp+1;             %末元素置1
      
      a(4*gn+n-gp-dp,:)=[];a(:,4*gn+n-gp-dp)=[]; % 刪edp 行和列
      a(7*gn,7*gn)=1;       gp=gp+1;             %末元素置1
   end
   if k==3,                                      % 無g+Q繞組
      a(4*gn+n-gp-dp,:)=[];a(:,4*gn+n-gp-dp)=[]; % 刪edp 行和列
      a(7*gn,7*gn)=1;       gp=gp+1;             %末元素置1

      a(5*gn+n-qp-dp-gp,:)=[]; a(:,5*gn+n-qp-dp-gp)=[];%刪edpp 行和列
      a(7*gn,7*gn)=1;       qp=qp+1;             %末元素置1
   end
   if k==4,                                      % 無D+g+Q繞組,
      a(3*gn+n-dp,:)=[];    a(:,3*gn+n-dp)=[];   % 刪eqpp 行和列
      a(7*gn,7*gn)=1;       dp=dp+1;             %末元素置1
      
      a(4*gn+n-gp-dp,:)=[];a(:,4*gn+n-gp-dp)=[]; % 刪edp 行和列
      a(7*gn,7*gn)=1;       gp=gp+1;             %末元素置1

      a(5*gn+n-qp-dp-gp,:)=[]; a(:,5*gn+n-qp-dp-gp)=[];%刪edpp 行和列
      a(7*gn,7*gn)=1;       qp=qp+1;             %末元素置1
   end
end
dq=dp+qp+gp;                                     %刪去的行列總數
a(7*gn-dq+1:7*gn,:)=[];                          %刪去最后dq行和列
a(:,7*gn-dq+1:7*gn)=[];
%====================================
%計算特征根,自然振蕩頻率,阻尼比eig0,wn,k
[v1,d]=eig(a);                  %左特征向量，特征根對角陣
eig0=diag(d);                   %特征根
%====================================
%計算右特征向量w和參與矩陣pwv
w=inv(v1)';                     %右特征向量
pwv=abs(w.*v1);                 %參與矩陣
%====================================
n=7*gn-dq;                     %以七階模型為基礎
pwvm=0.;
for k=1:n
   pwvm(1,k)=max(pwv(:,k));
   for kk=1:n
      if pwvm(1,k)==pwv(kk,k);
         pwvm(2,k)=kk;
      end
   end
end
pwvm=pwvm';
%=====================================
%機電模式相關比pmm和電氣模式相關比pee
n=7*gn-dq;                     %以七階模型為基礎
pmm=(sum(pwv(1:2*gn,:))./sum(pwv(2*gn+1:n,:)))';
format long,
pee=1./pmm;
%====================================
%pmm(n-2:n,1)=pmm(n-2:n,1)*0;
%=====================================
ewk=0;
for n=1:7*gn-dq                %以七階模型為基礎
   ewk(n,1)=n;                 %特征根序號
end
ewk(:,2)=eig0;                  %特征?
ewk(:,3)=-real(eig0)./abs(eig0); %特征根阻謀?
ewk(:,4)=abs(imag(eig0)/(2*3.1416));
ewk(:,5)=ewk(:,1);
ewk(:,6:7)=pwvm;
size(ewk);

%====================================
%模態分布|vki/vkj|
k=7*gn-dq;                 %以七階模型為基礎
for n=1:gn
   vv(n,1:k)=v1(n,1:k)./v1(1,1:k); % 以首行元素（1#機功角）為參考的模態分布矩陣
end
vvabs=abs(vv);             % 模態振幅相對值
vvangle=angle(vv)*180/pi;  % 模態相角相對于1#機的值
%===================================
%特征根對阻尼系數的靈敏度esD
for m=1:7*gn-dq            %以七階模型為基礎
   for n=1:gn
      as=zeros(7*gn-dq);   %以七階模型為基礎
      as(gn+n,gn+n)=1./tj(n,n);
      es=(w(:,m)'*as*v1(:,m));
      esD((m-1)*gn+n,1)=m;
      esD((m-1)*gn+n,2)=n;
      esD((m-1)*gn+n,3)=abs(es);
      esD((m-1)*gn+n,4)=angle(es)*180/pi;
   end
end
%========================================
%特征根對AVR增益ke的靈敏度eske
for m=1:7*gn-dq            %以七階模型為基礎
   for n=1:gn
      as=zeros(7*gn-dq);   %以七階模型為基礎
      as(5*gn+n-dq,1:gn)=-1./tj(n,n)*k5(n,:);
      as(5*gn+n-dq,3*gn+1-dq:4*gn-dq)=-1./tj(n,n)*k6(n,:);
      as(5*gn+n-dq,4*gn+1-dq:5*gn-dq)=-1./tj(n,n)*k7(n,:);
      es=(w(:,m)'*as*v1(:,m));
      eske((m-1)*gn+n,1)=m;
      eske((m-1)*gn+n,2)=n;
      eske((m-1)*gn+n,3)=abs(es);
      eske((m-1)*gn+n,4)=angle(es)*180/pi;
   end
end
%========================================
%特征根對AVR時間常數te的靈敏度este 
for m=1:7*gn-dq            %以七階模型為基礎
   for n=1:gn
      as=zeros(7*gn-dq);   %以七階模型為基礎
      as(5*gn+n-dq,1:gn)=(1./tj(n,n))^2*ke(n,n)*k5(n,:);
      as(5*gn+n-dq,3*gn+1-dq:4*gn-dq)=(1./tj(n,n))^2*ke(n,n)*k6(n,:);
      as(5*gn+n-dq,4*gn+1-dq:5*gn-dq)=(1./tj(n,n))^2*ke(n,n)*k7(n,:);
      as(5*gn+n-dq,5*gn+1-dq:6*gn-dq)=(1./tj(n,n))^2;
      es=(w(:,m)'*as*v1(:,m));
      este((m-1)*gn+n,1)=m;
      este((m-1)*gn+n,2)=n;
      este((m-1)*gn+n,3)=abs(es);
      este((m-1)*gn+n,4)=angle(es)*180/pi;
   end
end
%==========================================
%計算結果輸出
format short
%disp('發電機模型類型，阻尼系數')
%mt,D,
%disp('特征根序號，特征根, 機電相關比，阻尼比')
%disp('自然頻率，阻尼頻率，最大相關因子，強相關狀態變量序號')
%pmm,vtqp,ewk,

%disp('模態幅值分布，模態相角分布')
%vvabs
%vvangle
%disp('線性化系數矩陣')
%k1,k2,k3,k4,k5,k6,k7,k8
%k9,k10,k11,k12,k13,k14,k15
%pwv(:,1:6),pwv(:,7:12),pwv(:,13:18),
%pwv(:,19:24)%,pwv(:,25:30),pwv(:,31:36)
%format short
%bzm1;,
%lpqv1;,
%gpqv1;,
%gpm1;,
%nn;,
%gn;,
%ln;,
%yy;,
%ynl;,
%yt;,
%a;,
%v1;,
%w;,
%pee;,
%esD;,
%eske;,
%este;