fm_gams.m

电力系统的psat

end

X.labels = {cellstr(num2str([1:Bus.n]')), ...
	    {'V0','t0','Pg0','Qg0','Pl0','Ql0', ...
	     'Qgmax','Qgmin','Vmax','Vmin','ksw','kpv'}};
L.labels = {cellstr(num2str([1:Line.n]')), ...
	    {'g','b','g0','b0','Pijmax','Pjimax'}};
lambda.labels = {'lmin';'lmax';'omega';'line'};

S.name = 'S';
D.name = 'D';
X.name = 'X';
L.name = 'N';
lambda.name = 'lambda';

if type == 2
  Ch.val = zeros(Ypdp.len+1,1);
  Ch.val = [0;Ypdp.day(:,Ypdp.d)]*Ypdp.week(Ypdp.w)* ...
      Ypdp.year(Ypdp.y)/100/100/100;
  Ch.val(1) = Ch.val(2);
  for idx = 1:Ypdp.len+1
    Ch.labels{1,idx} = ['H',num2str(idx-1)];
  end
  %Ch.labels = strrep(strcat('H',cellstr(num2str([0:24]'))),' ','');
  Ch.name = 'Ch';
elseif type == 4
  Ch.val = [1;1];
  Ch.name = 'Ch';
  Ch.labels = {'H0';'H1'};
end

% ------------------------------------------------------------
% Launch GAMS solver
% ------------------------------------------------------------

control = int2str(method);
flow = int2str(GAMS.flow);
currentpath = pwd;
file = 'fm_gams';
if ~rem(type,2)
  file = [file,'2'];
end
if GAMS.libinclude
  file = [file,' ',GAMS.ldir];
end

%if clpsat.init
%  cd(Path.psat)
%end

switch control

 % ------------------------------------------------------------------
 case '1' % S I M P L E   A U C T I O N
 % ------------------------------------------------------------------

  if type == 1 % Single Period Auction
    [Ps,Pd,MCP,modelstat,solvestat] = psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV, ...
                       nBusref,control,flow,S,D,X);
    [Pij,Pji,Qg] = updatePF(Pd,Ps,iBQg);

  elseif ~rem(type,2) % Single/Multi Period Auction with UC
    [Ps,Pd,MCP,modelstat,solvestat] = psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV, ...
                       nBusref,nH,control,flow,S,D,X,Ch);
    numh = size(MCP,1);
    a = zeros(numh,Bus.n);
    V = zeros(numh,Bus.n);
    Qg = zeros(numh,Bus.n);
    Pij = zeros(numh,Line.n);
    Qij = zeros(numh,Line.n);
    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*ones(1,Bus.n);

  end

 % ------------------------------------------------------------------
 case '2' % M A R K E T   C L E A R I N G   M E C H A N I S M
 % ------------------------------------------------------------------

  if type == 1 % Single Period Auction
    [Ps,Pd,MCP,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);
    [Pij,Pji,Qg] = updatePF(Pd,Ps,iBQg);

  elseif ~rem(type,2) % Single/Multi Period Auction with UC
    [Ps,Pd,MCP,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);

    numh = size(MCP,1);
    a = zeros(numh,Bus.n);
    V = zeros(numh,Bus.n);
    Qg = zeros(numh,Bus.n);
    Pij = zeros(numh,Line.n);
    Qij = zeros(numh,Line.n);
    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);

  Psc_idx = Ps_idx;
  Psm_idx = zeros(Bus.n,Supply.n);
  iteration = 0;
  idx_gen = zeros(Supply.n,1);

  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);
    Psc_idx = Ps_idx;
    Psm_idx = zeros(Bus.n,Supply.n);
    iteration = iteration + 1;
    if iteration > 10
      fm_disp('* * * Maximum number of iteration with no convergence!')
      break
    end
    idx = find((1+lambdac+kg)*Ps > Psmax);
    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;
      idx = find(idx_gen);
      for iii = 1:length(idx),
        k = idx(iii);
        Psc_idx(Supply.bus(k),k) = 0;
        Psm_idx(Supply.bus(k),k) = 1;
      end
    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;
    PV_old = PV;
    SW_old = SW;
    PQ_old = PQ;
    Supply_old = Supply;
    Demand_old = Demand;
    tgd = Demand.con(:,4)./Demand.con(:,3);
    Demand.con = [];
    Demand.bus = [];
    Demand.n = 0;
    Supply.con = [];
    Supply.bus = [];
    Supply.n = 0;
    Bus_old = Bus;

    % defining voltage limits
    % --------------------------------------------------------
    Vmin = zeros(Bus.n,1);
    Vmax = zeros(Bus.n,1);
    Vmin(PQ.bus) = PQ.con(:,7);
    Vmax(PQ.bus) = PQ.con(:,6);
    Vmin(PV.bus) = PV.con(:,9);
    Vmax(PV.bus) = PV.con(:,8);
    Vmin(SW.bus) = SW.con(:,9);
    Vmax(SW.bus) = SW.con(:,8);
    Vmin(find(Vmin == 0)) = 0.8;
    Vmax(find(Vmax == 0)) = 1.2;

    fm_disp
    stop_opf = 0;
    %unfeas = 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_old.bus)/(1+lambdaci);
          deltaPs = -ml./mPceq(Supply_old.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_old.n,1);
          deltaPs = zeros(Supply_old.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 = Pd(i,:)' + delta*deltaPd.*Pd(i,:)';
        Psi = Ps(i,:)' + delta*deltaPs.*Ps(i,:)';
        Pdi = max(Pdi, Demand_old.con(:,6));
        Psi = max(Psi, Supply_old.con(:,5));