fm_cpf.m

电力系统分析计算程序

function lambda_crit = fm_cpf(fun)
% FM_CPF continuation power flow routines for computing nose curves
%        and determining critical points (saddle-node bifurcations)
%
% [LAMBDAC] = FM_CPF
%
%       LAMBDAC: loading paramter at critical or limit point
%
%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://www.uclm.es/area/gsee/Web/Federico
%
% Copyright (C) 2002-2008 Federico Milano

fm_var

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

lambda_crit = [];
lastwarn('')

%  Settings
if CPF.ilim
  ilim = CPF.flow;
else
  ilim = 0;
end
type = double(~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);
hopf = 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
% ----------------------------------------------------------------------

% disable conversion to impedance for PQ loads
forcepq = Settings.forcepq;
Settings.forcepq = 1;

nodyn = 0;

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

PV2PQ = Settings.pv2pq;
noDem = 0;
noSup = 0;
sp = ' * ';

Settings.pv2pq = 0;

Imax = getflowmax(Line,ilim);
switch ilim
  case 1, fw = 'I ';
  case 2, fw = 'P ';
  case 3, fw = 'S ';
end

if CPF.areaannualgrowth, growth(Areas,'area'), end
if CPF.regionannualgrowth, growth(Regions,'region'), 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 = findzero(Demand);
  if ~isempty(no_dpq),
    fm_disp
    for i = 1:length(no_dpq)
      fm_disp(['No power direction found in "Demand.con" for Bus ', ...
        Bus.names{Demand.bus(no_dpq(i))}])
    end
    fm_disp('Continuation load flow routine may have convergence problems.',2)
    fm_disp
  end
else
  noDem = 1;
  idx = findpos(PQ);
  Demand = add(Demand,dmdata(PQ,idx));
  PQ = pqzero(PQ,idx);
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 & ~Syn.n
  no_sp = findzero(Supply);
  if ~isempty(no_sp),
    fm_disp
    if length(no_sp) == Supply.n
      fm_disp(['No power directions found in "Supply.con" for all buses.'])
      if noDem
        fm_disp('Supply data will be ignored.')
        noSup = 1;
        Supply = remove(Supply,[1:Supply.n]);
        if qlim & type, SW = move2sup(SW); end
        %SW = move2sup(SW);
        %PV = move2sup(PV);
      else
        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
      end
    else
      for i = 1:length(no_sp)
        fm_disp(['No power direction found in "Supply.con" for Bus ', ...
          Bus.names{Supply.bus(no_sp(i))}])
      end
      fm_disp('Continuation power flow routine may have convergence problems.',2)
      fm_disp
    end
  end
else
  noSup = 1;
  Supply = remove(Supply,[1:Supply.n]);
  if qlim & type, SW = move2sup(SW); 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

Kjac = sparse(1,DAE.m+DAE.n+2);
Ljac = sparse(1,DAE.m+DAE.n+2);
kjac = sparse(1,DAE.m+DAE.n+1);
Kjac(1,DAE.n+SW.refbus) = 1;
kjac(1,DAE.n+SW.refbus) = 1;

% chose a PQ to display the CPF progress in the main window
PQidx = pqdisplay(PQ);
fm_snap('cleansnap')

if ~PQidx % absence of PQ buses
  PQidx = pqdisplay(Demand);
  if ~PQidx
    if ~perp
      perp = 1;
      fm_disp('* * * Switch to perpendicular intersection * * * ')
    end
    PQidx = 1;
  end
end

PQvdx = PQidx + Bus.n;

% ---------------------------------------------------------------------
% Continuation Power Flow Routine
% ---------------------------------------------------------------------

tic

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

l_vect = [];

lambda = CPF.linit;
DAE.lambda = lambda;
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 = [];

pqlim(PQ,vlim,0,0,0);
fm_out(0,0,0)
if nodyn
  DAE.Fx = 1;
  Varout.idx = Varout.idx+1;
end
fm_out(1,0,0)

% initial Jacobian Gk
DAE.Gk = sparse(DAE.m,1);
if type, Gkcall(Supply), end
if type, Gkcall(PV), end
Gkcall(SW)

DAE.kg = 0;
Gycall(Line)
if noSup
  Gyreactive(PV)
else
  Gycall(PV)
end
Gyreactive(SW)

% First predictor step
% ---------------------------------------------------------------

Jsign = 1;
count_qmin = 0;
kqmin = 1;
d_lambda = 0;
lambda_old = CPF.linit;
inc = zeros(DAE.n+DAE.m+2,1);
Ljac(end) = 1;

% critical point
y_crit = DAE.y;
x_crit = DAE.x;
k_crit = kg;
l_crit = lambda;

