fm_gams.m

电力系统分析计算程序

function fm_gams
% FM_GAMS initialize and call GAMS to solve
%         several kind of Market Clearing Mechanisms
%
% FM_GAMS
%
%GAMS settings are stored in the structure GAMS, with
%the following fields:
%
%      METHOD   1 -> simple auction
%               2 -> linear OPF (DC power flow)
%               3 -> nonlinear OPF (AC power flow)
%               4 -> nonlinear VSC-OPF
%               5 -> maximum loading condition
%               6 -> continuation OPF
%
%      TYPE     1 -> single period auction
%               2 -> multi period auction
%               3 -> VSC single period auction
%               4 -> VSC multi period auction
%
%see also FM_GAMS.GMS, FM_GAMSFIG and
%structures CPF and OPF for futher settings
%
%Author:    Federico Milano
%Date:      29-Jan-2003
%Update:    01-Feb-2003
%Update:    06-Feb-2003
%Version:   1.0.2
%
%E-mail:    [e-mail]
%Web-site:  http://www.uclm.es/area/gsee/Web/Federico
%
% Copyright (C) 2002-2008 Federico Milano

global DAE OPF CPF GAMS Bus File clpsat
global Path Settings Snapshot Varname
global PV PQ SW Line Shunt jay Varout
global Supply Demand Rmpl Rmpg Ypdp

[u,w] = system('gams');
if u
  fm_disp('GAMS is not properly installed on your system.',2)
  return
end

if ~autorun('PSAT-GAMS Interface',0)
  return
end

if DAE.n
  fm_disp(['Dynamic data are not supported within the PSAT-GAMS interface.'],2)
  return
end

if ~Supply.n,
  fm_disp(['Supply data have to be specified before in order to ', ...
	   'run PSAT-GAMS interface'],2)
  return
end

if ~Demand.n,
  if GAMS.basepg & ~clpsat.init
    Settings.ok = 0;
    uiwait(fm_choice(['Exclude (recommended) base generator powers?']))
    GAMS.basepg = ~Settings.ok;
  end
  noDem = 1;
  Demand = add(Demand,'dummy');
else
  noDem = 0;
end

length(Snapshot.y);
if ~GAMS.basepl
  buspl = Snapshot(1).Pl;
  busql = Snapshot(1).Ql;
  Bus.Pl(:) = 0;
  Bus.Ql(:) = 0;
  PQ = pqzero(PQ,'all');
end
if ~GAMS.basepg
  ploss = Snapshot(1).Ploss;
  Snapshot(1).Ploss = 0;
  buspg = Snapshot(1).Pg;
  busqg = Snapshot(1).Qg;
  Bus.Pg(:) = 0;
  Bus.Qg(:) = 0;
  SW = swzero(SW,'all');
  PV = pvzero(PV,'all');
end

fm_disp
fm_disp('---------------------------------------------------------')
fm_disp(' PSAT-GAMS Interface')
fm_disp('---------------------------------------------------------')
fm_disp

tic
method = GAMS.method;
modelstat = 0;
solvestat = 0;

type = GAMS.type;
omega = GAMS.omega;
if GAMS.method == 6 & GAMS.type ~= 1
  fm_disp(['WARNING: Continuation OPF can be run only with Single' ...
	   ' Period Auctions.'])
  fm_disp('Voltage Stability Constrained OPF will be solved.')
  method = 4;
end
if GAMS.method == 6 & GAMS.flow ~= 1
  fm_disp(['WARNING: Continuation OPF can be run only with Current ' ...
           'Limits.'])
  fm_disp('Current limits in transmission lines will be used.')
  GAMS.flow = 1;
end
if GAMS.type == 3 & GAMS.method ~= 4
  fm_disp(['WARNING: Pareto Set Single Period Auction can be run ' ...
           'only for VSC-OPF.'])
  fm_disp( '         Single Period Auction will be solved.')
  fm_disp
  type = 1;
end
if GAMS.type == 3 & length(GAMS.omega) == 1
  fm_disp(['WARNING: Weighting factor is scalar. ', ...
           'Single Period Auction will be solved.'])
  fm_disp
  type = 1;
