潮流样本.m

来自「潮流样本。能运用于电力系统潮流计算」· M 代码 · 共 658 行 · 第 1/2 页

   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;