fm_cpf.m

用于电力系统的一个很好的分析软件

      [Qmax_idx,Qmin_idx] = pvlim(PV);
      [Qswmax_idx,Qswmin_idx] = swlim(SW);
      if ~kqmin
        Qmin_idx = [];
        Qswmin_idx = [];
      end
    end

    [PQ,lim_v] = pqlim(PQ,vlim,sp,lambda,one);

    if lim_v
      sigma_corr = 1;
      proceed = 1;
      if stop > 1
        sigma_corr = -1;
        break
      end
    end

    if ilim
      [Fij,Fji] = flows(Line,ilim);
      Iij_idx = find(Fij > Imax & Iij);
      Iji_idx = find(Fji > Imax & Iji);
    end

    % determination of the initial loading factor in case
    % of infeasible underexcitation of generator at zero
    % load condition
    if ~iterazione & qlim & ~isempty(Qmin_idx) & count_qmin <= 5
      count_qmin = count_qmin+1;
      if count_qmin > 5
        fm_disp([sp,'There are generator Qmin limit violations at ' ...
          'the initial point.'])
        fm_disp([sp,'Generator Qmin limits will be disabled.'])
        kqmin = 0;
        lambda_old = CPF.linit;
        lambda = CPF.linit;
        sigma_corr = 1;
      else
        lambda_old = lambda_old + d_lambda;
        fm_disp([sp,'Initial loading parameter changed to ', ...
          fvar(lambda_old-one,4)])
      end
      proceed = 0;
      break
    end

    % Check for Hopf Bifurcations
    if DAE.n >= 2 & CPF.hopf & hopf
      As = DAE.Fx-DAE.Fy*(DAE.Gy\DAE.Gx);
      if DAE.n > 100
        opt.disp = 0;
        auto = eigs(As,20,'SR',opt);
      else
        auto = eig(full(As));
      end
      auto = round(auto/Settings.lftol)*Settings.lftol;
      hopf_idx = find(real(auto) > 0);
      if ~isempty(hopf_idx)
        hopf = 0;
        hopf_idx = find(abs(imag(auto(hopf_idx))) > 1e-5);
        if ~isempty(hopf_idx)
          lim_hb = 1;
          fm_disp([sp,'Hopf bifurcation encountered.'])
        end
        if stop
          sigma_corr = -1;
          break
        end
      end
    end

    if ~isempty(Iij_idx) & ilim
      Iij_idx = Iij_idx(1);
      Iij(Iij_idx) = 0;
      fm_disp([sp,fw,'from bus #',fvar(Line.fr(Iij_idx),4), ...
        ' to bus #',fvar(Line.to(Iij_idx),4), ...
        ' reached I_max at lambda = ',fvar(lambda-one,9)],1)
      sigma_corr = 1;
      proceed = 1;
      if stop > 1
        proceed = 1;
        sigma_corr = -1;
        break
      end
    end

    if ~isempty(Iji_idx) & ilim
      Iji_idx = Iji_idx(1);
      Iji(Iji_idx) = 0;
      fm_disp([sp,fw,'from bus #',fvar(Line.to(Iji_idx),4), ...
        ' to bus #',fvar(Line.fr(Iji_idx),4), ...
        ' reached I_max at lambda = ',fvar(lambda-one,9)],1)
      sigma_corr = 1;
      proceed = 1;
      if stop > 1
        proceed = 1;
        sigma_corr = -1;
        break
      end
    end

    if lambda < CPF.linit
      cpfmsg = [sp,'lambda is lower than initial value'];
      if iterazione > 5
        proceed = 0;
        break
      end
    end

    if abs(lambda-lambda_old) > 5*abs(d_lambda) & perp & iterazione
      fm_disp([sp,'corrector solution is too far from predictor value'])
      proceed = 0;
      break
    end
    if lambda > lambda_old & lambda < max(l_vect)
      fm_disp([sp,'corrector goes back increasing lambda'])
      proceed = 0;
      break
    end

    lim_q = (~isempty(Qmax_idx) | ~isempty(Qmin_idx) | ...
      (~isempty(Qswmax_idx) | ~isempty(Qswmin_idx))*type) & qlim;
    lim_i = (~isempty(Iij_idx) | ~isempty(Iji_idx)) & ilim;

    anylimit = (lim_q | lim_v | lim_i) & CPF.stepcut;

    if iter_corr > iter_max % no convergence
      if lambda_old < 0.5*max(l_vect)
        fm_disp([sp,'Convergence problems in the unstable curve'])
        lambda = -1;
        proceed = 1;
        break
      end
      if sigma_corr < 1e-3
        proceed = 1;
        break
      end
      if CPF.show
        fm_disp([sp,'Max # of iters. at corrector step.'])
        fm_disp([sp,'Reduced Variable Increments in ', ...
          'Corrector Step ', num2str(0.5*sigma_corr)])
      end
      cpfmsg = [sp,'Reached maximum number of iterations ', ...
        'for corrector step.'];
      proceed = 0;
      break
    elseif anylimit & sigma_corr > 0.11
      proceed = 0;
      break
    elseif lim_q
      if ~isempty(Qmax_idx)
        PQ = add(PQ,pqdata(PV,Qmax_idx(1),'qmax',sp,lambda,one,noSup));
        if noSup,
          PV = move2sup(PV,Qmax_idx(1));
        else
          PV = remove(PV,Qmax_idx(1));
        end
      elseif ~isempty(Qmin_idx)
        PQ = add(PQ,pqdata(PV,Qmin_idx(1),'qmin',sp,lambda,one,noSup));
        if noSup
          PV = move2sup(PV,Qmin_idx(1));
        else
          PV = remove(PV,Qmin_idx(1));
        end
      elseif ~isempty(Qswmax_idx)
        PQ = add(PQ,pqdata(SW,Qswmax_idx(1),'qmax',sp,lambda,one));
        SW = remove(SW,Qswmax_idx(1));
      elseif ~isempty(Qswmin_idx)
        PQ = add(PQ,pqdata(SW,Qswmin_idx(1),'qmin',sp,lambda,one));
        SW = remove(SW,Qswmin_idx(1));
      end
      lim_lib = 1;
      Qmax_idx = [];
      Qmin_idx = [];
      Qswmax_idx = [];
      Qswmin_idx = [];
      if ~iterazione
        d_lambda = 0;
        lambda_old = CPF.linit;
        lambda = CPF.linit;
        if perp
          inc(end) = 1e-5;
        else
          Ljac(end) = 1;
        end
      else
        lambda = lambda_old;
        sigma_corr = 1;
        proceed = 0;
        break
      end
    else
      proceed = 1;
      sigma_corr = 1;
      if stop & ksign < 0,
        if lim_lib
          fm_disp([sp,'Saddle-Node Bifurcation encountered.'])
        else
          fm_disp([sp,'Limit-Induced Bifurcation encountered.'])
        end
        sigma_corr = -1;
      end
      break
    end
  end

  switch proceed 
   case 1
    if lambda < 0 & iterazione > 1
      fm_disp([sp,'lambda < 0 at iteration ',num2str(iterazione)])
      break
    end
    l_vect = [l_vect; lambda];

    if CPF.show
      fm_disp(['Iteration = ',fvar(iterazione+1,5),'lambda =', ...
        fvar(lambda-one,9), '    kg =',fvar(kg,9)],1)
    end
    iterazione = iterazione + 1;

    fm_out(2,lambda,iterazione)
    
    fm_status('cpf','update',[lambda, DAE.y(PQvdx)],iterazione)

    if sigma_corr < tol, break, end
    sigma_corr = 1;
    if lambda > lambda_old
      y_crit = DAE.y;
      x_crit = DAE.x;
      k_crit = kg;
      l_crit = lambda;
    end
   case 0
    if abs(lambda-CPF.linit) < 0.001 & ...
          abs(lambda-lambda_old) <= 10*abs(d_lambda) & ...
          iterazione > 1
      %fm_disp([sp,'d_lambda = ',num2str([d_lambda, lambda])])
      break
    end
    DAE.y = y_old;
    DAE.x = x_old;
    kg = kg_old;
    lambda = lambda_old;
    sigma_corr = 0.1*sigma_corr;
    if sigma_corr < tol,
      if ~isempty(cpfmsg)
        fm_disp(cpfmsg)
      end
      if iterazione == 0
        fm_disp([sp,'Infeasible initial loading condition.'])
      else
        fm_disp([sp,'Convergence problem encountered.'])
      end
      break
    end
  end

  % stop routine
  % --------------------------------------------------------------

  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

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

  DAE.lambda = lambda;

  % update Jacobians

  fm_call('kg')

  if nodyn, DAE.Fx = 1; end

  if noSup
    Gyreactive(PV)
  else
    Gycall(PV);
  end
  Gyreactive(SW)

  if Bus.n > 250
    [L,U,P] = luinc([DAE.Fx,DAE.Fy,DAE.Fk;DAE.Gx,DAE.Gy,DAE.Gk;kjac],1e-6);
  else
    [L,U,P] = lu([DAE.Fx,DAE.Fy,DAE.Fk;DAE.Gx,DAE.Gy,DAE.Gk;kjac]);
  end
  dz_dl = -U\(L\(P*[DAE.Fl;DAE.Gl;0]));
  Jsign_old = Jsign;
  Jsign = signperm(P)*sign(prod(diag(U)));

  if lim_lib
    fm_snap('assignsnap','new','LIB',lambda)
  elseif lim_hb
    fm_snap('assignsnap','new','HB',lambda)
  end

  if iterazione == 1
    if noDem & lambda == 1
      fm_snap('assignsnap','start','OP',lambda)
    elseif ~noDem & lambda == 0
      fm_snap('assignsnap','start','OP',lambda)
    else
      fm_snap('assignsnap','start','Init',lambda)
    end
  end

  if Jsign ~= Jsign_old & Sflag & iterazione > 1
    ksign = -1;
    Sflag = 0;
    if ~lim_lib, fm_snap('assignsnap','new','SNB',lambda), end
  end

  Norm = norm(dz_dl,2);
  if Norm == 0, Norm = 1; end
  d_lambda = ksign*CPF.step/Norm;
  d_lambda = min(d_lambda,0.35);
  d_lambda = max(d_lambda,-0.35);
  if ksign > 0, d_lambda = max(d_lambda,0); end
  if ksign < 0, d_lambda = min(d_lambda,0); end
  d_z = d_lambda*dz_dl;
  inc = [d_z; d_lambda];
  if ~perp & iterazione
    a = inc(DAE.n+PQvdx);
    a = min(a,0.025);
    a = max(a,-0.025);
    inc(DAE.n+PQvdx) = a;
    Ljac(end) = 0;
    Ljac(DAE.n+PQvdx) = 1;
  end