end
if GAMS.type == 1 & length(GAMS.omega) > 1
  fm_disp(['WARNING: Weighting factor is a vector. ', ...
           'First omega entry will be used.'])
  fm_disp
  omega = omega(1);
end
if ~rem(GAMS.type,2) & ~Rmpg.n
  type = 1;
  fm_disp(['WARNING: No Ramping data were found. ', ...
           'Single Period Auction will be solved.'])
  fm_disp
end
if GAMS.type == 2 & Rmpg.n & ~Ypdp.n
  type = 4;
  fm_disp(['WARNING: No Power Demand Profile was found. Single ' ...
           'Period Auction with Unit Commitment will be solved.'])
  fm_disp
end

% resetting time vector in case of previous time simulations
if type == 3, Varout.t = []; end

switch method
 case 1, fm_disp(' Simple Auction')
 case 2, fm_disp(' Market Clearing Mechanism')
 case 3, fm_disp(' Standard OPF')
 case 4, fm_disp(' Voltage Stability Constrained OPF')
 case 5, fm_disp(' Maximum Loading Condition')
 case 6, fm_disp(' Continuation OPF')
 case 7, fm_disp(' Congestion Management')
end

switch type
 case 1, fm_disp(' Single-Period Auction')
 case 2, fm_disp(' Multi-Period Auction')
 case 3, fm_disp(' Pareto Set Single-Period Auction')
 case 4, fm_disp(' Single-Period Auction with Unit Commitment')
end

if (GAMS.flatstart | isempty(Snapshot)) & GAMS.method > 2
  DAE.y(Bus.a) = getzeros(Bus);
  DAE.y(Bus.v) = getones(Bus);
else
  DAE.y = Snapshot(1).y;
end

% ------------------------------------------------------------
% Parameter definition
% ------------------------------------------------------------

% dimensions
nBus = int2str(Bus.n);
nQg = int2str(PV.n+SW.n);
nBusref = int2str(SW.refbus);

[nSW,SW_idx,ksw] = gams(SW);
[nPV,PV_idx,kpv] = gams(PV);
[nLine,L,Li,Lj,Gh,Bh,Ghc,Bhc] = gams(Line,method);
[Gh,Bh,Ghc,Bhc] = gams(Shunt,method,Gh,Bh,Ghc,Bhc);
[nPd,Pd_idx,D] = gams(Demand);
[nPs,Ps_idx,S] = gams(Supply,type);
[nH,Ch] = gams(Ypdp,type);

% indexes
iBPs = Supply.bus;
iBPd = Demand.bus;
iBQg = [SW.bus; PV.bus];

% Fixed powers
Pg0 = Bus.Pg;
Pl0 = Bus.Pl;
Ql0 = Bus.Ql;

% Generator reactive powers and associated limits
Qg0 = getzeros(Bus);
Qg0(iBQg) = Bus.Qg(iBQg);
[Qgmax,Qgmin] = fm_qlim('all');

% Voltage limits
V0 = DAE.y(Bus.v);
t0 = DAE.y(Bus.a);
[Vmax,Vmin] = fm_vlim(1.5,0.2);

% ------------------------------------------------------------
% Data structures
% ------------------------------------------------------------

X.val = [V0,t0,Pg0,Qg0,Pl0,Ql0,Qgmax,Qgmin,Vmax,Vmin,ksw,kpv];
lambda.val = [GAMS.lmin(1),GAMS.lmax(1),GAMS.omega(1),GAMS.line];

X.labels = {cellstr(num2str(Bus.a)), ...
	    {'V0','t0','Pg0','Qg0','Pl0','Ql0', ...
	     'Qgmax','Qgmin','Vmax','Vmin','ksw','kpv'}};
lambda.labels = {'lmin','lmax','omega','line'};

X.name = 'X';
lambda.name = 'lambda';

% ------------------------------------------------------------
% 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 | ispc, 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.y(Bus.a)]';
      V(i,:) = [DAE.y(Bus.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.y(Bus.a)]';
      V(i,:) = [DAE.y(Bus.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