runpsat.m

电力系统的psat

function runpsat(varargin)
% RUNPSAT run PSAT routine for power system analysis
%
% RUNPSAT([FILE,[PATH]],[PERTFILE,[PERTPATH]],ROUTINE)
%
%   FILE:     string containing the PSAT data file (can be a
%             simulink model)
%   PATH:     string containing the absolute path of the data
%             file (default path is "pwd")
%   PERTFILE: string containing the PSAT perturbation file
%             (default is the empty string)
%   PERTPATH: string containing the absolute path of the
%             perturbation file (default is the empty string)
%   ROUTINE:  name of the routine to be launched:
%
%     General options:
%
%       'data'    => set data file
%       'pert'    => set perturbation file
%       'opensys' => open saved system
%       'savsys'  => save currenst system
%       'pfrep'   => write power flow solution
%       'eigrep'  => write eigenvalue report file
%       'pmurep'  => write PMU placement report file
%       'plot'    => plot TD results (Octave only)
%
%     Routines:
%
%       'pf'      => power flow
%       'cpf'     => continuation power flow
%       'snb'     => SNB computation (direct method)
%       'limit'   => LIB computation
%       'n1cont'  => N-1 contingency analysis
%       'opf'     => optimal power flow
%       'cpfatc'  => ATC computation through CPF analysis
%       'sensatc' => ATC computation through sensitivity
%                    analysis
%       'td'      => time domain simulation
%       'sssa'    => small signal stability analysis
%       'pmu'     => PMU placement
%       'gams'    => OPF through PSAT-GAMS interface
%       'uw'      => CPF through PSAT-UWPFLOW interface
%
%Author:    Federico Milano
%Date:      23-Feb-2004
%Version:   1.0.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.

fm_var

% last input is the routine type
routine = varargin{nargin};

% Simulink models are not supported on GNU/Octave
if Settings.octave & strcmp(routine,'data') & ...
      ~isempty(findstr(varargin{1},'.mdl'))
  fm_disp('Simulink models are not supported on GNU/Octave')
  return
end

% check if the data file has been changed
changedata = strcmp(routine,'data');
if nargin > 1
  changedata = changedata | ~strcmp(varargin{1},File.data);
end

if changedata, Settings.init = 0; end

% check inputs
switch nargin
 case 5
  File.data = varargin{1};
  Path.data = checksep(varargin{2});
  File.pert = varargin{3};
  Path.pert = checksep(varargin{4});
 case 4
  File.data = varargin{1};
  Path.data = checksep(varargin{2});
  File.pert = varargin{3};
  Path.pert = [pwd,filesep];
 case 3
  switch routine
   case 'data'
    File.data = varargin{1};
    Path.data = checksep(varargin{2});
   case 'pert'
    File.pert = varargin{1};
    Path.pert = checksep(varargin{2});
   case 'opensys'
    datafile = varargin{1};
    datapath = checksep(varargin{2});
   otherwise
    File.data = varargin{1};
    Path.data = checksep(varargin{2});
    File.pert = '';
    Path.pert = '';
  end
 case 2
  switch routine
   case 'data'
    File.data = varargin{1};
    Path.data = [pwd,filesep];
   case 'pert'
    File.pert = varargin{1};
    Path.pert = [pwd,filesep];
   case 'opensys'
    datafile = varargin{1};
    datapath = [pwd,filesep];
   otherwise
    File.data = varargin{1};
    Path.data = [pwd,filesep];
    File.pert = '';
    Path.pert = '';
  end
 case 1
  % nothing to do...
 otherwise
  error('Invalid number of arguments: check synthax...')
end

% remove extension from data file (only Matlab files)
if length(File.data) >= 2 & strcmp(routine,'data')
  if strcmp(File.data(end-1:end),'.m')
    File.data = File.data(1:end-2);
  end
end

% remove extension from perturbation file (only Matlab files)
if length(File.pert) >= 2 & strcmp(routine,'pert')
  if strcmp(File.pert(end-1:end),'.m')
    File.pert = File.pert(1:end-2);
  end
end

% set local path as data path to prevent undesired change
% of path within user defined functions
Path.local = Path.data;

% check if the data file is a Simulink model
File.data = strrep(File.data,'.mdl','(mdl)');

