fm_cpf.m

电力系统的psat

    DAE.gp = zeros(Bus.n,1);
    DAE.gq = zeros(Bus.n,1);
    if Settings.octave
      DAE.J11 = zeros(Bus.n,Bus.n);
      DAE.J21 = zeros(Bus.n,Bus.n);
      DAE.J12 = zeros(Bus.n,Bus.n);
      DAE.J22 = zeros(Bus.n,Bus.n);
    else
      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);
    end
  end
  if Settings.octave
    Gl = zeros(2*Bus.n,1);
    Gk = zeros(2*Bus.n,1);
  else
    Gl = sparse(2*Bus.n,1);
    Gk = sparse(2*Bus.n,1);
  end

  % component calls
  fm_call('kg');

  if SW.n+PV.n
    Gbus = Bus.n+[SW.bus;PV.bus];
    DAE.Jlfv(Gbus,:) = 0;
    DAE.Jlfv(:,Gbus) = 0;
    if Settings.octave
      DAE.Jlfv(Gbus,Gbus) = eye(PV.n+SW.n);
    else
      DAE.Jlfv(Gbus,Gbus) = speye(PV.n+SW.n);
    end
  end
  if SW.n
    DAE.Jlfv(:,SW.bus) = 0;
    DAE.Jlfv(SW.bus,:) = 0;
    if Settings.octave
      DAE.Jlfv(SW.bus,SW.bus) = eye(SW.n);
    else
      DAE.Jlfv(SW.bus,SW.bus) = speye(SW.n);
    end
    DAE.Fy(:,SW.bus) = 0;
    DAE.Gx(SW.bus,:) = 0;
    DAE.Fy(:,Bus.n+SW.bus) = 0;
    DAE.Gx(Bus.n+SW.bus,:) = 0;
  end
  if PV.n
    DAE.Fy(:,Bus.n+PV.bus) = 0;
    DAE.Gx(Bus.n+PV.bus,:) = 0;
  end
  if nodyn, DAE.Fx = 1; end

  % in case of floating power supply, base case power is
  % set for the reference bus
  if Supply.n
    for i = 1:SW.n
      k = SW.bus(i);
      DAE.gp(k) = DAE.gp(k) - SW.con(i,10);
    end
  end
  for i = 1:Demand.n
    k = Demand.bus(i);
    DAE.gp(k) = lambda*Demand.con(i,3) + DAE.gp(k);
    DAE.gq(k) = lambda*Demand.con(i,4) + DAE.gq(k);
    Gl(k,1) = Gl(k,1)+Demand.con(i,3);
    Gl(k+Bus.n,1) = Gl(k+Bus.n,1)+Demand.con(i,4);
  end
  for i = 1:Supply.n
    k = Supply.bus(i);
    DAE.gp(k) = DAE.gp(k) - lambda*Supply.con(i,3);
    Gl(k,1) = Gl(k,1)-Supply.con(i,3);
  end
  if Supply.n
    if type % distributed slack bus model
      %DAE.Jlfv(SW.bus,:) = DAE.Jlf(SW.bus,:);
      %DAE.Jlfv(:,SW.bus) = DAE.Jlf(:,SW.bus);
      for i = 1:Supply.n
        k = Supply.bus(i);
        DAE.gp(k) = DAE.gp(k) - kg*ksu(i)*Supply.con(i,3);
        Gk(k,1) = Gk(k,1)-ksu(i)*Supply.con(i,3);
      end
    else    % single slack bus model
      for i = 1:SW.n
        k = SW.bus(i);
        Psup = sum(Supply.con(find(Supply.bus == k),3));
        DAE.gp(k) = DAE.gp(k) - kg*ksw(i)*Psup;
        if Psup == 0,
          Psup = SW.con(i,10);
        end
        Gk(k,1) = -ksw(i)*Psup;
        DAE.Jlfv(SW.bus,:) = 0;
        if Settings.octave
          DAE.Jlfv(SW.bus,SW.bus) = eye(SW.n);
        else
          DAE.Jlfv(SW.bus,SW.bus) = speye(SW.n);
        end
      end
    end
  end
  if Syn.n
    DAE.f(pm_idx) = DAE.f(pm_idx)+ ...
        0.5*(lambda+kg-1)*Pm_syn(syn_idx)./Syn.con(syn_idx,18);
    Flambda(pm_idx) = 0.5*Pm_syn(syn_idx)./Syn.con(syn_idx,18);
    Fk(pm_idx) = 0.5*Pm_syn(syn_idx)./Syn.con(syn_idx,18);
    DAE.Jlfv(Settings.refbus,:) = 0;
    DAE.Jlfv(:,Settings.refbus) = 0;
    DAE.Jlfv(Settings.refbus,Settings.refbus) = 1;
    DAE.gp(Settings.refbus) = 0;
  end
  if PV.n
    DAE.gq(PV.bus) = 0;
    Gl(PV.bus+Bus.n) = 0;
  end
  if SW.n
    DAE.gq(SW.bus) = 0;
    Gl(SW.bus+Bus.n) = 0;
    Gl(SW.bus) = 0;
  end

  % predictor step
  % --------------------------------------------------------------
  if iterazione > 0
    [L,U,P] = lu([DAE.Fx,DAE.Fy,Fk;DAE.Gx,DAE.Jlfv,Gk; kjac]);
    dz_dl = -U\(L\(P*[Flambda;Gl;0]));
    Jsign_old = Jsign;
    Jsign = sign(prod(diag(U))*det(P));
    if Jsign ~= Jsign_old & Sflag
      ksign = -1; 
      Sflag = 0;
    end
    Norm = norm(dz_dl);
    if Norm == 0, Norm = 1; end
    d_lambda = ksign*CPF.step/Norm;
    d_lambda = min(d_lambda,0.5);
    d_lambda = max(d_lambda,-0.5);
    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;
    if lambda > lambda_old, W = d_z(DAE.n+1:end); end
    inc = [d_z; d_lambda];
    if ~perp
      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
  else
    W = ones(2*Bus.n+1,1);
    Jsign = sign(det([DAE.Fx,DAE.Fy,Fk;DAE.Gx,DAE.Jlfv,Gk; kjac]));
    d_lambda = 1e-5;
    lambda_old = CPF.linit;
    inc = zeros(DAE.n+2*Bus.n+2,1);
    if perp
      inc(end) = 1e-5;
    else
      d_lambda = 1e-5;
      Ljac(end) = 1;
    end
  end

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

  V_old = DAE.V;
  ang_old = DAE.a;
  x_old = DAE.x;
  kg_old = kg;
  lambda_old = lambda;
  corrector = 1;
  if PV.n
    inc(DAE.n+Bus.n+PV.bus) = 0;
  end
  if SW.n
    inc(DAE.n+Bus.n+SW.bus) = 0;
    inc(DAE.n+SW.bus) = 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

      if isempty(Line.Y)
        DAE.gp = zeros(Bus.n,1);
        DAE.gq = zeros(Bus.n,1);
        if Settings.octave
          DAE.J11 = zeros(Bus.n,Bus.n);
          DAE.J21 = zeros(Bus.n,Bus.n);
          DAE.J12 = zeros(Bus.n,Bus.n);
          DAE.J22 = zeros(Bus.n,Bus.n);
        else
          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);
        end
      end
      if Settings.octave
        Gl = zeros(2*Bus.n,1);
        Gk = zeros(2*Bus.n,1);
      else
        Gl = sparse(2*Bus.n,1);
        Gk = sparse(2*Bus.n,1);
      end

      % component calls
      fm_call('kg');

      if SW.n+PV.n
        Gbus = Bus.n+[SW.bus;PV.bus];
        DAE.Jlfv(Gbus,:) = 0;
        DAE.Jlfv(:,Gbus) = 0;
        if Settings.octave
          DAE.Jlfv(Gbus,Gbus) = eye(PV.n+SW.n);
        else
          DAE.Jlfv(Gbus,Gbus) = speye(PV.n+SW.n);
        end
      end
      if SW.n
        DAE.Jlfv(:,SW.bus) = 0;
        DAE.Jlfv(SW.bus,:) = 0;
        if Settings.octave
          DAE.Jlfv(SW.bus,SW.bus) = eye(SW.n);
        else
          DAE.Jlfv(SW.bus,SW.bus) = speye(SW.n);
        end
        DAE.Fy(:,SW.bus) = 0;
        DAE.Gx(SW.bus,:) = 0;
        DAE.Fy(:,Bus.n+SW.bus) = 0;
        DAE.Gx(Bus.n+SW.bus,:) = 0;
      end
      if PV.n
        DAE.Fy(:,Bus.n+PV.bus) = 0;
        DAE.Gx(Bus.n+PV.bus,:) = 0;
      end
      if nodyn, DAE.Fx = 1; end

      % In case of floating power supply, base case power
      % is set for the reference bus
      if Supply.n,
        for i = 1:SW.n,
          k = SW.bus(i);
          DAE.gp(k) = DAE.gp(k) - SW.con(i,10);
        end
      end
      for i = 1:Demand.n
        k = Demand.bus(i);
        DAE.gp(k) = lambda*Demand.con(i,3) + DAE.gp(k);
        DAE.gq(k) = lambda*Demand.con(i,4) + DAE.gq(k);
        Gl(k,1) = Gl(k,1)+Demand.con(i,3);
        Gl(k+Bus.n,1) = Gl(k+Bus.n,1)+Demand.con(i,4);
      end
      for i = 1:Supply.n
        k = Supply.bus(i);
        DAE.gp(k) = DAE.gp(k) - lambda*Supply.con(i,3);
        Gl(k,1) = Gl(k,1)-Supply.con(i,3);
      end
      if Supply.n
        if type % distributed slack bus
          %DAE.Jlfv(SW.bus,:) = DAE.Jlf(SW.bus,:);
          %DAE.Jlfv(:,SW.bus) = DAE.Jlf(:,SW.bus);
          for i = 1:Supply.n
            k = Supply.bus(i);
            DAE.gp(k) = DAE.gp(k) - kg*ksu(i)*Supply.con(i,3);
            Gk(k,1) = Gk(k,1)-ksu(i)*Supply.con(i,3);
          end
        else % single slack bus
          for i = 1:SW.n
            k = SW.bus(i);
            Psup = sum(Supply.con(find(Supply.bus == k),3));
            if Psup == 0,
              Psup = SW.con(i,10);
            end
            DAE.gp(k) = DAE.gp(k) - kg*ksw(i)*Psup;
            Gk(k,1) = -ksw(i)*Psup;
            DAE.Jlfv(SW.bus,:) = 0;
            if Settings.octave
              DAE.Jlfv(SW.bus,SW.bus) = eye(SW.n);
            else
              DAE.Jlfv(SW.bus,SW.bus) = speye(SW.n);
            end
          end
        end
      end
      if Syn.n
        DAE.f(pm_idx) = DAE.f(pm_idx)+ ...
            0.5*(lambda+kg-1)*Pm_syn(syn_idx)./Syn.con(syn_idx,18);
        Flambda(pm_idx,1) = 0.5*Pm_syn(syn_idx)./Syn.con(syn_idx,18);
        Fk(pm_idx) = 0.5*Pm_syn(syn_idx)./Syn.con(syn_idx,18);
        DAE.Jlfv(Settings.refbus,:) = 0;
        DAE.Jlfv(:,Settings.refbus) = 0;
        DAE.Jlfv(Settings.refbus,Settings.refbus) = 1;
        DAE.gp(Settings.refbus) = 0;
        %Syn.pm(syn_idx) = (lambda+kg-1)*Pm_syn(syn_idx);
      end

      if PV.n
        DAE.gq(PV.bus) = 0;
        Gl(PV.bus+Bus.n) = 0;
      end
      if SW.n
        DAE.gq(SW.bus) = 0;
        Gl(SW.bus+Bus.n) = 0;
        Gl(SW.bus) = 0;
        DAE.Jlfv(SW.bus(1),SW.bus(1)) = 1;
      end

      if perp
        Cinc = sigma_corr*inc;
        cont_eq = Cinc'*([DAE.x;DAE.a;DAE.V;kg;lambda]- ...
			 [x_old; ang_old; V_old;kg_old; ...
			  lambda_old]-Cinc);
	AC = [DAE.Fx, DAE.Fy, Fk, Flambda; DAE.Gx, ...
	      DAE.Jlfv, Gk, Gl; Cinc'; Kjac];
      else
        if iterazione > 0
	  cont_eq = DAE.V(PQidx) - ...
		    sigma_corr*inc(PQidx+DAE.n+Bus.n) - ...
		    V_old(PQidx);
        else
	  cont_eq = lambda - sigma_corr*d_lambda - lambda_old;
	end
        AC = [DAE.Fx, DAE.Fy, Fk, Flambda; DAE.Gx, ...
	      DAE.Jlfv, Gk, Gl; Ljac; Kjac];
      end

      inc_corr = -AC\[DAE.f; DAE.gp; DAE.gq; cont_eq; 0];
      DAE.x = DAE.x + inc_corr(1:DAE.n);
      DAE.a = DAE.a + inc_corr(1+DAE.n: Bus.n+DAE.n);
      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));
      %disp([DAE.a;DAE.x(Syn.delta)])
      %pause
    end

    % Generator reactive power computations
    if qlim
      DAE.gq = 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);
      end
      Bus.Ql = DAE.gq;
      Q = DAE.glfq;
      Bus.Qg = Q + Bus.Ql;
      if PV.n
        Qmax_idx = find(Bus.Qg(PV.bus) > PV.con(:,6));
        Qmin_idx = find(Bus.Qg(PV.bus) < PV.con(:,7));
      end
      if SW.n
        Qswmax_idx = find(Bus.Qg(SW.bus) > SW.con(:,6));
        Qswmin_idx = find(Bus.Qg(SW.bus) < SW.con(:,7));
      end
    end
    if vlim
      Vmax_idx = find(DAE.V(PQ.bus) > PQ.con(:,6) & PQVmax);
      Vmin_idx = find(DAE.V(PQ.bus) < PQ.con(:,7) & PQVmin);
    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

    if ~isempty(Vmax_idx) & vlim
      Vmax_idx = Vmax_idx(1);