fm_eigen.m

电力系统的psat

function fm_eigen(type)
% FM_EIGEN compute eigenvalues of static and dynamic Jacobian
%          matrices
%
% FM_EIGEN(TYPE,REPORT)
%       TYPE     1 -> Jlfd eigenvalues
%                2 -> Jlfv eigenvalues
%                3 -> Jlf eigenvalues
%                4 -> As eigenvalues
%       REPORT   1 -> create report file
%                0 -> no report file
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    05-Mar-2004
%Version:   1.1.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2005 Federico Milano
%
% This toolbox is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2.0 of the License, or
% (at your option) any later version.
%
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

global DAE Bus Settings Varname File Path
global Snapshot PV SW Fig Theme SSSA

switch type

 case 'matrix' % set matrix type

  mat = zeros(4,1);
  mat(1) = findobj(gcf,'Tag','Checkbox1');
  mat(2) = findobj(gcf,'Tag','Checkbox2');
  mat(3) = findobj(gcf,'Tag','Checkbox3');
  mat(4) = findobj(gcf,'Tag','Checkbox4');

  tplot = zeros(3,1);
  tplot(1) = findobj(gcf,'Tag','Radiobutton1');
  tplot(2) = findobj(gcf,'Tag','Radiobutton2');
  tplot(3) = findobj(gcf,'Tag','Radiobutton3');

  ca = find(mat == gcbo);
  vals = zeros(4,1);
  vals(ca) = 1;
  for i = 1:4
    set(mat(i),'Value',vals(i))
  end
  a = get(tplot(3),'Value');
  if ca == 4
    set(tplot(3),'Enable','on')
  else
    if a
      set(tplot(3),'Value',0)
      set(tplot(1),'Value',1)
      SSSA.map = 1;
    end
    set(tplot(3),'Enable','off')
  end
  if ca == 4 & SSSA.neig > DAE.n-1,
    set(findobj(Fig.eigen,'Tag','EditText1'),'String','1')
    SSSA.neig = 1;
  elseif ca < 4 & SSSA.neig > Bus.n-1,
    set(findobj(Fig.eigen,'Tag','EditText1'),'String','1')
    SSSA.neig = 1;
  end
  SSSA.matrix = ca;
  SSSA.report = [];

 case 'map' % set eigenvalue map type

  tplot = zeros(3,1);
  tplot(1) = findobj(gcf,'Tag','Radiobutton1');
  tplot(2) = findobj(gcf,'Tag','Radiobutton2');
  tplot(3) = findobj(gcf,'Tag','Radiobutton3');

  ca = find(tplot == gcbo);
  vals = zeros(3,1);
  vals(ca) = 1;
  for i = 1:3
    set(tplot(i),'Value',vals(i))
  end
  SSSA.map = find(vals);
  SSSA.report = [];

 case 'neig' % Set number of eigs to be computed

  if SSSA.matrix == 4
    amax = DAE.n;
  else
    amax = Bus.n;
  end

  number = get(gcbo,'String');
  try
    a = round(str2num(number));
    if a > 0 & a < amax
      SSSA.neig = a;
      SSSA.report = [];
    else
      set(gcbo,'String',num2str(SSSA.neig));
    end
  catch
    set(gcbo,'String',num2str(SSSA.neig));
  end

 case 'method' % Set method for eigenvalue computation

  t1 = findobj(Fig.eigen,'Tag','Radiobutton1');
  t3 = findobj(Fig.eigen,'Tag','Radiobutton2');

  a = get(gcbo,'Value');
  hedit = findobj(Fig.eigen,'Tag','EditText1');
  if a == 1
    set(hedit,'Enable','off')
    set(t3,'Enable','on')
  else
    set(hedit,'Enable','on')
    if get(t3,'Value')
      set(t3,'Value',0)
      set(t1,'Value',1)
      SSSA.map = 1;
    end
    set(t3,'Enable','off')
  end
  SSSA.method = a;
  SSSA.report = [];

 case 'runsssa'

  % check for data file
  if isempty(File.data)
    fm_disp('Set a data file before running eigenvalue analysis.',2)
    return
  end

  % check for initial power flow solution
  if ~Settings.init
    fm_disp('Solve base case power flow...')
    Settings.show = 0;
    fm_set('lf')
    Settings.show = 1;
    if ~Settings.init, return, end
  end

  meth = {'LM';'SM';'LR';'SR';'LI';'SI'};
  neig = SSSA.neig;
  if Settings.octave
    opts = 0;
    if SSSA.method-1
      fm_disp(['Only full eigenvalue analysis is supported on ' ...
               'Octave'])
    end
  else
    opts = SSSA.method-1;
  end
  uno = 0;
  tipo_mat = SSSA.matrix;
  tipo_plot = SSSA.map;
  opt.disp = 0;
  SSSA.report = [];
  pf = [];

  if isempty(Bus.n)
    fm_disp('No loaded system. Eigenvalue computation cannot be run.',2)
    return
  end

  % build eigenvalu names
  if (Settings.vs == 0), fm_idx(3), end

  % initialize report structures
  Header{1,1}{1,1} = 'EIGENVALUE REPORT';
  Header{1,1}{2,1} = ' ';
  Header{1,1}{3,1} = ['P S A T  ',Settings.version];
  Header{1,1}{4,1} = ' ';
  Header{1,1}{5,1} = 'Author:  Federico Milano, (c) 2002-2005';
  Header{1,1}{6,1} = 'e-mail:  [e-mail]';
  Header{1,1}{7,1} = 'website: http://thunderbox.uwaterloo.ca/~fmilano';
  Header{1,1}{8,1} = ' ';
  Header{1,1}{9,1} = ['File:  ', Path.data,strrep(File.data,'(mdl)','.mdl')];
  Header{1,1}{10,1} = ['Date:  ',datestr(now,0)];

  Matrix{1,1} = [];
  Cols{1,1} = '';
  Rows{1,1} = '';

  if tipo_mat == 4

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

    Fy = DAE.Fy;
    Gx = DAE.Gx;
    a = [PV.bus; SW.bus];
    b2 = 1:Bus.n;
    c = [];
    for i = 1:length(a)
      c = [c; find(b2 == a(i))];
    end
    if ~isempty(c)
      b2(c) = [];
    end
    a = [SW.bus];
    b1 = 1:Bus.n;
    c = [];
    for i = 1:length(a)
      c = [c; find(b1 == a(i))];
    end
    if ~isempty(c)
      b1(c) = [];
    end
    b = [b1, b2+Bus.n];
    if Settings.octave
      As = DAE.Fx - Fy(:,b)*([DAE.Jlfv(b,b)-1e-6*eye(length(b))]\ ...
                             Gx(b,:));
    else
      As = DAE.Fx - Fy(:,b)*([DAE.Jlfv(b,b)-1e-6*speye(length(b))]\ ...
                             Gx(b,:));
      %As = DAE.Fx - Fy(:,b)*(DAE.Jlfv(b,b)\Gx(b,:));
    end
    if tipo_plot == 3
      if Settings.octave
        As = (As+8*speye(DAE.n))/(As-8*eye(DAE.n));
      else
        As = (As+8*speye(DAE.n))/(As-8*speye(DAE.n));
      end
    end
    if opts
      auto = eigs(As,neig,meth{opts},opt);
    else
      [V, auto] = eig(full(As));
      auto = diag(auto);
    end

    auto = round(auto/Settings.lftol)*Settings.lftol;
    num_auto = length(auto);
    autor = real(auto);
    autoi = imag(auto);

    for i = 1:num_auto;
      names{i,1} = strrep(Varname.uautostate{i},'_As',' As ');
      names{i,1} = strrep(names{i,1},'eigenvalue','Eig');
    end
    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} = 'STATE MATRIX EIGENVALUES';
    Cols{2,1} = {'Eigevalue', 'Most Associated States', ...
                 'Real part','Imag. Part','Frequency'};
    Matrix{2,1} = zeros(num_auto,3);
    Matrix{2,1}(:,[1 2]) = [autor, autoi];

    for i = 1:num_auto;
      if imag(auto(i)) == 0 & ~opts
        [part, idxs] = max(pf(i,:));
        stat = Varname.ux{idxs};
        frec = 0;
      elseif ~opts
        [part, idxs] = sort(pf(i,:));
        stat = [Varname.ux{idxs(end)},', ',Varname.ux{idxs(end-1)}];
        frec = abs(imag(auto(i))/2/3.1416);
      else
        stat = '---';
        frec = 0;
      end
      Rows{2,1}{i,1} = names{i};
      Rows{2,1}{i,2} = stat;
      Matrix{2,1}(i,3) = frec;
    end
    if ~opts
      uno = fix(DAE.n/5);
      due = rem(DAE.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.ux{5*(k-1)+1:min(5*k,DAE.n)}};
        Rows{2+k} = names;
        Matrix{2+k,1} = pf(:,5*(k-1)+1:min(5*k,DAE.n));
      end
    end

  elseif tipo_mat == 3

    Jlf = DAE.Jlf;
    Jlf(SW.bus,:) = 0;
    Jlf(SW.bus+Bus.n,:) = 0;
    Jlf(:,SW.bus) = 0;
    Jlf(:,SW.bus+Bus.n) = 0;
    if Settings.octave
      Jlf(SW.bus,SW.bus) = Jlf(SW.bus,SW.bus) + eye(SW.n);
      Jlf(SW.bus+Bus.n,SW.bus+Bus.n) = Jlf(SW.bus+Bus.n,SW.bus+Bus.n) ...
          + 999*eye(SW.n);
    else
      Jlf(SW.bus,SW.bus) = Jlf(SW.bus,SW.bus) + speye(SW.n);
      Jlf(SW.bus+Bus.n,SW.bus+Bus.n) = Jlf(SW.bus+Bus.n,SW.bus+Bus.n) ...
          + 999*speye(SW.n);
    end
    Jlf(PV.bus+Bus.n,:) = 0;
    Jlf(:,PV.bus+Bus.n) = 0;
    if Settings.octave
      Jlf(PV.bus+Bus.n,PV.bus+Bus.n) = Jlf(PV.bus+Bus.n,PV.bus+Bus.n) ...
          + 999*eye(PV.n);
    else
      Jlf(PV.bus+Bus.n,PV.bus+Bus.n) = Jlf(PV.bus+Bus.n,PV.bus+Bus.n) ...
          + 999*speye(PV.n);
    end
    Jlfptheta = Jlf(1:Bus.n,1:Bus.n)+diag(Bus.n)*1e-5;
    elementinulli = find(diag(Jlfptheta == 0));
    if ~isempty(elementinulli)
      for i = 1:length(elementinulli)
        Jlfptheta(elementinulli(i),elementinulli(i)) = 1;
      end
    end
    Jlfr = (Jlf(Bus.n+1:2*Bus.n,Bus.n+1:2*Bus.n) - ...
            Jlf(Bus.n+1:2*Bus.n, 1:Bus.n)*...
      (Jlfptheta\Jlf(1:Bus.n,Bus.n+1:2*Bus.n)));
    if opts
      auto = eigs(Jlfr,neig,meth{opts},opt);
    else
      [V, auto] = eig(full(Jlfr));
      auto = diag(auto);
    end
    auto = round(auto/Settings.lftol)*Settings.lftol;
    num_auto = length(auto);
    autor = real(auto);
    autoi = imag(auto);

    for i = 1:num_auto
      names{i,1} = strrep(Varname.uautojlfr{i},'_Jlfr',' Jlf');
      names{i,1} = strrep(names{i,1},'eigenvalue','Eig');
    end
    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 STANDARD POWER JACOBIAN MATRIX';