runpsat.m

电力系统的psat

      if length(SW.con(1,:)) == 10
        SW.con = [SW.con, ones(SW.n,1)];
      end
      DAE.V(SW.bus) = SW.con(:,4);
      DAE.a(SW.bus) = SW.con(:,5);
    else
      fm_disp('Error: No slack bus found.',2)
      return
    end

    if (~isempty(PV.con))
      PV.n = length(PV.con(:,1));
      PV.bus = Bus.int(round(PV.con(:,1)));
      for i = 1:PV.n
        if length(find(PV.bus == PV.bus(i))) > 1
          fm_disp(['Error: More than one PV generator ', ...
                   'connected to the same bus.'],2)
          return
        end
      end
      DAE.V(PV.bus) = PV.con(:,5);
      if length(PV.con(1,:)) == 5
        PV.con = [PV.con, 999*ones(PV.n,1), ...
                  -999*ones(PV.n,1), 1.1*ones(PV.n,1), ...
                  0.9*ones(PV.n,1), ones(PV.n,1)];
      end
      if length(PV.con(1,:)) == 9; PV.con = [PV.con, ones(PV.n,1)]; end
    end

    % check of distributed slack bus model consistency
    if ~isempty(PV.con)
      idx = find(~[PV.con(:,10); SW.con(:,11)]);
    else
      idx = find(~SW.con(:,11));
    end
    if length(idx) == PV.n+SW.n, SW.con(:,11) = 1; end

    if (~isempty(PQ.con))
      PQ.n = length(PQ.con(:,1));
      PQ.bus = Bus.int(round(PQ.con(:,1)));
      for i = 1:PQ.n
        if length(find(PQ.bus == PQ.bus(i))) > 1
          fm_disp('Error: More than one PQ load connected to the same bus.',2)
          return
        end
      end
      if length(PQ.con(1,:)) == 5
        PQ.con = [PQ.con, 1.2*ones(PQ.n,1), ...
                  0.8*ones(PQ.n,1), zeros(PQ.n,1)];
      elseif length(PQ.con(1,:)) == 7
        PQ.con = [PQ.con, zeros(PQ.n,1)];
      end
      for i = 1:PQ.n
        if PQ.con(i,6) == 0, PQ.con(i,6) = 1.2; end
        if PQ.con(i,7) == 0, PQ.con(i,7) = 0.8; end
      end
      PQ.P0 = PQ.con(:,4);
      PQ.Q0 = PQ.con(:,5);
    end

    if ~isempty(Mn.con)
      Mn.n = length(Mn.con(:,1));
      Mn.bus = Bus.int(round(Mn.con(:,1)));
      Mn.init = find(~Mn.con(:,8));
    end

    % fix component parameters according to current bases
    if Settings.conv
      if SW.n
        SW.con(:,6)  = SW.con(:,6).*SW.con(:,2)/Settings.mva;
        SW.con(:,7)  = SW.con(:,7).*SW.con(:,2)/Settings.mva;
        SW.con(:,10) = SW.con(:,10).*SW.con(:,2)/Settings.mva;
      end
      if PQ.n
        PQ.con(:,4) = PQ.con(:,4).*PQ.con(:,2)/Settings.mva;
        PQ.con(:,5) = PQ.con(:,5).*PQ.con(:,2)/Settings.mva;
        PQ.P0 = PQ.con(:,4);
        PQ.Q0 = PQ.con(:,5);
      end
      if PV.n
        PV.con(:,4) = PV.con(:,4).*PV.con(:,2)/Settings.mva;
        PV.con(:,6) = PV.con(:,6).*PV.con(:,2)/Settings.mva;
        PV.con(:,7) = PV.con(:,7).*PV.con(:,2)/Settings.mva;
      end
      if Mn.n
        Mn.con(:,4) = Mn.con(:,4).*Mn.con(:,2)/Settings.mva;
        Mn.con(:,5) = Mn.con(:,5).*Mn.con(:,2)/Settings.mva;
      end
    end
    DAE.n = DAE.npf;

  end

  filedata = deblank(strrep(File.data,'(mdl)','_mdl'));
  if Settings.static % do not use dynamic components
    for i = 1:Comp.n
      comp_num = Comp.number{i};
      comp_name = strrep(comp_num,'.n','.con');
      comp_con = eval(['~isempty(',comp_name,')']);
      if comp_con & ~Comp.prop(i,6)
        eval([comp_name,' = [];']);
      end
    end
  end

  fm_spf
  SNB.init = 0;
  LIB.init = 0;
  CPF.init = 0;
  OPF.init = 0;

 case 'opf'  % solve optimal power flow
  fm_set('opf')
 case 'cpf'  % solve continuation power flow
  fm_cpf('main');
 case 'cpfatc'  % find ATC of the current system
  opftype = OPF.type;
  OPF.type = 4;
  fm_atc
  OPF.type = opftype;
 case 'sensatc'
  opftype = OPF.type;
  OPF.type = 5;
  fm_atc
  OPF.type = opftype;
 case 'n1cont'
  fm_n1cont;
 case 'td'   % solve time domain simulation
  fm_int
 case 'sssa' % solve small signal stability analyisis
  fm_eigen('runsssa')
 case 'snb'
  fm_snb
 case 'lib'
  fm_limit
 case 'pmu'
  fm_pmuloc;
 case 'pmurep'
  fm_pmurep;
 case 'gams' % solve OPF using the PSAT-GAMS interface
  fm_gams
 case 'uw'   % solve CPF using the PSAT-UWPFLOW interface
  fm_uwpflow('init')
  fm_uwpflow('uwrun')
 case 'plot'
  if ~Settings.octave
    fm_disp('This option is supported only on GNU/Octave')
    return
  end
  if isempty(Varout.t)
    fm_disp('No data is available for plotting')
    return
  end
  value = menu('Plot variables:','States','Voltages',['Active ' ...
                      'Powers'],'Reactive Powers',['Generator ' ...
                      'speeds'],'Generator angles');
  switch value
   case 1
    if ~DAE.n
      fm_disp('No dynamic component is loaded')
      return
    end
    octplot(DAE.n,Varout.t,Varout.x,Varname.fx)
   case 2
    if ~Bus.n
      fm_disp('No bus is present in the current network')
      return
    end
    octplot(DAE.n,Varout.t,Varout.V,Varname.fV)
   case 3
    if ~Bus.n
      fm_disp('No bus is present in the current network')
      return
    end
    octplot(DAE.n,Varout.t,Varout.p,Varname.fP)
   case 4
    if ~Bus.n
      fm_disp('No bus is present in the current network')
      return
    end
    octplot(DAE.n,Varout.t,Varout.q,Varname.fQ)
   case 5
    if ~Syn.n
      fm_disp('No synchronous generator is loaded')
      return
    end
    octplot(Syn.n,Varout.t,Varout.x(:,Syn.omega), ...
            Varname.fx(Syn.omega))
   case 6
    if ~Syn.n
      fm_disp('No synchronous generator is loaded')
      return
    end
    octplot(Syn.n,Varout.t,Varout.x(:,Syn.delta), ...
            Varname.fx(Syn.delta))
  end
 otherwise   % give an error message and exit
  error(['"',routine,'" is an invalid routine identifier.'])