end

fm_out(3,0,iterazione)
[lambda_crit, idx_crit] = max(Varout.t);
if isnan(lambda_crit), lambda_crit = lambda; end
if isempty(lambda_crit), lambda_crit = lambda; end
if CPF.show
  fm_disp([sp,'Maximum Loading Parameter lambda_max = ', ...
    fvar(lambda_crit-one,9)],1)
end

%  Visualization of results
% --------------------------------------------------------------
if nodyn
  DAE.n = 0;
  Varout.idx = Varout.idx-1;
end
Settings.lftime = toc;
Settings.iter = iterazione;

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

if ~isempty(DAE.y)
  if CPF.show
    fm_disp(['Continuation Power Flow completed in ', ...
      num2str(toc),' s'],1);
  end
  DAE.g = zeros(DAE.m,1);
  fm_call('load');
  glambda(Demand,lambda)
  % load powers
  Bus.Pl = DAE.g(Bus.a);
  Bus.Ql = DAE.g(Bus.v);
  % gen powers
  fm_call('gen')
  Bus.Qg = DAE.g(Bus.a);
  Bus.Pg = DAE.g(Bus.v);
  DAE.g = err_corr*ones(DAE.m,1);
  if (Settings.showlf == 1 & CPF.show) | Fig.stat
    SDbus = [Supply.bus;Demand.bus];
    report = cell(1,1);
    report{1,1} = ['Lambda_max = ', fvar(lambda_crit,9)];
    %for i = 1:length(dl_dp)
    %  report{2+i,1} = ['d lambda / d P ', ...
    %		       Bus.names{SDbus(i)},' = ', ...
    %		       fvar(dl_dp(i),9)];
    %end
    fm_stat(report);
  end
  if CPF.show & Fig.plot
    fm_plotfig
  end
end

%  Reset of SW, PV and PQ structures
Settings.forcepq = forcepq;
PQ = restore(PQ);
PV = restore(PV);
SW = restore(SW);
Demand = restore(Demand);
Supply = restore(Supply);

if CPF.show & Fig.main
  set(Fig.main,'Pointer','arrow');
  Settings.xlabel = ['Loading Parameter ',char(92),'lambda (p.u.)'];
  Settings.tf = 1.2*lambda_crit;
end

fm_status('cpf','close')

Settings.pv2pq = PV2PQ;
CPF.lambda = lambda_crit;
CPF.kg = kg;

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

% --------------------------------
function s = signperm(P)
% --------------------------------

[i,j,p] = find(sparse(P));
idx = find(i ~= j);
q = P(idx,idx);
s = det(q);