fm_gams.m

电力系统的psat

        Pdi = min(Pdi, Demand_old.con(:,5));
        Psi = min(Psi, Supply_old.con(:,4));

        % CPF step
        % ------------------------------------------------------
        if GAMS.basepl
          PQ.con(:,4) = PQ_old.con(:,4);
          PQ.con(:,5) = PQ_old.con(:,5);
          PV.con(:,4) = PV_old.con(:,4);
          SW.con(:,10) = SW_old.con(:,10);
          Snapshot(1).Pg = snappg;
        else
          PQ.con(:,4) = zeros(PQ_old.n,1);
          PQ.con(:,5) = zeros(PQ_old.n,1);
          PV.con(:,4) = zeros(PV_old.n,1);
          SW.con(:,10) = zeros(SW_old.n,1);
          Snapshot(1).Pg = zeros(Bus.n,1);
        end

        for kk = 1:Demand_old.n
          jj = find(PQ.bus == Demand_old.bus(kk));
          if ~isempty(jj)
            PQ.con(jj,4) = PQ.con(jj,4) + Pdi(kk);
            PQ.con(jj,5) = PQ.con(jj,5) + Pdi(kk).*tgd(kk);
          end
        end
        for kk = 1:Supply_old.n
          jj = find(PV.bus == Supply_old.bus(kk));
          if ~isempty(jj)
            PV.con(jj,4) = PV.con(jj,4) + Psi(kk);
          end
          jj = find(SW.bus == Supply_old.bus(kk));
          if ~isempty(jj)
            Snapshot(1).Pg(SW.bus(jj)) = ...
                Snapshot(1).Pg(SW.bus(jj)) + Psi(kk);
            %SW.con(jj,10) = SW.con(jj,10) + Psi(kk);
          end
        end

        DAE.V = Vci;
        DAE.a = aci;
        PV.con(:,5) = DAE.V(PV.bus);
        SW.con(:,4) = 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

        if isempty(CPF.lambda)
          fm_disp([' * CPF solution: <empty>'])
        elseif isnan(CPF.lambda)
          fm_disp([' * CPF solution: <NaN>'])
        else
          fm_disp([' * CPF solution: ',num2str(CPF.lambda-1)])
        end

        if isnan(CPF.lambda)
          stop_opf = 1;
          break
        end
        if isempty(CPF.lambda)
          stop_opf = 1;
          break
        end
        if CPF.lambda ~= lambdaci
          CPF.lambda = CPF.lambda - 0.995;
        end

        if CPF.lambda < lambdaci & abs(ml) <= 1e-5
          ml = 0;
          CPF.lambda = lmin+1e-5;
        end
        if CPF.lambda < lmin & abs(ml) > 1e-5
          fm_disp([' * Decrease Delta Ps and Delta Pd'])
          delta = 0.5*delta;
          if delta < 5e-8
            fm_disp([' * CPF method cannot find a higher lambda'])
            stop_opf = 1;
            break
          end
          repeat_cpf = 1;
        else
          repeat_cpf = 0;
        end

        % maximum lambda increment
        if (CPF.lambda - lmin) > 0.025 % & (abs(ml) > 1e-5 | CPF.lambda > 0.6)
          fm_disp(['lambda critical = ',num2str(CPF.lambda)])
          fm_disp(['Limit lambda increment to threshold (0.025)'])
          CPF.lambda = lmin + 0.025;
        end

        % stopping criterion
        % ------------------------------------------------------
        if last_point
          fm_disp('Reached maximum lambda.')
          if CPF.lambda > lmin
            fm_disp('Desired maximum lambda is not critical.')
          end
          stop_opf = 1;
          break
        end
        if i >= CPF.nump
          fm_disp('Reached maximum # of continuation steps.')
          stop_opf = 1;
          break
        end
        if CPF.lambda >= GAMS.lmax
          CPF.lambda = GAMS.lmax;
          last_point = 1;
        end
        if CPF.lambda == 0
          fm_disp('Base case solution is likely unfeasible.')
          stop_opf = 1;
          break
        end
        if abs(lmin-CPF.lambda) < 1e-5
          %fm_disp(['|lambda(i+1) - lambda(i)| = ', ...
          %         num2str(abs(lmin-CPF.lambda))])
          fm_disp('Lambda increment is lower than the desired tolerance.')
          stop_opf = 1;
          break
        elseif ~repeat_cpf
          if abs(ml) < 1e-5
            lmin = CPF.lambda+0.001;
            lmax = CPF.lambda+0.001;
          else
            lmin = CPF.lambda;
            lmax = CPF.lambda;
          end
          break
        end
      end
      %end
      if stop_opf, break, end
    end

    % restore original data and settings
    % --------------------------------------------------------
    DAE.V = Snapshot(1).V;
    DAE.a = Snapshot(1).ang;
    Snapshot(1).Pg = snappg;
    Bus.Pg = Bus_old.Pg;
    Bus.Qg = Bus_old.Qg;
    Bus.Pl = Bus_old.Pl;
    Bus.Ql = Bus_old.Ql;
    DAE.J11 = J11old;
    DAE.J12 = J12old;
    DAE.J21 = J21old;
    DAE.J22 = J22old;
    PV = PV_old;
    SW = SW_old;
    PQ = PQ_old;
    CPF = CPF_old;
    CPF.init = 4;
    Supply = Supply_old;
    Demand = Demand_old;

    Varout.t = [];
    Varout.x = [];
    Varout.V = [];
    Varout.ang = [];
    Varout.p = [];
    Varout.q = [];
    Varout.prflow = [];
    Varout.qrflow = [];
    Varout.psflow = [];
    Varout.qsflow = [];
    Varout.Pm = [];
    Varout.Vf = [];

    fm_disp
    % uncomment to plot [dP/d lambda] instead of [P]
    %Ps = DPs;
    %Pd = DPd;

  else

    fm_disp('Continuation OPF not implemented yet...')
    cd(currentpath)
    return

  end