if ~isempty(findstr(File.data,'(mdl)'))
  filedata = deblank(strrep(File.data,'(mdl)','_mdl'));
  if exist(filedata) ~= 2 | clpsat.refreshsim | strcmp(routine,'data')
    check = fm_sim;
    if ~check, return, end
  end
end

% launch PSAT computations
switch routine
 case 'data' % set data file
  % checking the consistency of the data file
  if Settings.octave
    localpath = pwd;
    cd(Path.data)
    check = exist(File.data);
    cd(localpath)    
  else
    check = exist([Path.data,File.data]);
  end
  if check ~= 2 & check ~= 4
    fm_disp(['Warning: The selected file is not valid or not in the ' ...
             'current folder!'])
  end
 case 'pert' % set perturbation file
  if Settings.octave
    localpath = pwd;
    cd(Path.pert)
    check = exist(File.pert);
    cd(localpath)    
  else
    check = exist([Path.pert,File.pert]);
  end
  % checking the consistency of the pert file
  if check ~= 2
    fm_disp(['Warning: The selected file is not valid or not in the ' ...
             'current folder!'])
  else
    localpath = pwd;
    cd(Path.pert)
    if Settings.hostver >= 6
      Hdl.pert = str2func(File.pert);
    else
      Hdl.pert = File.pert;
    end
    cd(localpath)
  end
 case 'opensys'
  fm_set('opensys',datafile,datapath)
  Settings.init = 0;
 case 'savesys'
  fm_set('savesys')
 case 'log'
  fm_text(1)
 case 'pfrep'
  fm_report
 case 'eigrep'
  fm_eigen('report')
 case 'pf'   % solve power flow

  if isempty(File.data)
    fm_disp('Set a data file before running Power Flow.',2)
    return
  end

  %comps = Syn.n+Mn.n+Svc.n+Sofc.n+Dfig.n+Cswt.n+Pl.n+Fl.n+Thload.n;

  if clpsat.readfile | Settings.init == 0
    fm_inilf
    filedata = [File.data,'  '];
    filedata = strrep(filedata,'@ ','');

    if ~isempty(findstr(filedata,'(mdl)')) & clpsat.refreshsim
      filedata1 = File.data(1:end-5);
      open_sys = find_system('type','block_diagram');
      donotclose = 0;
      for i = 1:length(open_sys)
        if strcmp(open_sys{i},filedata1)
          donotclose = 1;
        break
        end
      end
      if donotclose
        if strcmp(get_param(filedata1,'Dirty'),'on') | ...
              str2num(get_param(filedata1,'ModelVersion')) > Settings.mv
          check = fm_sim;
          if ~check,
            cd(Path.data)
            return,
          end
        end
      end
    end
    cd(Path.data)
    filedata = deblank(strrep(filedata,'(mdl)','_mdl'));
    a = exist(filedata);
    clear(filedata)
    if a == 2,
      lasterr('');
      try,
        fm_disp('Load data from file...')
        eval(filedata);
      catch,
        fm_disp(lasterr),
        fm_disp(['Something wrong with the data file "',filedata,'"']),
        return
      end
    else,
      fm_disp(['File "',filedata,'" not found or not an m-file'],2)
    end
    cd(Path.local)
    Settings.init = 0;
  end

  if Settings.init

    SW.con = SW.store;
    PV.con = PV.store;
    PQ.con = PQ.store;
    Mn.con = Mn.store;

    % count buses of different types
    if (~isempty(SW.con))
      SW.n = length(SW.con(:,1));
      SW.bus = Bus.int(round(SW.con(:,1)));
      for i = 1:SW.n
        if length(find(SW.bus == SW.bus(i))) > 1
          fm_disp(['Error: More than one slack generator ', ...
                   'connected to the same bus.'],2)
          return
        end
      end
      if length(SW.con(1,:)) < 8
        SW.con = [SW.con(:,1),100*ones(SW.n,1), ...
                  ones(SW.n,1),SW.con(:,[2:end])];
      end
      if length(SW.con(1,:)) == 5
        SW.con = [SW.con, 999*ones(SW.n,1), -999*ones(SW.n,1), ...
                  1.1*ones(SW.n,1), 0.9*ones(SW.n,1)];
      end