fm_cpf.m

电力系统的psat

function [lambda_crit, dl_dp] = fm_cpf(fun)
% FM_CPF continuation power flow routines for computing nose curves
%        and determining critical points (saddle-node bifurcations)
%
% [LAMBDAC] = FM_CPF
% [LAMBDAC,DL_DP] = FM_CPF
%
%       LAMBDAC: loading paramter at critical or limit point
%
%                                            d lambda |
%       DL_DP: sensitivity coefficients      -------- |
%                                               d p   |c
%
%see also CPF structure and FM_CPFFIG
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    16-Sep-2003
%Version:   1.0.1
%
%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

%global etaR

if ~autorun('Continuation Power Flow'), return, end

%  Settings
if CPF.ilim
  ilim = CPF.flow;
else
  ilim = 0;
end
type = ~CPF.sbus;
if type & SW.n == 1 & Supply.n == 1
  if Supply.bus ~= SW.bus
    fm_disp(' * * ')
    fm_disp('Only one Supply found. Single slack bus model will be used.')
    type = 0;
  end
end
stop = CPF.type - 1;
vlim = CPF.vlim;
qlim = CPF.qlim;
perp = ~(CPF.method - 1);

if strcmp(fun,'atc') | strcmp(fun,'gams')
  one = 1;
else
  one = 0;
end

%if PV.n+SW.n == 1
%  type = 0;
%end

if CPF.show
  fm_disp
  if type
    fm_disp('Continuation Power Flow - Distribuited Slack Bus')
  else
    fm_disp('Continuation Power Flow - Single Slack Bus')
  end
  fm_disp(['Data file "',Path.data,File.data,'"'])
  fm_disp
  if perp
    fm_disp('Corrector Method:         Perpendicular Intersection')
  else
    fm_disp('Corrector Method:         Local Parametrization')
  end
  if vlim
    fm_disp('Check Voltage Limits:     Yes')
  else
    fm_disp('Check Voltage Limits:     No')
  end
  if qlim
    fm_disp('Check Generator Q Limits: Yes')
  else
    fm_disp('Check Generator Q Limits: No')
  end
  if ilim
    fm_disp('Check Flow Limits:        Yes')
  else
    fm_disp('Check Flow Limits:        No')
  end
  fm_disp
end

% Initializations of vectors and matrices of state and algebraic
% functions and variables
% ----------------------------------------------------------------------
%warning off

% make a copy of PQ, PV and SW structures
PQold = PQ;
PVold = PV;
SWold = SW;

% disable conversion to impedance for PQ loads
PQ.con(:,8) = 0;

nodyn = 0;

ksw = 0;
kpv = 0;
if SW.n, ksw = SW.con(:,11); end
if PV.n, kpv = PV.con(:,10); end
% fix unconsistent participation factors
if sum(kpv) == 0,
  kpv = ones(PV.n,1);
end

if Syn.n
  Pm_syn = Syn.pm;
  pm_idx = Syn.omega;
  syn_idx = [1:Syn.n]';
  if Tg.n
    pm_idx(Tg.syn) = [];
    syn_idx(Tg.syn) = [];
  end
end

% initialization of the state vector and Jacobian matices
if DAE.n > 0
  DAE.f = ones(DAE.n,1);  % state variable time derivatives
  DAE.x = ones(DAE.n,1);  % state variables
  fm_xfirst;
  if Settings.octave
    DAE.Fx = zeros(DAE.n,DAE.n);   % Dx(f)
    DAE.Fy = zeros(DAE.n,2*Bus.n); % Dy(f)
    DAE.Gx = zeros(2*Bus.n,DAE.n); % Dx(g)
    Flambda = zeros(DAE.n,1);
    Fk = zeros(DAE.n,1);
  else
    DAE.Fx = sparse(DAE.n,DAE.n);   % Dx(f)
    DAE.Fy = sparse(DAE.n,2*Bus.n); % Dy(f)
    DAE.Gx = sparse(2*Bus.n,DAE.n); % Dx(g)
    Flambda = sparse(DAE.n,1);
    Fk = sparse(DAE.n,1);
  end
