fm_gams.m

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

    for i = 1:numh
      [Piji,Pjii,Qgi] = updatePF(Pd(i,:)',Ps(i,:)',iBQg);
      a(i,:) = DAE.a';
      V(i,:) = DAE.V';
      Pij(i,:) = Piji';
      Pji(i,:) = Pjii';
      Qg(i,iBQg) = Qgi';
    end
    ro = MCP;

  end

 % ------------------------------------------------------------------
 case '3' % S T A N D A R D   O P T I M A L   P O W E R   F L O W
 % ------------------------------------------------------------------

  if type == 1 % Single Period Auction

    [Ps,Pd,V,a,Qg,ro,Pij,Pji,mV,mFij,mFji,modelstat,solvestat] = ...
        psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV,nBusref,control, ...
             flow,Gh,Bh,Li,Lj,Ps_idx,Pd_idx,S,D,X,L);
    NCP = compNCP(V,a,mV,mFij,mFji);

  elseif ~rem(type,2) % Single/Multi Period Auction with UC

    [Ps,Pd,V,a,Qg,ro,Pij,Pji,mV,mFij,mFji,modelstat,solvestat] = ...
        psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV,nBusref,nH,control, ...
             flow,Gh,Bh,Li,Lj,Ps_idx,Pd_idx,S,D,X,L,Ch);

    NCP = zeros(length(Ps(:,1)),Bus.n);
    for i = 1:length(Ps(:,1))
      NCPi = compNCP(V(i,:)',a(i,:)',mV(i,:)',mFij(i,:)',mFji(i,:)');
      NCP(i,:) = NCPi';
    end

  end

 % ------------------------------------------------------------------
 case '7' % C O N G E S T I O N   M A N A G E M E N T
 % ------------------------------------------------------------------

  %lambda_values = [0.0:0.01:0.61];
  %n_lambda = length(lambda_values);
  %GAMS.dpgup = zeros(Supply.n,n_lambda);
  %GAMS.dpgdw = zeros(Supply.n,n_lambda);
  %GAMS.dpdup = zeros(Demand.n,n_lambda);
  %GAMS.dpddw = zeros(Demand.n,n_lambda);
  %for i = 1:length(lambda_values)
  %lambda.val(1) = lambda_values(i);

  iteration = 0;

  idx_gen = zeros(Supply.n,1);
  Psc_idx = Ps_idx;
  Psm_idx = zeros(Bus.n,Supply.n);
 
  while 1
    [Ps,Pd,dPSup,dPSdw,dPDup,dPDdw,V,a,Qg,ro,Pij,Pji,mV,mFij,mFji, ...
     lambdac,kg,Vc,ac,Qgc,Pijc,Pjic,Pceq,lambdam,modelstat,solvestat] = ...
        psatgams('fm_cong',nBus,nLine,nPs,nPd,nBusref,nSW,nPV,control,flow, ...
                 Gh,Bh,Ghc,Bhc,Li,Lj,Ps_idx,Psc_idx,Psm_idx,Pd_idx, ...
                 SW_idx,PV_idx,S,D,X,L,lambda);

    iteration = iteration + 1;
    if iteration > 10
      fm_disp('* * * Maximum number of iteration with no convergence!')
      break
    end
    idx = psupper(Supply,(1+lambdac+kg)*Ps);
    if sum(idx_gen(idx)) == length(idx)
      fm_disp(['* * * iter = ',num2str(iteration), ...
               ', #viol., ',num2str(length(find(idx_gen))), ...
               ' lambda = ', num2str(lambdac), ...
               ' kg = ', num2str(kg)])
      break
    else
      % loop until there are no violations of power supply limits
      idx_gen(idx) = 1;
      fm_disp(['* * * iter = ',num2str(iteration),', #viol. = ', ...
               num2str(length(idx)),', lambda = ', ...
               num2str(lambdac),' kg = ', num2str(kg)])
      drawnow;
      Psc_idx = psidx(Supply,~idx_gen);
      Psm_idx = psidx(Supply,idx_gen);
      
    end
  end

  %GAMS.dpgup(:,i) = dPSup;
  %GAMS.dpgdw(:,i) = dPSdw;
  %GAMS.dpdup(:,i) = dPDup;
  %GAMS.dpddw(:,i) = dPDdw;
  %if ~rem(i,10), disp(['Current lambda = ',num2str(lambda_values(i))]),end
  %end
  %GAMS.lvals = lambda_values;
  NCP = compNCP(V,a,mV,mFij,mFji);


 % ------------------------------------------------------------------
 case '4' % V O L T A G E   S T A B I L I T Y   C O N S T R A I N E D
          % O P T I M A L   P O W E R    F L O W
 % ------------------------------------------------------------------

  if type == 1 % Single Period Auction
    [Ps,Pd,V,a,Qg,ro,Pij,Pji,mV,mFij,mFji, ...
     lambdac,kg,Vc,ac,Qgc,Pijc,Pjic,Pceq,modelstat,solvestat] = ...
        psatgams(file,nBus,nLine,nPs,nPd,nBusref,nSW,nPV,control,flow, ...
                 Gh,Bh,Ghc,Bhc,Li,Lj,Ps_idx,Pd_idx,SW_idx,PV_idx,S,D,X,L,lambda);
    NCP = compNCP(V,a,mV,mFij,mFji);

  elseif ~rem(type,2) % Single/Multi Period Auction with UC

    [Ps,Pd,V,a,Qg,ro,Pij,Pji,mV,mFij,mFji, ...
     lambdac,kg,Vc,ac,Qgc,Pijc,Pjic,Pceq,modelstat,solvestat] = ...
        psatgams(file,nBus,nLine,nPs,nPd,nBusref,nH,control,flow, ...
                 Gh,Bh,Ghc,Bhc,Li,Lj,Ps_idx,Pd_idx,S,D,X,L,lambda,Ch);

    NCP = zeros(length(lambdac),Bus.n);
    for i = 1:length(lambdac)
      NCPi = compNCP(V(i,:)',a(i,:)',mV(i,:)',mFij(i,:)',mFji(i,:)');
      NCP(i,:) = NCPi';
    end

  elseif type == 3 % Pareto Set Single Period Auction

    fm_disp
    for i = 1:length(omega)
      fm_disp(sprintf(' VSC-OPF #%d, %3.1f%% - omega: %5.4f', ...
		      i,100*i/length(omega),omega(i)))
      lambda.val = [GAMS.lmin(1),GAMS.lmax(1),omega(i),GAMS.line];

      [Psi,Pdi,Vi,ai,Qgi,roi,Piji,Pjii,mV,mFij,mFji, ...
       lambdaci,kgi,Vci,aci,Qgci,Pijci,Pjici,Pceq,modelstat,solvestat] = ...
          psatgams(file,nBus,nLine,nPs,nPd,nBusref,nSW,nPV,control,flow, ...
                   Gh,Bh,Ghc,Bhc,Li,Lj,Ps_idx,Pd_idx,SW_idx,PV_idx, ...
                   S,D,X,L,lambda);
      gams_mstat(modelstat)
      gams_sstat(solvestat)

      Ps(i,:) = Psi';
      Pd(i,:) = Pdi';
      V(i,:) = Vi';
      a(i,:) = ai';
      Qg(i,:) = Qgi';
      ro(i,:) = roi';
      Pij(i,:) = Piji';
      Pji(i,:) = Pjii';
      lambdac(i,:) = lambdaci';
      kg(i,:) = kgi';
      Vc(i,:) = Vci';
      ac(i,:) = aci';
      Qgc(i,:) = Qgci';
      Pijc(i,:) = Pijci';
      Pjic(i,:) = Pjici';
      NCPi = compNCP(Vi,ai,mV,mFij,mFji);
      NCP(i,:) = NCPi';

    end
    fm_disp
  end

 % ------------------------------------------------------------------
 case '5' % M A X I M U M   L O A D I N G   C O N D I T I O N
 % ------------------------------------------------------------------

  if type == 1 % Single Period Auction

    [Ps,Pd,V,a,Qg,ro,Pij,Pji,lambdac,kg,modelstat,solvestat] = ...
        psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV,nBusref, ...
		 control,flow,Gh,Bh,Li,Lj,Ps_idx,Pd_idx, ...
                 SW_idx,PV_idx,S,D,X,L);

  elseif  ~rem(type,2) % Single/Multi Period Auction with UC

    [Ps,Pd,V,a,Qg,ro,Pij,Pji,lambdac,kg,modelstat,solvestat] = ...
        psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV,nBusref,nH, ...
		 control,flow,Gh,Bh,Li,Lj,Ps_idx,Pd_idx,SW_idx, ...
                 PV_idx,S,D,X,L,Ch);

  end

 % ------------------------------------------------------------------
 case '6' % C O N T I N U A T I O N
          % O P T I M A L   P O W E R    F L O W
 % ------------------------------------------------------------------

  initial_time = clock;
  if type == 1 % single period OPF, no discrete variables

    % number of steps
    i = 0;
    last_point = 0;
    % initial lambda = 0. Base case has to be feasible
    lmin = 0;
    lmax = 0;
    Lambda = lmin;
    % save actual CPF settings
    CPF_old = CPF;
    %CPF.nump = 50;
    CPF.show = 0;
    CPF.type = 3;
    CPF.sbus = 0;
    CPF.vlim = 1;
    CPF.ilim = 1;
    CPF.qlim = 1;
    CPF.init = 0;
    CPF.step = 0.25;
    control = '6';

    % save actual power flow data
    % ------------------------------------------------------
    snappg = Snapshot(1).Pg;
    %Snapshot(1).Pg = [];
    J11old = DAE.J11;
    J12old = DAE.J12;
    J21old = DAE.J21;
    J22old = DAE.J22;
    Bus_old = Bus;

    % defining voltage limits
    [Vmax,Vmin] = fm_vlim(1.2,0.8);

    fm_disp
    stop_opf = 0;

    while 1

      % OPF step
      % ------------------------------------------------------
      i = i + 1;
      fm_disp(sprintf('Continuation OPF #%d, lambda_c = %5.4f', ...
		      i,lmin))
      lambda.val = [lmin,lmax,0,GAMS.line];

      % call GAMS
      [Psi,Pdi,Vi,ai,Qgi,roi,Piji,Pjii,mV,mFij,mFji, ...
       lambdaci,kgi,Vci,aci,Qgci,Pijci,Pjici,mPceq,ml,modelstat,solvestat] = ...
        psatgams(file,nBus,nLine,nPs,nPd,nBusref,nSW,nPV,control,flow, ...
             Gh,Bh,Ghc,Bhc,Li,Lj,Ps_idx,Pd_idx,SW_idx,PV_idx,S,D,X,L,lambda);
      gams_mstat(modelstat)
      gams_sstat(solvestat)

      Lambda(i,1) = lambdaci;
      Ps(i,:) = Psi';
      Pd(i,:) = Pdi';
      V(i,:) = Vi';
      a(i,:) = ai';
      Qg(i,:) = Qgi';
      ro(i,:) = roi';
      Pij(i,:) = Piji';
      Pji(i,:) = Pjii';
      lambdac(i,:) = lambdaci;
      kg(i,:) = kgi;
      Vc(i,:) = Vci';
      ac(i,:) = aci';
      Qgc(i,:) = Qgci';
      Pijc(i,:) = Pijci';
      Pjic(i,:) = Pjici';
      NCPi = compNCP(Vi,ai,mV,mFij,mFji);
      NCP(i,:) = NCPi';
      ML(i,1) = ml;

      % check consistency of the solution (LMP > 0)
      if modelstat > 3 %min(abs(roi)) < 1e-5

        fm_disp('Unfeasible OPF solution. Discarding last solution.')
        Lambda(end) = [];
        Ps(end,:) = [];
        Pd(end,:) = [];
        V(end,:) = [];
        a(end,:) = [];
        Qg(end,:) = [];
        ro(end,:) = [];
        Pij(end,:) = [];
        Pji(end,:) = [];
        lambdac(end) = [];
        kg(end) = [];
        Vc(end,:) = [];
        ac(end,:) = [];
        Qgc(end,:) = [];
        Pijc(end,:) = [];
        Pjic(end,:) = [];
        NCP(end,:) = [];
        ML(end) = [];
        lambdaci = lambdac(end,:);
        break

      end

      % ------------------------------------------------------
      % Bid variations to allow loading parameter increase
      %
      %                      d mu_Pceq_i
      % D P_i = -sign(P_i)  -------------
      %                      d mu_lambda
      %
      % where:
      %
      % P_i = power bid i
      % mu_Pceq_i = Lagrangian multiplier of critical PF eq. i
      % mu_lambda = Lagrangian multiplier of lambda
      % ------------------------------------------------------

      delta = 0.05;

      while 1

        if abs(ml) > 1e-5
          deltaPd =  ml./mPceq(Demand.bus)/(1+lambdaci);
          deltaPs = -ml./mPceq(Supply.bus)/(1+lambdaci+kgi);
          delta_max = norm([deltaPs; deltaPd]);
          if delta_max == 0, delta_max = 1; end
          deltaPd = deltaPd/delta_max;
          deltaPs = deltaPs/delta_max;
        else
          deltaPd = zeros(Demand.n,1);
          deltaPs = zeros(Supply.n,1);
        end
        %ml
        %mPceq
        %delta_max = max(norm([deltaPs; deltaPd]));
        %if delta_max == 0, delta_max = 1; end
        DPs(i,:) = deltaPs'/Settings.mva;
        DPd(i,:) = deltaPd'/Settings.mva;

        Pdi = pdbound(Demand,Pd(i,:)' + delta*deltaPd.*Pd(i,:)');
        Psi = psbound(Supply,Ps(i,:)' + delta*deltaPs.*Ps(i,:)');

        % CPF step
        % ------------------------------------------------------
        if GAMS.basepl
          PQ = pqreset(PQ,'all');
          PV = pvreset(PV,'all');
          SW = swreset(SW,'all');
          Snapshot(1).Pg = snappg;
        else
          PQ = pqzero(PQ,'all');
          PV = pvzero(PV,'all');
          SW = swzero(SW,'all');
          Snapshot(1).Pg = zeros(Bus.n,1);
        end

        Demand = pset(Demand,Pdi);
        pqsum(Demand,1);
        Supply = pset(Supply,Psi);
        pgsum(Supply,1);
        swsum(Supply,1);

        DAE.V = Vci;
        DAE.a = aci;
        PV = setvg(PV,'all',DAE.V(PV.bus));
        SW = setvg(SW,'all',DAE.V(SW.bus));
        DAE.x = Snapshot(1).x;
        Bus.Pg = Bus_old.Pg;
        Bus.Qg = Bus_old.Qg;
        Bus.Pl = Bus_old.Pl;
        Bus.Ql = Bus_old.Ql;

        % avoid aborting CPF routine due to limits
        % ------------------------------------------------------

        % voltage limits
        idx = find(abs(DAE.V-Vmax) < CPF.tolv | DAE.V > Vmax);
        if ~isempty(idx)
          DAE.V(idx) = Vmax(idx)-1e-6-CPF.tolv;
        end
        idx = find(abs(DAE.V-Vmin) < CPF.tolv | DAE.V < Vmin);
        if ~isempty(idx)
          DAE.V(idx) = Vmin(idx)+1e-6+CPF.tolv;
        end

        CPF.kg = 0;
        CPF.lambda = 1; %lambdaci + 1;
        CPF.linit = 1+lambdaci*0.25;
        CPF.init = 0;

        % set contingency for CPF analysis
        if GAMS.line
          Line.con = line_cont;
          fm_y;
        end

        % ---------------------------------------------
        % call continuation power flow routine
        fm_cpf('gams');
        %CPF.lambda = CPF.lambda + 1;
        % ---------------------------------------------

        % reset admittance line
        if GAMS.line
          Line.Y = Y_orig;
          Line.con = line_orig;
        end