fm_eigen.m

这是一个很适合研究和学习用的电力系统仿真软件

    Jlfvr = (DAE.Jlfv(Bus.n+1:2*Bus.n,Bus.n+1:2*Bus.n) ...
             - DAE.Jlfv(Bus.n+1:2*Bus.n, 1:Bus.n)*...
      (DAE.Jlfv(1:Bus.n,1:Bus.n)\DAE.Jlfv(1:Bus.n,Bus.n+1:2* ...
                                          Bus.n)));
    Jlfvr(SW.bus,SW.bus) = Jlfvr(SW.bus,SW.bus) + 998* ...
        eye(SW.n);
    Jlfvr(PV.bus,PV.bus) = Jlfvr(PV.bus,PV.bus) + 998* ...
        eye(PV.n);
    if opts
      auto = eigs(Jlfvr,neig,meth{opts},opt);
    else
      [V, auto] = eig(full(Jlfvr));
      auto = diag(auto);
    end
    auto = round(auto/Settings.lftol)*Settings.lftol;
    num_auto = length(auto);
    autor = real(auto);
    autoi = imag(auto);
    names = cellstr(strcat('Eig Jlfv #',num2str([1:num_auto]')));

    if ~opts
      W = inv(V);
      WtV = sum(abs(W').*abs(V));
      pf = abs(W).*abs(V');
      for i = 1:length(auto), pf(i,:) = pf(i,:)/WtV(i); end
      %pf = abs(vect'.*(inv(vect'))');
    end
    Header{2,1} = 'EIGENVALUES OF THE COMPLETE POWER JACOBIAN MATRIX';
    Cols{2,1} = {'Eigevalue', 'Real part', ...
                   'Imaginary Part'};
    Rows{2,1} = names;
    Matrix{2,1} = [autor, autoi];

    if ~opts
      uno = fix(Bus.n/5);
      due = rem(Bus.n,5);
      if due > 0, uno = uno + 1; end
      for k = 1:uno
        Header{2+k,1} = 'PARTICIPATION FACTORS (Euclidean norm)';
        Cols{2+k} = {'  ',Varname.bus{5*(k-1)+1:min(5*k,Bus.n)}};
          Rows{2+k} = names;
          Matrix{2+k,1} = pf(:,5*(k-1)+1:min(5*k,Bus.n));
      end
    end

  elseif tipo_mat ==  1

    if DAE.n == 0
      fm_disp('Since no dynamic component is loaded, Jlfd = Jlfv.',2)
      DAE.Jlfd = DAE.Jlfv;
    else
      Fx_mod = DAE.Fx+diag(-1e-5*ones(DAE.n,1));
      DAE.Jlfd = (DAE.Jlfv - DAE.Gx*(Fx_mod\DAE.Fy));
    end
    Jlfdr = (DAE.Jlfd(Bus.n+1:2*Bus.n,Bus.n+1:2*Bus.n) ...
             - DAE.Jlfd(Bus.n+1:2*Bus.n, 1:Bus.n)*...
      (DAE.Jlfd(1:Bus.n,1:Bus.n)\DAE.Jlfd(1:Bus.n,Bus.n+1:2* ...
                                          Bus.n)));
    Jlfdr(SW.bus,SW.bus) = Jlfdr(SW.bus,SW.bus) + 998* ...
        speye(SW.n);
    Jlfdr(PV.bus,PV.bus) = Jlfdr(PV.bus,PV.bus) + 998* ...
        speye(PV.n);

    if opts
      auto = eigs(Jlfdr,neig,meth{opts},opt);
    else
      [V, auto] = eig(full(Jlfdr));
      auto = diag(auto);
    end
    auto = round(auto/Settings.lftol)*Settings.lftol;
    num_auto = length(auto);
    autor = real(auto);
    autoi = imag(auto);
    names = cellstr(strcat('Eig Jlfd #',num2str([1:num_auto]')));

    if ~opts
      W = inv(V);
      WtV = sum(abs(W').*abs(V));
      pf = abs(W).*abs(V');
      for i = 1:length(auto), pf(i,:) = pf(i,:)/WtV(i); end
    end
    Header{2,1} = 'EIGENVALUES OF THE DYNAMIC POWER JACOBIAN MATRIX';
    Cols{2,1} = {'Eigevalue', 'Real part', ...
                 'Imaginary Part'};
    Rows{2,1} = names;
    Matrix{2,1} = [autor, autoi];

    if ~opts
      uno = fix(Bus.n/5);
      due = rem(Bus.n,5);
      if due > 0, uno = uno + 1; end
      for k = 1:uno
        Header{2+k,1} = 'PARTICIPATION FACTORS (Euclidean norm)';
        Cols{2+k} = {'  ',Varname.bus{5*(k-1)+1:min(5*k,Bus.n)}};
        Rows{2+k} = names;
        Matrix{2+k,1} = pf(:,5*(k-1)+1:min(5*k,Bus.n));
      end
    end

  end

  auto_neg  = find(autor < 0);
  auto_pos  = find(autor > 0);
  auto_real = find(autoi == 0);
  auto_comp = find(autoi < 0);
  auto_zero = find(autor == 0);

  num_neg  = length(auto_neg);
  num_pos  = length(auto_pos);
  num_real = length(auto_real);
  num_comp=length(auto_comp);
  num_zero = length(auto_zero);

  if Fig.eigen

    hdl = zeros(8,1);
    hdl(1) = findobj(Fig.eigen,'Tag','Text3');
    hdl(2) = findobj(Fig.eigen,'Tag','Text4');
    hdl(3) = findobj(Fig.eigen,'Tag','Text5');
    hdl(4) = findobj(Fig.eigen,'Tag','Text6');
    hdl(5) = findobj(Fig.eigen,'Tag','Axes1');
    hdl(6) = findobj(Fig.eigen,'Tag','Listbox1');
    hdl(7) = findobj(Fig.eigen,'Tag','Text1');
    hdl(8) = findobj(Fig.eigen,'Tag','Text2');

    set(hdl(1),'String',num2str(num_pos));
    set(hdl(2),'String',num2str(num_neg));
    set(hdl(3),'String',num2str(num_comp));
    set(hdl(4),'String',num2str(num_zero));
    set(hdl(7),'String',num2str(DAE.n));
    set(hdl(8),'String',num2str(Bus.n));

    autovalori = cell(length(autor),1);
    if num_auto < 10
      d = '';
      e = '';
    elseif num_auto < 100
      d = ' ';
      e = '';
    elseif num_auto < 1000
      d = ' ';
      e = ' ';
    else
      d = ' ';
      e = '  ';
    end
    for i = 1:length(autor)
      if autor(i)>=0
        a = ' ';
      else
        a = '';
      end
      if autoi(i)>=0
        c = '+';
      else
        c = '-';
      end
      if i < 10
        f1 = [d,e];
      elseif i < 100
        f1 = e;
      else
        f1 = '';
      end
      if tipo_plot == 3
        autovalori{i,1} = ['|',char(181),'(A)| #',num2str(i), ...
                           f1, ' ', fvar(abs(auto(i)),9)];
      else
        autovalori{i,1} = [char(181),'(A) #',num2str(i),f1, ...
                           ' ',a,num2str(autor(i)),' ',c, ...
                           ' j',num2str(abs(autoi(i)))];
      end
    end
    set(hdl(6),'String',autovalori,'Value',1)
  end

  Header{3+uno,1} = 'STATISTICS';
  Cols{3+uno} = '';
  Rows{3+uno} = '';
  Matrix{3+uno,1} = [];

  if tipo_mat < 4
    Rows{3+uno}{1,1} = 'NUMBER OF BUSES';
    Matrix{3+uno,1}(1,1) = Bus.n;
  else
    Rows{3+uno}{1,1} = 'DYNAMIC ORDER';
    Matrix{3+uno,1}(1,1) = DAE.n;
  end
  Rows{3+uno}{2,1} = '# OF EIGS WITH Re(mu) < 0';
  Matrix{3+uno,1}(2,1) = num_neg;
  Rows{3+uno}{3,1} = '# OF EIGS WITH Re(mu) > 0';
  Matrix{3+uno,1}(3,1) = num_pos;
  Rows{3+uno}{4,1} = '# OF REAL EIGS';
  Matrix{3+uno,1}(4,1) = num_real;
  Rows{3+uno}{5,1} = '# OF COMPLEX PAIRS';
  Matrix{3+uno,1}(5,1) = num_comp;
  Rows{3+uno}{6,1} = '# OF ZERO EIGS';
  Matrix{3+uno,1}(6,1) = num_zero;

  % save eigenvalues and participation factors in SSSA structure
  SSSA.eigs = auto;
  SSSA.pf = pf;

  if Fig.eigen, axes(hdl(5)); fm_eigen('graph')
  end

  SSSA.report.Matrix = Matrix;
  SSSA.report.Header = Header;
  SSSA.report.Cols = Cols;
  SSSA.report.Rows = Rows;

 case 'report'

  if SSSA.matrix == 4 & DAE.n == 0
    fm_disp('No dynamic component loaded. State matrix is not defined',2)
    return
  end

  if isempty(SSSA.report), fm_eigen('runsssa'), end

  % writing data...
  fm_write(SSSA.report.Matrix,SSSA.report.Header, ...
           SSSA.report.Cols,SSSA.report.Rows)
 
 case 'graph'

  hgca = gca;

  if isempty(SSSA.eigs)
    fm_eigen('runsssa')
  end

  axes(hgca)
  autor = real(SSSA.eigs);
  autoi = imag(SSSA.eigs);
  num_auto = length(SSSA.eigs);
  idx = find(autor > -990);
  if ~isempty(idx)
    autor = autor(idx);
    autoi = autoi(idx);
  end

  if Fig.eigen
    switch SSSA.map
     case 1
      idxn = find(autor < 0);
      idxz = find(autor == 0);
      idxp = find(autor > 0);
      if SSSA.matrix == 4
        hdle = plot(autor(idxn), autoi(idxn),'bx', ...
                    autor(idxz), autoi(idxz),'go', ...
                    autor(idxp), autoi(idxp),'rx');
      else
        hdle = plot(autor(idxn), autoi(idxn),'rx', ...
                    autor(idxz), autoi(idxz),'go', ...
                    autor(idxp), autoi(idxp),'bx');
      end
      hold on
      plot([0,0],ylim,':k');
      plot(xlim,[0,0],':k');
      set(hdle,'MarkerSize',8);
      xlabel('Real');
      ylabel('Imag');
      hold off
      set(hgca,'Tag','Axes1')
     case 2
      surf(real(SSSA.pf))
      set(hgca,'XLim',[1 num_auto],'YLim',[1 num_auto]);
      view(0,90);
      box('on')
      ylabel('Eigenvalues');
      if SSSA.matrix == 4
        xlabel('State Variables');
      else
        xlabel('Buses');
      end
      shading('interp')
      colormap('summer');
      title('Participation Factors')
      set(hgca,'Tag','Axes1');
     case 3
      t = 0:0.01:2*pi+0.01;
      x = cos(t);
      y = sin(t);
      plot(x,y,'k:')
      hold on
      idxn = find(autor < 1);
      idxz = find(autor == 1);
      idxp = find(autor > 1);
      hdle = plot(autor(idxn), autoi(idxn),'bx', ...
                  autor(idxz), autoi(idxz),'go', ...
                  autor(idxp), autoi(idxp),'rx');
      set(hdle,'MarkerSize',8);
      xlabel('Real');
      ylabel('Imag');
      xlim(1.1*xlim);
      ylim(1.1*ylim);
      plot([0,0],1.1*ylim,':k');
      plot(1.1*xlim,[0,0],':k');
      hold off
      set(hgca,'Tag','Axes1');
    end
    set(hgca,'Color',Theme.color11)
  end

end