fm_eigen.m

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

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-2006 Federico Milano

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;
  opts = SSSA.method-1;
  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(2), 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-2006';
  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

    b1 = setdiff([1:Bus.n],SW.bus');
    b2 = setdiff([1:Bus.n],[PV.bus; SW.bus]');
    b = [b1, b2+Bus.n];
    As = DAE.Fx - DAE.Fy(:,b)*([DAE.Jlfv(b,b)-1e-6*speye(length(b))]\DAE.Gx(b,:));
    As = DAE.Fx - DAE.Fy*(DAE.Jlfv\DAE.Gx);
    if tipo_plot == 3
      As = (As+8*speye(DAE.n))/(As-8*speye(DAE.n));
    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);
    names = cellstr(strcat('Eig As #',num2str([1:num_auto]')));

    if ~opts
      W = pinv(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.uvars{idxs};
        frec = 0;
      elseif ~opts
        [part, idxs] = sort(pf(i,:));
        stat = [Varname.uvars{idxs(end)},', ',Varname.uvars{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.uvars{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;
    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);
    Jlf(PV.bus+Bus.n,:) = 0;
    Jlf(:,PV.bus+Bus.n) = 0;
    Jlf(PV.bus+Bus.n,PV.bus+Bus.n) = Jlf(PV.bus+Bus.n,PV.bus+Bus.n) ...
        + 999*eye(PV.n);
    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);
    names = cellstr(strcat('Eig Jlfr #',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 STANDARD 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 == 2