while err_max > tol

  cpfmsg = [];

  % corrector step
  % ---------------------------------------------------------------

  lim_hb = 0;
  lim_lib = 0;
  y_old = DAE.y;
  x_old = DAE.x;
  kg_old = kg;
  lambda_old = lambda;
  corrector = 1;
  if PV.n, inc(DAE.n+getbus(PV,'v')) = 0; end
  if SW.n, inc(DAE.n+getbus(SW,'v')) = 0; end

  while corrector

    modinc = norm(inc);
    if Fig.main
      if ~get(Fig.main,'UserData'), break, end
    end

    iter_corr = 0;
    err_corr = 1;

    while err_corr > tol

      if (iter_corr > iter_max), break, end
      if Fig.main
        if ~get(Fig.main,'UserData'), break, end
      end

      DAE.lambda = lambda;
      DAE.Fl = sparse(DAE.n,1);
      DAE.Gl = sparse(DAE.m,1);
      DAE.Gk = sparse(DAE.m,1);

      % call component functions

      fm_call('kg')

      if nodyn, DAE.Fx = 1; end

      glambda(Demand,lambda)
      glambda(Supply,lambda,type*kg)
      glambda(Syn,lambda,kg)
      glambda(Tg,lambda,kg)

      Glcall(Pl)
      Glcall(Mn)
      Glcall(Demand)
      Glcall(Supply)
      Glcall(Syn)
      Glcall(Tg)
      Glcall(Wind)
      Flcall(Mot)

      if type, Gkcall(Supply), end

      if noSup
        glambda(PV,lambda,type*kg)
        glambda(SW,lambda,kg)
        greactive(PV)
        if type, Gkcall(PV), end
        Glcall(PV)
        Gyreactive(PV)
        Glcall(SW)
      else
        gcall(PV);
        Gycall(PV)
        glambda(SW,1,kg)
      end

      Glreac(PV)
      Fxcall(PV)
      Gkcall(SW)
      greactive(SW)
      Gyreactive(SW)
      Glreac(SW)
      Fxcall(SW,'onlyq')
      
      if perp*iterazione
        Cinc = sigma_corr*inc;
        %cont_eq = Cinc'*([DAE.x;DAE.y;kg;lambda]- ...
        %  [x_old; y_old; kg_old; lambda_old]-Cinc);
        Cinc(end-1) = 0;
        cont_eq = Cinc'*([DAE.x;DAE.y;0;lambda]- ...
                         [x_old; y_old; 0; lambda_old]-Cinc);
        inc_corr = -[DAE.Fx, DAE.Fy, DAE.Fk, DAE.Fl; DAE.Gx, DAE.Gy, ...
                     DAE.Gk, DAE.Gl; Cinc'; Kjac]\[DAE.f; DAE.g; ...
                            cont_eq; DAE.y(SW.refbus)];
        if strcmp(lastwarn,['Matrix is singular to working ' ...
                            'precision.'])
          Cinc(end) = 0;
          cont_eq = Cinc'*([DAE.x;DAE.y;0;0]- ...
                           [x_old; y_old; 0; 0]-Cinc);
          inc_corr = -[DAE.Fx, DAE.Fy, DAE.Fk, DAE.Fl; DAE.Gx, DAE.Gy, ...
                       DAE.Gk, DAE.Gl; Cinc'; Kjac]\[DAE.f; DAE.g; ...
                              cont_eq; DAE.y(SW.refbus)];
        end
      else
        if iterazione
          cont_eq = DAE.y(PQvdx)-sigma_corr*inc(PQvdx+DAE.n)-y_old(PQvdx);
        else
          cont_eq = lambda - sigma_corr*d_lambda - lambda_old;
        end
        inc_corr = -[DAE.Fx, DAE.Fy, DAE.Fk, DAE.Fl; DAE.Gx, DAE.Gy, ...
          DAE.Gk, DAE.Gl; Ljac; Kjac]\[DAE.f; DAE.g; ...
          cont_eq; DAE.y(SW.refbus)];
      end

      DAE.x = DAE.x + inc_corr(1:DAE.n);
      DAE.y = DAE.y + inc_corr(1+DAE.n:DAE.m+DAE.n);
      kg = kg + inc_corr(end-1);
      lambda = lambda + inc_corr(end);
      iter_corr = iter_corr + 1;
      err_corr = max(abs(inc_corr));
    end

    % Generator reactive power computations
    if qlim
      DAE.g = zeros(DAE.m,1);
      fm_call('load');
      glambda(Demand,lambda)
      Bus.Ql = DAE.g(Bus.v);
      fm_call('gen');
      glambda(Demand,lambda)
      Bus.Qg = DAE.g(Bus.v);
      DAE.g = zeros(DAE.m,1);