fm_gams.m

电力系统分析计算程序

 % ------------------------------------------------------------------

  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 = [];
    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 = getzeros(Bus);
        end

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

        DAE.y(Bus.a) = aci;
        DAE.y(Bus.v) = Vci;
        PV = setvg(PV,'all',DAE.y(PV.vbus));
        SW = setvg(SW,'all',DAE.y(SW.vbus));
        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
        Vbus = DAE.y(Bus.v);
        idx = find(abs(Vbus-Vmax) < CPF.tolv | Vbus > Vmax);
        if ~isempty(idx)
          DAE.y(idx+Bus.n) = Vmax(idx)-1e-6-CPF.tolv;
        end
        idx = find(abs(Vbus-Vmin) < CPF.tolv | Vbus < Vmin);
        if ~isempty(idx)
          DAE.y(idx+Bus.n) = 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
          status = Line.u(GAMS.line);
          Line = setstatus(Line,GAMS.line,0);
        end

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

        % reset admittance line
        if GAMS.line
          Line = setstatus(Line,GAMS.line,status);
        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.y = Snapshot(1).y;
    Snapshot(1).Pg = snappg;
    Bus.Pg = Bus_old.Pg;
    Bus.Qg = Bus_old.Qg;
    Bus.Pl = Bus_old.Pl;
    Bus.Ql = Bus_old.Ql;
    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.fr));
  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 = fm_strjoin('PS_',{Bus.names{Supply.bus}}');
  Varname.uvars = [Varname.uvars;fm_strjoin('PD_',{Bus.names{Demand.bus}}')];
  Varname.uvars = [Varname.uvars;fm_strjoin('PG_',{Bus.names{:}}')];
  Varname.uvars = [Varname.uvars;fm_strjoin('PL_',{Bus.names{:}}')];
  Varname.uvars = [Varname.uvars;fm_strjoin('Pay_S_',{Bus.names{:}}')];
  Varname.uvars = [Varname.uvars;fm_strjoin('Pay_D_',{Bus.names{:}}')];
  Varname.uvars = [Varname.uvars;fm_strjoin('theta_',{Bus.names{:}}')];
  Varname.uvars = [Varname.uvars;fm_strjoin('V_',{Bus.names{:}}')];
  Varname.uvars = [Varname.uvars;fm_strjoin('Qg_',{Bus.names{iBQg}}')];
  if GAMS.method > 2 & GAMS.method ~= 5
    Varname.uvars = [Varname.uvars;fm_strjoin('LMP_',{Bus.names{:}}')];
    Varname.uvars = [Varname.uvars;fm_strjoin('NCP_',{Bus.names{:}}')];
  elseif GAMS.method == 2
    Varname.uvars = [Varname.uvars;fm_strjoin('LMP_',{Bus.names{:}}')];
  elseif GAMS.method == 5
    Varname.uvars = [Varname.uvars;fm_strjoin(bslash,'rho_',{Bus.names{:}}')];
  else
    Varname.uvars = [Varname.uvars;{'MCP'}];
  end
  Varname.uvars = [Varname.uvars;fm_strjoin(flow,Lf,'-',Lt)];
  Varname.uvars = [Varname.uvars;fm_strjoin(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 = fm_strjoin('P_{S',{Bus.names{Supply.bus}}','}');
  Varname.fvars = [Varname.fvars;fm_strjoin('P_{D',{Bus.names{Demand.bus}}','}')];
  Varname.fvars = [Varname.fvars;fm_strjoin('P_{G',{Bus.names{:}}','}')];
  Varname.fvars = [Varname.fvars;fm_strjoin('P_{L',{Bus.names{:}}','}')];
  Varname.fvars = [Varname.fvars;fm_strjoin('Pay_{S',{Bus.names{:}}','}')];
  Varname.fvars = [Varname.fvars;fm_strjoin('Pay_{D',{Bus.names{:}}','}')];
  Varname.fvars = [Varname.fvars;fm_strjoin(bslash,'theta_{',{Bus.names{:}}','}')];
  Varname.fvars = [Varname.fvars;fm_strjoin('V_{',{Bus.names{:}}','}')];
  Varname.fvars = [Varname.fvars;fm_strjoin('Q_{g',{Bus.names{iBQg}}','}')];
  if GAMS.method > 2  & GAMS.method ~= 5
    Varname.fvars = [Varname.fvars;fm_strjoin('LMP_{',{Bus.names{:}}','}')];
    Varname.fvars = [Varname.fvars;fm_strjoin('NCP_{',{Bus.names{:}}','}')];
  elseif GAMS.method == 2
    Varname.fvars = [Varname.fvars;fm_strjoin('LMP_{',{Bus.names{:}}','}')];
  elseif GAMS.method == 5
    Varname.fvars = [Varname.fvars;fm_strjoin(bslash,'rho_',{Bus.names{:}}')];
  else
    Varname.fvars = [Varname.fvars;{'MCP'}];
  end
  Varname.fvars = [Varname.fvars;fm_strjoin(flow,'{',Lf,'-',Lt,'}')];
  Varname.fvars = [Varname.fvars;fm_strjoin(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;fm_strjoin('thetac_',{Bus.names{:}}')];
    Varname.uvars = [Varname.uvars;fm_strjoin('Vc_',{Bus.names{:}}')];
    Varname.uvars = [Varname.uvars;fm_strjoin('Qgc_',{Bus.names{iBQg}}')];
    Varname.uvars = [Varname.uvars;fm_strjoin(flow,'c',Lf,'-',Lt)];
    Varname.uvars = [Varname.uvars;fm_strjoin(flow,'c',Lt,'-',Lf)];

    Varname.fvars = [Varname.fvars;{[bslash,'lambda_c'];'k_g_c'}];
    Varname.fvars = [Varname.fvars;fm_strjoin(bslash,'theta_{c',{Bus.names{:}}','}')];
    Varname.fvars = [Varname.fvars;fm_strjoin('V_{c',{Bus.names{:}}','}')];
    Varname.fvars = [Varname.fvars;fm_strjoin('Q_{gc',{Bus.names{iBQg}}','}')];
    Varname.fvars = [Varname.fvars;fm_strjoin(flow,'{c',Lf,'-',Lt,'}')];
    Varname.fvars = [Varname.fvars;fm_strjoin(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';