else  % no dynamic components
  nodyn = 1;
  DAE.n = 1;
  DAE.f = 0;
  DAE.x = 0;
  DAE.Fx = 1;
  if Settings.octave
    DAE.Fy = zeros(1,2*Bus.n);
    DAE.Gx = zeros(2*Bus.n,1);
  else
    DAE.Fy = sparse(1,2*Bus.n);
    DAE.Gx = sparse(2*Bus.n,1);
  end
  Flambda = 0;
  Fk = 0;
end

noDem = 0;
noSup = 0;
sp = ' * ';

if ilim
  tps = Line.con(:,11).*exp(jay*Line.con(:,12)*pi/180);
  r = Line.con(:,8);
  rx = Line.con(:,9);
  chrg = Line.con(:,10)/2;
  z = r + jay*rx;
  y = 1./z;
  g = real(y);
  b = imag(y);
  Imax = Line.con(:,12+ilim);
  switch ilim
   case 1, fw = 'I ';
   case 2, fw = 'P ';
   case 3, fw = 'S ';
  end
  Imax(find(Imax <= 0)) = 99999;  % no limit for undefined limit
end

% the slack bus active power is the power injection obtained
% by the power flow
for i = 1:SW.n,
  k = SW.bus(i);
  length(Snapshot.V); % do not remove this line!!! 
                      % It fixes an Octave bug
  if ~isempty(Snapshot(1).Pg)
    pg = Snapshot(1).Pg;
    if pg(k) ~= 0
      SW.con(i,10) = pg(k);
    end
  end
end

% if no Demand.con is found, the load direction
% is assumed to be the one of the PQ load
if Demand.n
  no_dpq = find(Demand.con(:,3) == 0 & Demand.con(:,4) == 0);
  if ~isempty(no_dpq),
    fm_disp
    for i = 1:length(no_dpq)
      fm_disp(['No power direction found in "Demand.con" for Bus ', ...
               Varname.bus{Demand.bus(no_dpq(i))}])
    end
    fm_disp('Continuation load flow routine may have convergence problems.',2)
    fm_disp
  end
else
  noDem = 1;
  Ppos = find(PQ.con(:,4) > 0);
  Demand.con = zeros(length(Ppos),14);
  Demand.n = length(Ppos);
  Demand.bus = PQ.bus(Ppos);
  Demand.con(:,[1 2 3 4]) = PQ.con(Ppos,[1 2 4 5]);
  Demand.con(:,14) = 1;
  PQ.con(Ppos,[4 5]) = 0;
end

% if no Supply.con is found, the load direction
% is assumed to be the one of the PV or the SW generator
if Supply.n
  no_sp = find(Supply.con(:,3) == 0);
  if ~isempty(no_sp),
    fm_disp
    if length(no_sp) == Supply.n
      fm_disp(['No power directions found in "Supply.con" for all buses.'])
      fm_disp('Remove "Supply" components or set power directions.')
      fm_disp('Continuation power flow interrupted',2)
      if CPF.show
        set(Fig.main,'Pointer','arrow');
      end
      return
    else
      for i = 1:length(no_sp)
        fm_disp(['No power direction found in "Supply.con" for Bus ', ...
                 Varname.bus{Supply.bus(no_sp(i))}])
      end
      fm_disp('Continuation power flow routine may have convergence problems.',2)
      fm_disp
    end
  end
