fm_cpf.m

一个较好的MATLAB潮流程序

	break
      end
    end

    if ~isempty(Vmin_idx) & vlim
      Vmin_idx = Vmin_idx(1);
      PQVmin(Vmin_idx) = 0;
      fm_disp([sp,'PQ load at bus #',fvar(PQ.bus(Vmin_idx),4), ...
	       ' reached V_min at lambda = ',fvar(lambda-one,9)])
      sigma_corr = 1;
      proceed = 1;
      if stop > 1
	proceed = 1;
	sigma_corr = -1;
	break
      end
    end

    % determination of the initial loading factor in case
    % of infeasible underexcitation of generator at zero
    % load condition
    if iterazione == 0 & qlim & ~isempty(Qmin_idx)
      lambda_old = lambda_old + 0.1;
      fm_disp([sp,'Initial loading parameter changed to ', ...
	       fvar(lambda_old-one,4)])
      proceed = 0;
      break
    end

    % Check for Hopf Bifurcations
    if DAE.n >= 2 & CPF.hopf
      As = DAE.Fx-DAE.Fy*(DAE.Jlfv\DAE.Gx);
      if Settings.octave
        auto = eig(As);
      else
        opt.disp = 0;
        auto = eigs(As,2,'SR',opt);
      end
      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'];
      proceed = 0;
      break
    end

    % check lambda step
    if abs(lambda-lambda_old) > 10*abs(d_lambda) & perp
      disp([sp,'corrector solution is too far from predictor increment'])
      proceed = 0;
      break
    end

    if lambda > lambda_old & lambda < max(l_vect)
      disp([sp,'corrector goes back increasing lambda'])
      proceed = 0;
      break
    end

    anylimit = ~(isempty(Qmax_idx) & isempty(Qmin_idx) & ...
                 isempty(Qswmax_idx) & isempty(Qswmin_idx) & ...
                 isempty(Vmax_idx) & isempty(Vmin_idx) & ...
                 isempty(Iij_idx) & isempty(Iji_idx)) & 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.01
      proceed = 0;
      break
    elseif ~isempty(Qmax_idx) & qlim
      Qmax_idx = Qmax_idx(1);
      PV.con(Qmax_idx,[4 6]) = -PV.con(Qmax_idx,[4 6]);
      a = find(PQ.bus == PV.bus(Qmax_idx));
      if isempty(a)
        PQ.con = [PQ.con; [PV.con(Qmax_idx,[1 2 3 4 6 8 9]), 0]];
        PQ.n = PQ.n + 1;
        PQ.bus = [PQ.bus;PV.bus(Qmax_idx)];
        PQVmax = [PQVmax; 1];
        PQVmin = [PQVmin; 1];
      else
        PQ.con(a,[4 5]) = PQ.con(a,[4 5]) + PV.con(Qmax_idx,[4 6]);
        PQ.con(a,6) = min(PQ.con(a,6),PV.con(Qmax_idx,8));
        PQ.con(a,7) = max(PQ.con(a,7),PV.con(Qmax_idx,9));
      end
      PV.n = PV.n - 1;
      PV.con(Qmax_idx,:) = [];
      fm_disp([sp,'PV gen. at bus #', fvar(PV.bus(Qmax_idx),4), ...
	       ' reached Q_max at lambda = ',fvar(lambda-one,9)],1)
      PV.bus(Qmax_idx) = [];
      Qmax_idx = [];
      proceed = 0;
      break
    elseif ~isempty(Qmin_idx) & qlim
      Qmin_idx = Qmin_idx(1);
      PV.con(Qmin_idx,[4 7]) = -PV.con(Qmin_idx,[4 7]);
      a = find(PQ.bus == PV.bus(Qmin_idx));
      if isempty(a)
        PQ.con = [PQ.con; [PV.con(Qmin_idx,[1 2 3 4 7 8 9]), 0]];
        PQ.n = PQ.n + 1;
        PQ.bus = [PQ.bus;PV.bus(Qmin_idx)];
        PQVmax = [PQVmax; 1];
        PQVmin = [PQVmin; 1];
      else
        PQ.con(a,[4 5]) = PQ.con(a,[4 5]) + PV.con(Qmin_idx,[4 7]);
        PQ.con(a,6) = min(PQ.con(a,6),PV.con(Qmin_idx,8));
        PQ.con(a,7) = max(PQ.con(a,7),PV.con(Qmin_idx,9));
      end
      PV.n = PV.n - 1;
      PV.con(Qmin_idx,:) = [];
      fm_disp([sp,'PV gen. at bus #', fvar(PV.bus(Qmin_idx)',4), ...
	       ' reached Q_min at lambda = ',fvar(lambda-one,9)],1)
      PV.bus(Qmin_idx) = [];
      Qmin_idx = [];
      proceed = 0;
      break
    elseif ~isempty(Qswmax_idx) & qlim & type
      Qswmax_idx = Qswmax_idx(1);
      SW.con(Qswmax_idx,[10 6]) = -SW.con(Qswmax_idx,[10 6]);
      a = find(PQ.bus == SW.bus(Qswmax_idx));
      if isempty(a)
        PQ.con = [PQ.con; [SW.con(Qswmax_idx,[1 2 3 10 6 8 9]), 0]];
        PQ.n = PQ.n + 1;
        PQ.bus = [PQ.bus;SW.bus(Qswmax_idx)];
        PQVmax = [PQVmax; 1];
        PQVmin = [PQVmin; 1];
      else
        PQ.con(a,[4 5]) = PQ.con(a,[4 5]) + SW.con(Qswmax_idx,[10 6]);
        PQ.con(a,6) = min(PQ.con(a,6),SW.con(Qswmax_idx,8));
        PQ.con(a,7) = max(PQ.con(a,7),SW.con(Qswmax_idx,9));
      end
      SW.n = SW.n - 1;
      SW.con(Qswmax_idx,:) = [];
      fm_disp([sp,'SW gen. at bus #', fvar(SW.bus(Qswmax_idx)',4), ...
	       ' reached Q_max at lambda = ',fvar(lambda-one,9)],1)
      SW.bus(Qswmax_idx) = [];
      Qswmax_idx = [];
      proceed = 0;
      break
    elseif ~isempty(Qswmin_idx) & qlim & type
      Qswmin_idx = Qswmin_idx(1);
      SW.con(Qswmin_idx,[10 7]) = -SW.con(Qswmin_idx,[10 7]);
      a = find(PQ.bus == SW.bus(Qswmin_idx));
      if isempty(a)
        PQ.con = [PQ.con; [SW.con(Qswmin_idx,[1 2 3 10 7 8 9]), 0]];
        PQ.n = PQ.n + 1;
        PQ.bus = [PQ.bus;SW.bus(Qswmin_idx)];
        PQVmax = [PQVmax; 1];
        PQVmin = [PQVmin; 1];
      else
        PQ.con(a,[4 5]) = PQ.con(a,[4 5]) + SW.con(Qswmin_idx,[10 7]);
        PQ.con(a,6) = min(PQ.con(a,6),SW.con(Qswmin_idx,8));
        PQ.con(a,7) = max(PQ.con(a,7),SW.con(Qswmin_idx,9));
      end
      SW.n = SW.n - 1;
      SW.con(Qswmin_idx,:) = [];
      fm_disp([sp,'SW gen. at bus #', fvar(SW.bus(Qswmin_idx)',4), ...
	       ' reached Q_min at lambda = ',fvar(lambda-one,9)],1)
      SW.bus(Qswmin_idx) = [];
      Qswmin_idx = [];
      proceed = 0;
      break
    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
    if lambda < 0 & iterazione > 1
      fm_disp([sp,'lambda < 0 at iter. ',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;
    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(PQidxplot)],iterazione)

    if sigma_corr < tol, break, end
    sigma_corr = 1;
  else
    if abs(lambda-CPF.linit) < tol & iterazione > 1
      break
    end
    DAE.V = V_old;
    DAE.a = ang_old;
    DAE.x = x_old;
    kg = kg_old;
    lambda = lambda_old;
    sigma_corr = 0.5*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
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

% Sensitivity Coefficients
% ------------------------------------------------------------------
if Settings.octave
  Dpfc = zeros(2*Bus.n+1,Demand.n+Supply.n);
else
  Dpfc = sparse(2*Bus.n+1,Demand.n+Supply.n);
end
for i = 1:Supply.n
  k = Supply.bus(i);
  Dpfc(k,i) = -lambda;
end
for i = 1:Demand.n
  k = Demand.bus(i);
  Dpfc(k,Supply.n+i) = lambda;
  Dpfc(end,i) = kg;
end

try
  dl_dp = full(-W'*Dpfc/(W'*[Gl;0]))';
catch
  dl_dp = NaN;
end

%  Visualization of results
% --------------------------------------------------------------
if nodyn, DAE.n = 0; end
Settings.lftime = toc;
Settings.iter = iterazione;
DAE.V = Varout.V(idx_crit,:)';
DAE.a = Varout.ang(idx_crit,:)';
if DAE.n, DAE.x = Varout.x(idx_crit,:)'; end
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');
  for i = 1:Demand.n
    k = Demand.bus(i);
    DAE.gq(k) = lambda*Demand.con(i,4) + DAE.gq(k);
    DAE.gp(k) = lambda*Demand.con(i,3) + DAE.gp(k);
  end
  % 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
end

%  Reset of SW, PV and PQ structures
%if ~strcmp(fun,'snb')  
PQ = PQold;
PV = PVold;
SW = SWold;
if noDem
  Demand.con = [];
  Demand.n = 0;
  Demand.bus = [];
end
if noSup
  Supply.con = [];
  Supply.n = 0;
  Supply.bus = [];
end
if Syn.n
  Syn.pm = Pm_syn;
end
%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
CPF.lambda = lambda_crit;
CPF.kg = kg;

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