xg61183.m

6机系统k1-k18模型控制理论运用。

%6机系统k1-k18模型控制理论运用
%求取x_g，1+x_g，稳定裕量
%调入数据文件(电科院6机系统)
clear
load zhilu1.mat -ascii;   %支路数据文件                   
load fuhe1.mat -ascii;   %负荷数据文件                   
load fadian1.mat -ascii;  %发电机数据文件
load gpm1.mat -ascii;  %发电机参数文件
%load fhtx.mat -ascii; %负荷静态特性数组文件                   
%-----------------------------------------
%发电机的类型参数
mt1=[5 5 5 5 5 5];
%mt1=[4 4 4 4 4 4];
%mt1=[1 1 1 1 1 1];
%-----------------------------------------
nm=max(zhilu1(:,2:3)); nn=max(nm');%网络节点总数
bnn=size(zhilu1);      bn=bnn(1,1);%支路数
lnn=size(fuhe1);      ln=lnn(1,1);%负荷节点数
gnn=size(fadian1);     gn=gnn(1,1);%发电机节点数
%==================BX型PQ分解法=====================
%形成bp,不考虑充电电容及变压器非标准变比对导纳矩阵的影响
%但要考虑支路电阻的影响
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,不考虑支路电阻的影响,
%但要考虑支路充电电容及变压器的非标准变比的影响
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   %非标准变比侧          
         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   %非标准变比侧
         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);
%---------------------------------------
%求取导纳矩阵
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   %非标准变比侧          
         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   %非标准变比侧
         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节点的个数
v=ones(nn,1);  %电压
pg=zeros(nn,1);%发电机有功
qg=zeros(nn,1);%发电机无功
pl=zeros(nn,1);%负荷有功
ql=zeros(nn,1);%负荷无功
o=zeros(nn,1); %电压相角
%----------------------------------
%根据节点类型对bp，bpp进行修正原则是:将bp，bpp通取为n阶
%PO迭代:bp中对平衡节点所在的行和列,除对角线元素取作1外,其他元素均取为0.
%其相应的右端元素opv(i)也取为0(见后面计算);
%QV迭代:bpp中对平衡节点和PV节点所在的行和列,除对角线元素取作1外,其他元素均取为0.
%其相应的右端元素oqv(i)也取为0(见后面计算).
%取PQ节点的有功和无功功率和PV节点的有功,电压
for l=1:gnn;%发电机
   if fadian1(l,2)==0&fadian1(l,3)==0 %平衡节点
      f=fadian1(l,1);
      lxl=f;
      v(f)=1.00;
      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节点
      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节点的节点号
   else
      f=fadian1(l,1);
      pg(f)=fadian1(l,2);
      qg(f)=fadian1(l,3);
   end
end
for l=1:lnn %负荷
   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);
%--------------------------------------
while(max(abs(opv))>0.000001|max(abs(oqv))>0.000001)
   %---------PO迭代----------------
   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);%计算常数项
       av=av+v(i);
     end
     opv(lxl)=0;%位于后面计算方便,令其维数为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);%计算常数项
    end
    for i=1:pv
     l=pvj(i);
     oqv(l)=0;%对PV节点令其右端项为0
    end
   oqv(lxl)=0;
   ov=inv(bpp)*oqv;
   v=v-ov;%修正幅值
 end 
 %----------下一次迭代-------- 
end
%=================潮流已经收敛======================2002.3.28
%计算平衡节点的有功,无功
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节点的无功功率
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
%支路数据文件的衔接
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
%负荷数据的衔接
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
%发电机节点数据的衔接
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 %平衡节点
      gpqv(i,2)=1;
   end
   for l=1:pv
      if pvj(l)==abs(fadian1(i,1));%PV节点
         gpqv(i,2)=-1;
      end
   end
end
%=========================衔接部分完成===========================

%原始网络导纳矩阵yy
%================================
%恒定阻抗负荷模型下的网络导纳矩阵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;
%--------------------------------
%消去中间节点后的发电机节点间输电导纳矩阵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; %发电机节点（含无穷大母线）间的导纳矩阵
%--------------------------------
if mt1(1,1)==0            %有无穷大母线时
   yt=yt1(2:gn,2:gn);     %发电机节点间导纳矩阵
   yt0=yt1(2:gn,1);       %无穷大母线与发电机节点间的互导纳向量
   u0=gpqv(1,7);          %无穷大母线电压幅值
   gpqv(1,:)=[];
   mt=mt1(1,2:gn);
   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));
%-------------------------------
%根据发电机模型类型修正发电机参数
xq1=xq;
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;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;
%====================================
%==========k1-k18模型多机系统小扰动稳定的解耦分析x_========================
x_w=100*pi;
one=eye(gn);
m=0;
for l=0.01:0.001:10%频率
    x_x=j*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;
   for k=1:gn
      x_a1=x_a;
      x_a1(k,:)=[];
      x_a1(:,k)=[];%形成A杠
      x_ai=x_a(:,k);
      x_ai(k)=[];%形成Ai-i
      x_k(:,k)=-inv(x_a1)*x_ai;%koi是列向量
   end
   for k=1:gn
      x_hi=x_h(k,:);
      x_hi(k)=[];%形成Hi-i
      x_hk(k,1)=x_hi*x_k(:,k);%单机是个元素，所有机是列向量
      x_h1(k,k)=x_h(k,k)+x_hk(k,1);%单机是个元素，所有机是对角矩阵H杠
      
   end
   %计算开环传递函数 
   m=m+1;
   x_g0=x_w*x_h1/(tj*x_x*x_x);
   x_g(:,m)=diag(x_g0);
end
%1+x_g
[n,m]=size(x_g);
for k=1:n
   for l=1:m
      x_g1(k,l)=1+x_g(k,l);
   end
end
%稳定裕量
x_m=0.01:0.001:10;
for k=1:n
   [gm(k),pm(k),wg(k),wp(k)]=margin(abs(x_g(k,:)),angle(x_g(k,:))*180/pi,x_m);
   gm(k)=20*log10(gm(k));
end