fm_gams.m

这是一个很适合研究和学习用的电力系统仿真软件

        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 = restore(PV);
    SW = restore(SW);
    PQ = pqreset(PQ,'all');

    CPF = CPF_old;
    CPF.init = 4;

    Varout.t = [];
    Varout.vars = [];

    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;
TPQ = MVA*totp(PQ);

% character for backslash
bslash = char(92);

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 + TPQ*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 + TPQ;
    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

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

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

  switch GAMS.method
   case 3 % OPF

    Varout.vars = [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

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

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

    Varout.vars = [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

    Varname.uvars = [Varname.uvars;{'lambda_c';'kg_c'}];
    Varname.fvars = [Varname.fvars;{bslash,'lambda_c';'k_g_c'}];
    Varout.vars = [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
    Varout.vars = [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)'];
    Varout.t = Lambda';
  elseif type == 2 % Multi Period Auction
    Varout.vars = Varout.vars([2:end],:);
    Settings.xlabel = 'hour [h]';
    Varout.t = [1:Ypdp.len]';
  elseif type == 3 % Pareto Set Single Period Auction
    Settings.xlabel = [bslash,'omega (weighting factor)'];
    Varout.t = GAMS.omega';
  end
  
  Varout.idx = [1:length(Varout.vars(1,:))];

  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

  Demand = restore(Demand);

  if ~GAMS.basepl
    Bus.Pl = buspl;
    Bus.Ql = busql;
    PQ = pqreset(PQ,'all');
  end
  if ~GAMS.basepg
    Snapshot(1).Ploss = ploss;
    Bus.Pg = buspg;
    Bus.Qg = busqg;
    PV = pvreset(PV,'all');
  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 = pset(Demand,Pd);
Supply = pset(Supply,Ps);

if GAMS.method == 4 | GAMS.method == 6
  DAE.V = Vc;
  DAE.a = ac;
  fm_lf(1)
  glfpc = DAE.glfp;
  glfqc = DAE.glfq;
end

if GAMS.method >= 3
  DAE.V = V;
  DAE.a = a;
  fm_lf(1)
  Qg = Qg(iBQg);
end

if GAMS.method == 1
  ro = MCP*ones(Bus.n,1);
end
if GAMS.method == 2
  [rows,cols] = size(MCP);
  if rows == 1,
    ro = MCP';
  else
    ro = MCP;
  end
end

Qgmin = Qgmin(iBQg);
Qgmax = Qgmax(iBQg);
if GAMS.basepl
  PG = full((sparse(iBPs,1,Ps,Bus.n,1)+Bus.Pg)*MVA);
  PL = full((sparse(iBPd,1,Pd,Bus.n,1)+Bus.Pl)*MVA);
else
  PG = full(sparse(iBPs,1,Ps,Bus.n,1)*MVA);