end

% -------------------------------------------------------------------
% Output
% -------------------------------------------------------------------
MVA = Settings.mva;

if Settings.octave
  bslash = '';
else
  bslash = char(92);
end

if GAMS.method == 6, type = 3; end

if type == 2 | type == 3

  switch GAMS.flow
   case 0, flow = 'I_';
   case 1, flow = 'I_';
   case 2, flow = 'P_';
   case 3, flow = 'S_';
  end

  Lf = cellstr(num2str(Line.from));
  Lt = cellstr(num2str(Line.to));

  TD = MVA*sum(Pd')';
  if type == 2
    TD = MVA*sum(Pd,2);
    TTL = TD + MVA*sum(PQ.con(:,4))*Ch.val;
    TL = MVA*sum(Ps')' + MVA*sum(Bus.Pg)*Ch.val - TTL;
    TBL = TL - MVA*Snapshot(1).Ploss*Ch.val;
    for i = 1:size(Ps,1)
      PG(i,:) = full(sparse(1,iBPs,Ps(i,:),1,Bus.n)+Ch.val(i)*Bus.Pg')*MVA;
    end
    for i = 1:size(Pd,1)
      PL(i,:) = full(sparse(1,iBPd,Pd(i,:),1,Bus.n)+Ch.val(i)*Bus.Pl')*MVA;
    end
  elseif type == 3
    TTL = TD + MVA*sum(PQ.con(:,4));
    TL = MVA*sum(Ps')' + MVA*sum(Bus.Pg) - TTL;
    TBL = TL - MVA*Snapshot(1).Ploss;
    for i = 1:size(Ps,1)
      PG(i,:) = full(sparse(1,iBPs,Ps(i,:),1,Bus.n)+Bus.Pg')*MVA;
    end
    for i = 1:size(Pd,1)
      PL(i,:) = full(sparse(1,iBPd,Pd(i,:),1,Bus.n)+Bus.Pl')*MVA;
    end
  end
  PayS = -PG.*ro;
  PayD = PL.*ro;
  ISO = sum(PayS')'+sum(PayD')';
  if GAMS.method == 4 | GAMS.method == 6
    MLC = TTL.*(1+lambdac);
  elseif GAMS.method == 5
    MLC = TTL.*lambdac;
  end

  OPF.uname = strcat('PS_',{Varname.bus{Supply.bus}}');
  OPF.uname = [OPF.uname;strcat('PD_',{Varname.bus{Demand.bus}}')];
  OPF.uname = [OPF.uname;strcat('PG_',{Varname.bus{:}}')];
  OPF.uname = [OPF.uname;strcat('PL_',{Varname.bus{:}}')];
  OPF.uname = [OPF.uname;strcat('Pay_S_',{Varname.bus{:}}')];
  OPF.uname = [OPF.uname;strcat('Pay_D_',{Varname.bus{:}}')];
  OPF.uname = [OPF.uname;strcat('theta_',{Varname.bus{:}}')];
  OPF.uname = [OPF.uname;strcat('V_',{Varname.bus{:}}')];
  OPF.uname = [OPF.uname;strcat('Qg_',{Varname.bus{iBQg}}')];
  if GAMS.method > 2 & GAMS.method ~= 5
    OPF.uname = [OPF.uname;strcat('LMP_',{Varname.bus{:}}')];
    OPF.uname = [OPF.uname;strcat('NCP_',{Varname.bus{:}}')];
  elseif GAMS.method == 2
    OPF.uname = [OPF.uname;strcat('LMP_',{Varname.bus{:}}')];
  elseif GAMS.method == 5
    OPF.uname = [OPF.uname;strcat(bslash,'rho_',{Varname.bus{:}}')];
  else
    OPF.uname = [OPF.uname;{'MCP'}];
  end
  OPF.uname = [OPF.uname;strcat(flow,Lf,'-',Lt)];
  OPF.uname = [OPF.uname;strcat(flow,Lt,'-',Lf)];
  OPF.uname = [OPF.uname;{'Total Demand';'TTL';'Total Losses'; ...
                      'Total Bid Losses';'IMO Pay'}];
  if GAMS.method >= 4
    OPF.uname = [OPF.uname;{'MLC'}];
    OPF.uname = [OPF.uname;{'ALC'}];
  end

  OPF.fname = strcat('P_{S',{Varname.bus{Supply.bus}}','}');
  OPF.fname = [OPF.fname;strcat('P_{D',{Varname.bus{Demand.bus}}','}')];
  OPF.fname = [OPF.fname;strcat('P_{G',{Varname.bus{:}}','}')];
  OPF.fname = [OPF.fname;strcat('P_{L',{Varname.bus{:}}','}')];
  OPF.fname = [OPF.fname;strcat('Pay_{S',{Varname.bus{:}}','}')];
  OPF.fname = [OPF.fname;strcat('Pay_{D',{Varname.bus{:}}','}')];
  OPF.fname = [OPF.fname;strcat(bslash,'theta_{',{Varname.bus{:}}','}')];
  OPF.fname = [OPF.fname;strcat('V_{',{Varname.bus{:}}','}')];
  OPF.fname = [OPF.fname;strcat('Q_{g',{Varname.bus{iBQg}}','}')];
  if GAMS.method > 2  & GAMS.method ~= 5
    OPF.fname = [OPF.fname;strcat('LMP_{',{Varname.bus{:}}','}')];
    OPF.fname = [OPF.fname;strcat('NCP_{',{Varname.bus{:}}','}')];
  elseif GAMS.method == 2
    OPF.fname = [OPF.fname;strcat('LMP_{',{Varname.bus{:}}','}')];
  elseif GAMS.method == 5
    OPF.fname = [OPF.fname;strcat(bslash,'rho_',{Varname.bus{:}}')];
  else
    OPF.fname = [OPF.fname;{'MCP'}];
  end
  OPF.fname = [OPF.fname;strcat(flow,'{',Lf,'-',Lt,'}')];
  OPF.fname = [OPF.fname;strcat(flow,'{',Lt,'-',Lf,'}')];
  OPF.fname = [OPF.fname;{'Total Demand';'TTL';'Total Losses'; ...
                   'Total Bid Losses';'IMO Pay'}];
  if GAMS.method >= 4
    OPF.fname = [OPF.fname;{'MLC'}];
    OPF.fname = [OPF.fname;{'ALC'}];
  end

  switch GAMS.method
   case 3 % OPF

    OPF.varout = [Ps*MVA,Pd*MVA,PG,PL,PayS,PayD,a,V,Qg(:,iBQg)*MVA, ...
                  ro,NCP,Pij,Pji,TD,TTL,TL,TBL,ISO];

   case {4,6} % VSC-OPF

    OPF.uname = [OPF.uname;{'lambda_c';'kg_c'}];
    OPF.uname = [OPF.uname;strcat('thetac_',{Varname.bus{:}}')];
    OPF.uname = [OPF.uname;strcat('Vc_',{Varname.bus{:}}')];
    OPF.uname = [OPF.uname;strcat('Qgc_',{Varname.bus{iBQg}}')];
    OPF.uname = [OPF.uname;strcat(flow,'c',Lf,'-',Lt)];
    OPF.uname = [OPF.uname;strcat(flow,'c',Lt,'-',Lf)];

    OPF.fname = [OPF.fname;{[bslash,'lambda_c'];'k_g_c'}];
    OPF.fname = [OPF.fname;strcat(bslash,'theta_{c',{Varname.bus{:}}','}')];
    OPF.fname = [OPF.fname;strcat('V_{c',{Varname.bus{:}}','}')];
    OPF.fname = [OPF.fname;strcat('Q_{gc',{Varname.bus{iBQg}}','}')];
    OPF.fname = [OPF.fname;strcat(flow,'{c',Lf,'-',Lt,'}')];
    OPF.fname = [OPF.fname;strcat(flow,'{c',Lt,'-',Lf,'}')];

    OPF.varout = [Ps*MVA,Pd*MVA,PG,PL,PayS,PayD,a,V,Qg(:,iBQg)*MVA, ...
                  ro,NCP,Pij,Pji,TD,TTL,TL,TBL,ISO,MLC,MLC-TTL,lambdac,kg, ...
                  ac,Vc,Qgc(:,iBQg)*MVA,Pijc,Pjic];

   case 5 % MLC

    OPF.uname = [OPF.uname;{'lambda_c';'kg_c'}];
    OPF.fname = [OPF.fname;{bslash,'lambda_c';'k_g_c'}];
    OPF.varout = [Ps*MVA,Pd*MVA,PG,PL,PayS,PayD,a,V,Qg(:,iBQg)*MVA, ...
                  ro,Pij,Pji,TD,TTL,TL,TBL,ISO,MLC,lambdac,kg];

   otherwise % SA and MCM
    OPF.varout = [Ps*MVA,Pd*MVA,PG,PL,PayS,PayD,a,V,Qg(:,iBQg)*MVA, ...
                  MCP,Pij,Pji,TD,TTL,TL,TBL,ISO];

  end

  if GAMS.method == 6 % Continuation OPF
    Settings.xlabel = [bslash,'lambda (loading parameter)'];
    GAMS.hours = Lambda';
  elseif type == 2 % Multi Period Auction
    OPF.varout = OPF.varout([2:end],:);
    Settings.xlabel = 'hour [h]';
    GAMS.hours = [1:Ypdp.len]';
  elseif type == 3 % Pareto Set Single Period Auction
    Settings.xlabel = [bslash,'omega (weighting factor)'];
    GAMS.hours = GAMS.omega';
  end

  fm_disp(' ---------------------------------------------------------------')
  fm_disp([' Check file ',Path.psat,'fm_gams.lst for GAMS report.'])
  if strcmp(control,'6')
    fm_disp([' PSAT-GAMS Optimization Routine completed in ', ...
             num2str(etime(clock,initial_time)),' s'])
  else
    fm_disp([' PSAT-GAMS Optimization Routine completed in ',num2str(toc),' s'])
  end

  if noDem
    Demand.con = [];
    Demand.n = 0;
    Demand.bus = [];
  end

  if ~GAMS.basepl
    Bus.Pl = buspl;
    Bus.Ql = busql;
    PQ.con(:,[4 5]) = pqcon;
  end
  if ~GAMS.basepg
    Snapshot(1).Ploss = ploss;
    Bus.Pg = buspg;
    Bus.Qg = busqg;
    PV.con(:,4) = pvcon;
  end

  % restore original bus power injections
  Bus.Pg = Snapshot(1).Pg;
  Bus.Qg = Snapshot(1).Qg;
  Bus.Pl = Snapshot(1).Pl;
  Bus.Ql = Snapshot(1).Ql;

  return

end

if type == 4

  Ps = Ps(2,:)';
  Pd = Pd(2,:)';
  V = V(2,:)';
  a = a(2,:)';
  Qg = Qg(2,iBQg)';
  Pij = Pij(2,:)';
  Pji = Pji(2,:)';
  if GAMS.method <= 2
    MCP = MCP(2,:);
  end
  if GAMS.method >= 3
    ro = ro(2,:)';
    if GAMS.method ~= 5
      NCP = NCP(2,:)';
    end
  end
  if GAMS.method == 4 | GAMS.method == 6
    Vc = Vc(2,:)';
    ac = ac(2,:)';
    Qgc = Qgc(2,iBQg)';
    Pijc = Pijc(2,:)';
    Pjic = Pjic(2,:)';
  end
  if GAMS.method >= 4
    lambdac = lambdac(2);
    kg = kg(2);
  end
end

Demand.con(:,7) = Pd;
Supply.con(:,6) = Ps;

if GAMS.method == 4 | GAMS.method == 6
  DAE.V = Vc;
  DAE.a = ac;