end

% -------------------------------------------------------------
% restore SW, PV, PQ and Mn structures
function [SW,PV,PQ,Mn] = restorepf(datafile,datapath)

global Settings

localpath = pwd;

% Slack variables
SW.con = [];
SW.n = 0;
SW.bus = [];

% PV variables
PV.con = [];
PV.n = 0;
PV.bus = [];
PV.pq.con = [];
PV.pq.n = 0;
PV.pq.bus = [];

% PQ  variables
PQ.con = [];
PQ.n = 0;
PQ.bus = [];
PQ.P0 = [];
PQ.Q0 = [];

% Monomial Load  variables
Mn.con = [];
Mn.n = 0;
Mn.bus = [];
Mn.init = [];

filedata = [datafile,'  '];
filedata = strrep(filedata,'@ ','');

if ~isempty(findstr(filedata,'(mdl)'))
  filedata1 = datafile(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(datapath),
        return,
      end
    end
  end
end
cd(datapath)
filedata = deblank(strrep(filedata,'(mdl)','_mdl'));
a = exist(filedata);
clear(filedata)

if a == 2,
  lasterr('');
  try,
    fm_disp('Reload static 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(localpath)

% ----------------------------------------------------------------
function octplot(n,x,y,s,xl)

global Settings

plot(x,y(:,1),['1;',s{1},';'])
hold on
for i = 2:n
  FMT = [num2str(rem(i-1,6)+1),';',s{i},';'];
  plot(x,y(:,i),FMT)
end
xlabel(Settings.xlabel)
hold off

% ----------------------------------------------------------------
function string = checksep(string)

if ~strcmp(string(end),filesep)
  string = [string,filesep];
end