fm_cpf.m

基于PSAT 软件的多目标最优潮流计算用于中小型电力系统的分析和管理

      DAE.V = DAE.V + inc_corr(DAE.n+Bus.n+1: DAE.n+2*Bus.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.gq = zeros(Bus.n,1);
      fm_call('pq');
      glambda(Demand,lambda)
      Bus.Ql = DAE.gq;
      Q = DAE.glfq;
      Bus.Qg = Q + Bus.Ql;
      [Qmax_idx,Qmin_idx] = pvlim(PV,Bus.Qg(PV.bus));
      [Qswmax_idx,Qswmin_idx] = swlim(SW,Bus.Qg(SW.bus));
      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
      VV = DAE.V.*exp(jay*DAE.a);
      switch ilim
       case 1
        Fij = abs(((VV(Line.from) - tps.*VV(Line.to)).*y + ...
		   VV(Line.from).*(jay*chrg))./(tps.*conj(tps)));
        Fji = abs((VV(Line.to) - VV(Line.from)./tps).*y + ...
		  VV(Line.to).*(jay*chrg));
       case 2
        Fij = abs(real(VV(Line.from).*conj(((VV(Line.from) - ...
	      tps.*VV(Line.to)).*y + ...
	      VV(Line.from).*(jay*chrg))./(tps.*conj(tps)))));
        Fji = abs(real(VV(Line.to).*conj((VV(Line.to) - ...
	      VV(Line.from)./tps).*y + VV(Line.to).*(jay*chrg))));
       case 3
        Fij = abs(VV(Line.from).*conj(((VV(Line.from) - ...
	      tps.*VV(Line.to)).*y + VV(Line.from).* ...
	      (jay*chrg))./(tps.*conj(tps))));
        Fji = abs(VV(Line.to).*conj((VV(Line.to) - ...
	      VV(Line.from)./tps).*y + VV(Line.to).*(jay*chrg)));
      end
      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
      As = DAE.Fx-DAE.Fy*(DAE.Jlfv\DAE.Gx);
      opt.disp = 0;
      auto = eigs(As,2,'SR',opt);
      auto = round(auto/Settings.lftol)*Settings.lftol;
      hopf_idx = find(real(auto) > 0);
      if ~isempty(hopf_idx)
        if imag(auto(hopf_idx)) ~= 0
          fm_disp([sp,'Hopf bifurcation encountered.'])
        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.from(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.from(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
      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,
        fm_disp([sp,'Saddle-node Bifurcation or Limit ', ...
		 'Induced Bifurcation encountered.'])
        sigma_corr = -1;
      end
      break
    end
  end

  if proceed == 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;
    %etaR(iterazione) = min(diag(DAE.J22-DAE.J21*(DAE.J11\DAE.J12)));
    if Syn.n
      Syn.pm(syn_idx) = (lambda+kg)*Pm_syn(syn_idx);
    end
    fm_out(2,lambda,iterazione)
    if Syn.n
      Syn.pm(syn_idx) = Pm_syn(syn_idx);
    end

    fm_status('cpf','update',[lambda, DAE.V(PQidx)],iterazione)

    if sigma_corr < tol, break, end
    sigma_corr = 1; 
    if lambda > lambda_old
      V_crit = DAE.V;
      a_crit = DAE.a;
      x_crit = DAE.x;
      k_crit = kg;
      l_crit = lambda;
    end
  elseif proceed == 0
    if abs(lambda-CPF.linit) < 0.001 & ...
          abs(lambda-lambda_old) <= abs(d_lambda) & ...
          iterazione > 1
      fm_disp([sp,'d_lambda = ',num2str([d_lambda, lambda])])
      break
    end
    DAE.V = V_old;
    DAE.a = ang_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.J11 = sparse(Bus.n,Bus.n);
  DAE.J21 = sparse(Bus.n,Bus.n);
  DAE.J12 = sparse(Bus.n,Bus.n);
  DAE.J22 = sparse(Bus.n,Bus.n);
  fm_lf(2);
  if noSup
    Gyreactive(PV)
  else
    Gycall(PV);
  end
  Gyreactive(SW)
  DAE.Jlfv = [DAE.J11, DAE.J12; DAE.J21, DAE.J22];
  if Bus.n > 250
    [L,U,P] = luinc([DAE.Fx,DAE.Fy,Fk;DAE.Gx,DAE.Jlfv,DAE.Gk;kjac],1e-6);
  else
    [L,U,P] = lu([DAE.Fx,DAE.Fy,Fk;DAE.Gx,DAE.Jlfv,DAE.Gk;kjac]);
  end
  dz_dl = -U\(L\(P*[Flambda;DAE.Gl;0]));
  Jsign_old = Jsign;
  Jsign = signperm(P)*sign(prod(diag(U)));
  if Jsign ~= Jsign_old & Sflag & iterazione > 1
    ksign = -1;
    Sflag = 0;
  end
  Norm = norm(dz_dl,inf);
  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+Bus.n+PQidx);
    a = min(a,0.025);
    a = max(a,-0.025);
    inc(DAE.n+Bus.n+PQidx) = a;
    Ljac(end) = 0;
    Ljac(DAE.n+Bus.n+PQidx) = 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.V = V_crit;
DAE.a = a_crit;
DAE.x = x_crit;
kg = k_crit;
lambda = l_crit;

if ~isempty(DAE.V)
  if CPF.show
    fm_disp(['Continuation Power Flow completed in ', ...
	     num2str(toc),' s'],1);
  end
  maxgqerr = max(abs(DAE.gq));
  maxgperr = max(abs(DAE.gp));
  DAE.gq = zeros(Bus.n,1);
  DAE.gp = zeros(Bus.n,1);
  fm_call('pq');
  glambda(Demand,lambda)
  % load powers
  Bus.Pl = DAE.gp;
  Bus.Ql = DAE.gq;
  % total bus injections
  Vc = DAE.V.*exp(jay*DAE.a);
  S = Vc.*conj(Line.Y*Vc);
  % generated powers
  Bus.Qg = imag(S) + Bus.Ql;
  Bus.Pg = real(S) + Bus.Pl;
  DAE.gq = maxgqerr*ones(Bus.n,1);
  DAE.gp = maxgperr*ones(Bus.n,1);
  if (Settings.showlf == 1 & CPF.show) | Fig.stat
    SDbus = [Supply.bus;Demand.bus];
    report = cell(2+length(dl_dp),1);
    report{1,1} = ['Lambda_max = ', fvar(lambda_crit,9)];
    report{2,1} = ' ';
    for i = 1:length(dl_dp)
      report{2+i,1} = ['d lambda / d P ', ...
		       Varname.bus{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
PQ = restore(PQ);
PV = restore(PV);
SW = restore(SW);
if noDem, Demand = restore(Demand); end
if noSup, Supply = restore(Supply); end

if Syn.n
  Syn.pm = Pm_syn;
end
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')

%warning on
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);