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://www.uclm.es/area/gsee/Web/Federico
%
% Copyright (C) 2002-2008 Federico Milano

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

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

  uno = 0;
  tipo_mat = SSSA.matrix;
  tipo_plot = SSSA.map;
  SSSA.report = [];

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

  % build eigenvalue 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-2008';
  Header{1,1}{6,1} = 'e-mail:  [e-mail]';
  Header{1,1}{7,1} = 'website: http://www.uclm.es/area/gsee/Web/Federico';
  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

    As = DAE.Fx - DAE.Fy*(DAE.Gy\DAE.Gx) - 1e-6*speye(DAE.n);
    if tipo_plot == 3
      As = (As+8*speye(DAE.n))/(As-8*speye(DAE.n));
    end

    [auto,autor,autoi,num_auto,pf] = compute_eigs(As);

    names = cellstr(fm_strjoin('Eig As #',num2str([1:num_auto]')));

    Header{2,1} = 'STATE MATRIX EIGENVALUES';
    Cols{2,1} = {'Eigevalue', 'Most Associated States', ...
                 'Real part','Imag. Part','Pseudo-Frequency','Frequency'};
    Matrix{2,1} = zeros(num_auto,4);
    Matrix{2,1}(:,[1 2]) = [autor, autoi];

    for i = 1:num_auto;
      if autoi(i) == 0
        [part, idxs] = max(pf(i,:));
        stat = Varname.uvars{idxs};
        pfrec = 0;
        frec = 0;
      else
        [part, idxs] = sort(pf(i,:));
        stat = [Varname.uvars{idxs(end)},', ',Varname.uvars{idxs(end-1)}];
        pfrec = abs(imag(auto(i))/2/3.1416);
        frec = abs(auto(i))/2/3.1416;
      end
      Rows{2,1}{i,1} = names{i};
      Rows{2,1}{i,2} = stat;
      Matrix{2,1}(i,3) = pfrec;
      Matrix{2,1}(i,4) = frec;
    end

    [uno,Header,Cols,Rows,Matrix] = ...
        write_pf(pf,DAE.n,Varname.uvars,names,Header,Cols,Rows,Matrix);

  elseif tipo_mat == 3

    Jlf = build_gy(Line);
    if DAE.m > 2*Bus.n
      x3 = 1:2*Bus.n;
      Jlf = Jlf(x3,x3);
    end
    x1 = Bus.a;
    x2 = Bus.v;
    vbus = [getbus(SW,'a');getbus(SW,'v');getbus(PV,'v')];
    Jlf(vbus,:) = 0;
    Jlf(:,vbus) = 0;
    Jlf = Jlf + sparse(vbus,vbus,999,2*Bus.n,2*Bus.n);
    Jlfptheta = Jlf(x1,x1)+1e-5*speye(Bus.n,Bus.n);
    elementinulli = find(diag(Jlfptheta == 0));
    if ~isempty(elementinulli)
      for i = 1:length(elementinulli)
        Jlfptheta(elementinulli(i),elementinulli(i)) = 1;
      end
    end
    Jlfr = Jlf(x2,x2) - Jlf(x2,x1)*(Jlfptheta\Jlf(x1,x2));

    [auto,autor,autoi,num_auto,pf] = compute_eigs(Jlfr);

    names = cellstr(fm_strjoin('Eig Jlfr #',num2str([1:num_auto]')));

    Header{2,1} = 'EIGENVALUES OF THE STANDARD POWER JACOBIAN MATRIX';
    Cols{2,1} = {'Eigevalue', 'Most Associated Bus', 'Real part', ...
                 'Imaginary Part'};
    Rows{2,1} = names;
    Matrix{2,1} = [autor, autoi];

    for i = 1:num_auto;
      if autoi(i) == 0
        [part, idxs] = max(pf(i,:));
        stat = Bus.names{idxs};
      else
        [part, idxs] = sort(pf(i,:));
        stat = [Bus.names{idxs(end)},', ',Bus.names{idxs(end-1)}];
      end
      Rows{2,1}{i,1} = names{i};
      Rows{2,1}{i,2} = stat;
    end

    [uno,Header,Cols,Rows,Matrix] = ...
        write_pf(pf,Bus.n,Bus.names,names,Header,Cols,Rows,Matrix);

  elseif tipo_mat == 2

    x1 = Bus.a;
    x2 = Bus.v;
    if DAE.m > 2*Bus.n
      x3 = 1:2*Bus.n;
      x4 = 2*Bus.n+1:DAE.m;
      Gy = DAE.Gy(x3,x3)-DAE.Gy(x3,x4)*(DAE.Gy(x4,x4)\DAE.Gy(x4,x3));
    else
      Gy = DAE.Gy;
    end
    Jlfvr = Gy(x2,x2)-Gy(x2,x1)*(Gy(x1,x1)\Gy(x1,x2));
    vbus = [getbus(SW);getbus(PV)];
    Jlfvr = Jlfvr + sparse(vbus,vbus,998,Bus.n,Bus.n);

    [auto,autor,autoi,num_auto,pf] = compute_eigs(Jlfvr);

    names = cellstr(fm_strjoin('Eig Jlfv #',num2str([1:num_auto]')));

    Header{2,1} = 'EIGENVALUES OF THE COMPLETE POWER JACOBIAN MATRIX';
    Cols{2,1} = {'Eigevalue', 'Most Associated Bus', 'Real part', ...
                   'Imaginary Part'};
    Rows{2,1} = names;
    Matrix{2,1} = [autor, autoi];

    for i = 1:num_auto;
      if autoi(i) == 0
        [part, idxs] = max(pf(i,:));
        stat = Bus.names{idxs};
      else
        [part, idxs] = sort(pf(i,:));
        stat = [Bus.names{idxs(end)},', ',Bus.names{idxs(end-1)}];
      end
      Rows{2,1}{i,1} = names{i};
      Rows{2,1}{i,2} = stat;
    end

    [uno,Header,Cols,Rows,Matrix] = ...
        write_pf(pf,Bus.n,Bus.names,names,Header,Cols,Rows,Matrix);

  elseif tipo_mat ==  1