ljxsh61183xg1.m

MATLAB潮流计算

      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;xq=xq1;
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绕组,
      tqp=0;tqpp=0;tdpp=0;
      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
   %pwv(:,k)=pwv((:,k)/pwvm(1,k);
end
pwvm=pwvm';
for k=1:n
   pwvm(k,3)=eig0(k);
end
%找出和特征值次相关的变量位置及其相关度
%注意此时变量是以变量为一簇的
pwvm(:,5)=(min(pwv))';
for i=1:n
   for l=1:n
      wox=pwv(l,i);
      if wox<pwvm(i,1)&wox>pwvm(i,5);
         pwvm(i,5)=wox;
         pwvm(i,4)=l;
      end
   end
end

%=====================================
%机电模式相关比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;%（1特征值编号、2特征值、3阻尼比、4频率、5强相关度、
          %6强相关变量位置、7次相关度、8次相关变量位置）
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:9)=pwvm;
%====================-------------------------------======================
format short

%==========k1-k18模型多机系统小扰动稳定的解耦分析x_========================
load 61183xg1w.mat x_wr;
x_w=100*pi;
one=eye(gn,gn);
for l=1:gn  %低频振荡频率0.2-2.5
   x_x=j*x_wr(l)*2*pi; %初值
   %各控制器的传递函数
   Ge=ke/(one+te*x_x);%励磁系统
   Ggov=0;%原动机系统
   Gpss=0;%PSS系统
   %计算简化二阶的系统阵
   Gd=inv(one+x_x*tqpp);
   Gq=inv(one+x_x*tdpp);
   %
   Gtd=inv((one+x_x*tqp)*Gd);%
   Gtq=inv((one+x_x*tdp)*Gq);%
   %Gdpp=inv(one+Gd*k11)*Gd;
   Gdpp=inv(one+Gd*(k11+Gtd*k17))*Gd;%
   Gqpp=inv(one+Gq*(k9+Gtq*(k3+Ge*k6)))*Gq;%
   %Gqd=k15;
   Gqd=k15+Gtd*k16;%
   Gdq=k14+Gtq*(k13+Ge*k7);%
   %Gdo=k10;
   Gdo=k10+Gtd*k18;%
   Gqo=k8+Gtq*(k4+Ge*k5);%
   Gw=Gpss*Ge*Gtq;
   L_s=inv(one-Gdpp*Gqd*Gqpp*Gdq);
   L_d=D+Ggov+k2*L_s*Gqpp*Gw-k12*L_s*Gdpp*Gqd*Gqpp*Gw;
   L_o=k1+k2*L_s*(Gqpp*Gdq*Gdpp*Gdo-Gqpp*Gqo)+k12*L_s*(Gdpp*Gqd*Gqpp*Gqo-Gdpp*Gdo);

   x_sw=eye(gn)*x_x/x_w;%s*inv(wn)对角阵
   x_h=L_o+L_d*x_sw;
   x_xI=x_x*eye(gn);%初值形成对角阵
   x_a=(x_xI+inv(tj)*L_d)*x_sw+inv(tj)*L_o;
   x_a1=x_a;
   x_a1(l,:)=[];
   x_a1(:,l)=[];%形成A杠
   x_ai=x_a(:,l);
   x_ai(l)=[];%形成Ai-i
   x_k=-inv(x_a1)*x_ai;%koi 是列向量
   x_hi=x_h(l,:);
   x_hi(l)=[];%形成Hi-i 行相量
   x_hk=x_hi'.*x_k;%在fl下各机对第l机的互数--列向量
      
   x_sij(:,l)=real(x_hk);%互同步力矩系数 在fl下所有机对第l机的互数
   x_dij(:,l)=imag(x_hk)*x_w/x_x*j;%互阻尼力矩系数 在fl下所有机对第l机的互数
   x_sii(l,l)=real(x_h(l,l));%自同步力矩系数
   x_dii(l,l)=imag(x_h(l,l))*x_w/x_x*j;%自阻尼力矩系数
end
%计算力矩系数矩阵
for k=1:gn
   h=1;
   while h<gn+1
      if h==k %取自力矩系数
         x_si(h,k)=x_sii(h,k);
         x_di(h,k)=x_dii(h,k);
      elseif h>k %此时，互力矩矩阵元素行比力矩系数矩阵少1
         x_si(h,k)=x_sij(h-1,k);
         x_di(h,k)=x_dij(h-1,k);
      else   %h<k，此时可以取相应的元素
         x_si(h,k)=x_sij(h,k);
         x_di(h,k)=x_dij(h,k);
      end
      h=h+1;
   end
end
x_si(gn+1,:)=sum(x_sij);%各机总的互力矩系数
x_di(gn+1,:)=sum(x_dij);%各机总的互力矩系数
x_si(gn+2,:)=x_si(gn+1,:)+diag(x_sii)';%各机综合力矩系数
x_di(gn+2,:)=x_di(gn+1,:)+diag(x_dii)';%各机综合力矩系数