else
  noSup = 1;
  if PV.n
    Ppos = find(PV.con(:,4) > 0);
    Supply.con = zeros(length(Ppos),14);
    Supply.n = length(Ppos);
    Supply.bus = PV.bus(Ppos);
    Supply.con(:,[1 2 3]) = PV.con(Ppos,[1 2 4]);
    Supply.con(:,14) = 1;
    PV.con(Ppos,4) = 0;
  end
  if SW.n
    Ppos = find(SW.con(:,10) ~= 0);
    Supply.n = Supply.n+length(Ppos);
    Supply.bus = [Supply.bus; SW.bus(Ppos)];
    Supply.con(end+[1:length(Ppos)],[1 2 3]) = SW.con(Ppos,[1 2 10]);
    Supply.con((end-length(Ppos)+1):end,14) = 1;
    SW.con(Ppos,10) = 0;
  end
end

% Newton-Raphson routine settings
iter_max = Settings.lfmit;
iterazione = 0;
tol =  CPF.tolc;
tolf = CPF.tolf;
tolv = CPF.tolv;
err_max = tol+1;
proceed = 1;
sigma_corr = 1;
if DAE.n == 0, DAE.n = 1; end

if Settings.octave
  Kjac = zeros(1,2*Bus.n+DAE.n+2);
  Ljac = zeros(1,2*Bus.n+DAE.n+2);
  kjac = zeros(1,2*Bus.n+DAE.n+1);
else
  Kjac = sparse(1,2*Bus.n+DAE.n+2);
  Ljac = sparse(1,2*Bus.n+DAE.n+2);
  kjac = sparse(1,2*Bus.n+DAE.n+1);
end

if Supply.n
  Kjac(1,DAE.n+SW.bus) = 1;
else
  Kjac(1,end-1) = 1;
end

if Supply.n
  kjac(DAE.n+SW.bus) = 1;
else
  kjac(1,end) = 1;
end

% A PQ bus is chosen for visualizing the CPF progress in the main window
PQidx = 0;
for i = 1:PQ.n
  a = find(PQ.bus(i) == PV.bus);
  b = find(PQ.bus(i) == SW.bus);
  if isempty(a) & isempty(b)
    PQidx = PQ.bus(i);
    break
  end
end
if ~PQidx % absence of PQ buses
  for i = 1:Demand.n
    a = find(Demand.bus(i) == PV.bus);
    b = find(Demand.bus(i) == SW.bus);
    if isempty(a) & isempty(b)
      PQidx = Demand.bus(i);
      break
    end
  end
  if ~PQidx
    if ~perp
      perp = 1; 
      fm_disp('* * * Switch to perpendicular intersection * * * ')      
    end
    PQidx = 1;
  end
end

% chose bus to be monitored during CPF analysis
PQidxplot = PQidx;

% distributed slack bus coefficients
ksu = zeros(Supply.n,1);
for i = 1:PV.n
  idx = find(Supply.bus == PV.bus(i));
  if idx
    ksu(idx) = PV.con(i,10);
  end
end
for i = 1:SW.n
  idx = find(Supply.bus == SW.bus(i));
  if idx
    ksu(idx) = SW.con(i,11);
  end
end

% *********************************************************************
% Continuation Power Flow Routine
% *********************************************************************

tic

fm_status('cpf','init',12,{'m'},{'-'},{Theme.color11},[0 1.2])

l_vect = [];

lambda = CPF.linit;
lambda_crit = CPF.linit;
kg = 0;
lambda_old = -1;
Jsign = 0;
Jsign_old = 0;
Sflag = 1;
ksign = 1;
Qmax_idx = [];
Qmin_idx = [];
Vmax_idx = [];
Vmin_idx = [];
Iij_idx = [];
Iji_idx = [];
Iij = ones(Line.n,1);
Iji = ones(Line.n,1);
Qswmax_idx = [];
Qswmin_idx = [];
PQVmax = ones(PQ.n,1);
PQVmin = ones(PQ.n,1);
fm_out(0,0,0)
fm_out(1,0,0)

while err_max > tol

  cpfmsg = [];

  if iterazione >= CPF.nump & ~strcmp(fun,'gams')
    fm_disp([sp,'Reached maximum number of points.'])
    break
  end
  if Fig.main
    if ~get(Fig.main,'UserData'), break, end
  end

  if isempty(Line.Y)