fm_gams.m

电力系统分析计算程序

  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.y(Bus.a) = ac;
  DAE.y(Bus.v) = Vc;
  Line = gcall(Line);
  glfpc = Line.p;
  glfqc = Line.q;
end

if GAMS.method >= 3
  DAE.y(Bus.a) = a;
  DAE.y(Bus.v) = V;
  Line = gcall(Line);
  Qg = Qg(iBQg);
end

if GAMS.method == 1
  ro = MCP*getones(Bus);
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);
  PL = full(sparse(iBPd,1,Pd,Bus.n,1)*MVA);
end
QG = full(sparse(iBQg,1,Qg,Bus.n,1)*MVA);
QL = full((sparse(iBPd,1,Pd.*tanphi(Demand),Bus.n,1)+Bus.Ql)*MVA);
PayS = -ro(Bus.a).*PG;
PayD = ro(Bus.a).*PL;
ISOPay = -sum(ro(Bus.a).*Line.p*MVA);

if (Settings.showlf | GAMS.show) & clpsat.showopf

  fm_disp
  fm_disp(' Power Supplies')
  fm_disp(' ---------------------------------------------------------------')
  [Psmax,Psmin] = plim(Supply);
  if GAMS.method == 7
    fm_disp({'Bus','Ps','Ps max','Ps min','dPs_up','dPs_dw'})
    fm_disp({'<i>','[MW]','[MW]','[MW]','[MW]','[MW]'})
    fm_disp([getidx(Bus,Supply.bus),Ps*MVA,Psmax*MVA,Psmin*MVA,dPSup*MVA,dPSdw*MVA])
  else
    fm_disp({'Bus','Ps','Ps max','Ps min'})
    fm_disp({'<i>','[MW]','[MW]','[MW]'})
    fm_disp([getidx(Bus,Supply.bus),Ps*MVA,Psmax*MVA,Psmin*MVA])
  end
  fm_disp
  fm_disp(' Power Demands')
  fm_disp(' ---------------------------------------------------------------')
  [Pdmax,Pdmin] = plim(Demand);
  if GAMS.method == 7
    fm_disp({'Bus','Pd','Pd max','Pd min','dPd_up','dPd_dw'})
    fm_disp({'<i>','[MW]','[MW]','[MW]'})
    fm_disp([getidx(Bus,Demand.bus),Pd*MVA,Pdmax*MVA,Pdmin*MVA,dPDup*MVA,dPDdw*MVA])
  else
    fm_disp({'Bus','Pd','Pd max','Pd min'})
    fm_disp({'<i>','[MW]','[MW]','[MW]'})
    fm_disp([getidx(Bus,Demand.bus),Pd*MVA,Pdmax*MVA,Pdmin*MVA])
  end
  fm_disp
  fm_disp(' Generator Reactive Powers')
  fm_disp(' ---------------------------------------------------------------')
  if GAMS.method == 4 | GAMS.method == 6
    fm_disp({'Bus','Qg','Qgc','Qg max','Qg min'})
    fm_disp({'<i>','[MVar]','[MVar]','[MVar]','[MVar]'})
    fm_disp([getidx(Bus,iBQg),Qg*MVA,Qgc(iBQg)*MVA,Qgmax*MVA,Qgmin*MVA])
  else
    fm_disp({'Bus','Qg','Qg max','Qg min'})
    fm_disp({'<i>','[MVar]','[MVar]','[MVar]'})
    fm_disp([getidx(Bus,iBQg),Qg*MVA,Qgmax*MVA,Qgmin*MVA])
  end
  fm_disp
  fm_disp(' Power Flow Solution')
  fm_disp([' ----------------------------------------------------' ...
           '-----------'])
  fm_disp({'Bus','V','theta','PG','PL','QG','QL'})
  fm_disp({'<i>','[p.u.]','[rad]','[MW]','[MW]','[MVar]','[MVar]'})
  fm_disp([getidx(Bus,0),DAE.y(Bus.v),DAE.y(Bus.a),PG,PL,QG,QL])

  fm_disp
  fm_disp(' Prices and Pays')
  fm_disp([' ----------------------------------------------------' ...
           '-----------'])
  if GAMS.method == 3 | GAMS.method == 4 | GAMS.method == 6
    fm_disp({'Bus','LMP','NCP','Pay S','Pay D'})
    fm_disp({'<i>','[$/MWh]','[$/MWh]','[$/h]','[$/h]'})
    fm_disp([getidx(Bus,0),ro(Bus.a), NCP, PayS, PayD])
  else
    fm_disp({'Bus','LMP','Pay S','Pay D'})
    fm_disp({'<i>','[$/MWh]','[$/h]','[$/h]'})
    fm_disp([getidx(Bus,0),ro(Bus.a), PayS, PayD])
  end
  if GAMS.method == 4 | GAMS.method == 6
    fm_disp
    fm_disp(' "Critical" Power Flow Solution')
    fm_disp(' ---------------------------------------------------------------')
    fm_disp({'Bus','Vc','thetac','PGc','PLc','QGc','QLc'})
    fm_disp({'<i>','[p.u.]','[rad]','[MW]','[MW]','[MVar]', ...
             '[MVar]'})
    PG = (1+lambdac+kg)*PG;
    PL = (1+lambdac)*PL;
    QL = (1+lambdac)*QL;
    fm_disp([getidx(Bus,0),Vc,ac,PG,PL,Qgc*MVA,QL])
  end
  fm_disp
  if GAMS.flow
    fm_disp(' Flows on Transmission Lines')
    fm_disp(' ---------------------------------------------------------------')
    switch GAMS.flow
     case 1,
      fm_disp({'From Bus','To Bus','Iij','Iijmax', ...
               'Iij margin','Iji','Ijimax','Iji margin'},1)
     case 2,
      fm_disp({'From Bus','To Bus','Pij','Pijmax', ...
               'Pij margin','Pji','Pjimax','Pji margin'},1)
     case 3,
      fm_disp({'From Bus','To Bus','Sij','Sijmax', ...
               'Sij margin','Sji','Sjimax','Sji margin'},1)
    end
    fm_disp({'<i>','<j>','[p.u.]','[p.u.]', ...
	     '[p.u.]','[p.u.]','[p.u.]','[p.u.]'})
    fm_disp([Line.fr, Line.to,Pij, ...
             L.val(:,5),abs((-abs(Pij)+L.val(:,5))), ...
             Pji,L.val(:,5),abs((-abs(Pji)+L.val(:,5)))])
    fm_disp
  else
    fm_disp('Flow limits are disabled.')
  end
  if GAMS.method == 4 | GAMS.method == 6
    fm_disp(' Flows on Transmission Lines of the "Critical" System')
    fm_disp(' ---------------------------------------------------------------')
    switch GAMS.flow
     case 1,
      fm_disp({'From Bus','To Bus','Iijc','Iijcmax', ...
               'Iijc margin','Ijic','Ijicmax','Ijic margin'})
     case 2,
      fm_disp({'From Bus','To Bus','Pijc','Pijcmax', ...
               'Pijc margin','Pjic','Pjicmax','Pjic margin'})
     case 3,
      fm_disp({'From Bus','To Bus','Sijc','Sijcmax', ...
               'Sijc margin','Sjic','Sjicmax',['Sjic margin']})
    end
    fm_disp({'<i>','<j>','[p.u.]','[p.u.]', ...
	     '[p.u.]','[p.u.]','[p.u.]','[p.u.]'})
    if GAMS.flow
      fm_disp([Line.fr, Line.to,Pijc, ...
               L.val(:,5),abs((-abs(Pijc)+L.val(:,5))), ...
               Pjic,L.val(:,5),abs((-abs(Pjic)+L.val(:,5)))])
      fm_disp
    end
  end

  fm_disp
  fm_disp(' Totals')
  fm_disp(' ---------------------------------------------------------------')
  if GAMS.method >= 4,
    fm_disp([' omega = ',num2str(omega(1))])
    fm_disp([' lambda_c = ',num2str(lambdac),' [p.u.]'])
    fm_disp([' kg = ',num2str(kg),' [p.u.]'])
  end
  if GAMS.method == 1
    fm_disp([' Market Clearing Price = ',num2str(MCP),' [$/MWh]'])
  end

  total_loss = 1e-5*round(sum(Line.p)*1e5)*MVA;
  bid_loss = 1e-5*round((sum(Line.p)-Snapshot(1).Ploss)*1e5)*MVA;

  fm_disp([' Total Losses = ',num2str(total_loss),' [MW]'])
  fm_disp([' Bid Losses = ',num2str(bid_loss),' [MW]'])
  fm_disp([' Total demand = ',num2str(sum(Pd)*MVA),' [MW]'])
  fm_disp([' Total Transaction Level = ', ...
           fvar(sum(Pd)*MVA+TPQ,8),' [MW]']);
  if GAMS.method == 4 | GAMS.method == 6
    fm_disp([' Maximum Loading Condition = ', ...
             fvar((1+lambdac)*(sum(Pd)*MVA+TPQ),8),' [MW]']);
    fm_disp([' Available Loading Capability = ', ...
             fvar(lambdac*(sum(Pd)*MVA+TPQ),8),' [MW]']);
  end
  if GAMS.method == 5
    fm_disp([' Maximum Loading Condition = ', ...
             fvar(lambdac*(sum(Pd)*MVA+TPQ),8),' [MW]']);
  end
  fm_disp([' IMO Pay = ',num2str(ISOPay),' [$/h]']);
  fm_disp

end
fm_disp(' ---------------------------------------------------------------')
fm_disp([' Check file ',Path.psat,'fm_gams.lst for GAMS report.'])
gams_mstat(modelstat)
gams_sstat(solvestat)
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 = restore(Demand); end

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;

% ===============================================================
function [Pij,Pji,Qg] = updatePF(Pd,Ps,iBQg)
% Power FLow Solution with the current simple auction solution
% ===============================================================

global Settings Bus Line PQ PV SW Demand Supply GAMS

Busold = Bus;
Demand = pset(Demand,Pd);
Supply = pset(Supply,Ps);

pg = SW.pg;
pqsum(Demand,1);
pgsum(Supply,1);

show_old = Settings.show;
Settings.show = 0;
Settings.locksnap = 1;
fm_spf
Settings.locksnap = 0;
Settings.show = show_old;
[Pij,Pji] = flows(Line,max(GAMS.flow,1));
Qg = Bus.Qg(iBQg);
Bus = Busold;
SW = setpg(SW,'all',pg);
PQ = restore(PQ);
PV = pvreset(PV,'all');

% ==========================================================================
function NCP = compNCP(V,a,mV,mFij,mFji)
% Nodal Congestion Prices
% ==========================================================================

global DAE SW GAMS Line Bus

yold = DAE.y;
Gyold = DAE.Gy;

DAE.y(Bus.a) = a;
DAE.y(Bus.v) = V;
Gycall(Line)
fm_setgy(SW.refbus)
[Fij,Jij,Fji,Jji] = fjh2(Line,max(GAMS.flow,1));
dH_dtV = Jij'*mFij + Jji'*mFji + [getzeros(Bus);mV];
dH_dtV(SW.refbus,:) = 0;
NCP = DAE.Gy'\dH_dtV;
NCP = NCP(Bus.a);

DAE.y = yold;
DAE.Gy = Gyold;

% ==========================================================================
function varargout = psatgams(varargin)
% PSAT-GAMS interface
% ==========================================================================

global Settings Path

% writing GAMS input data
%---------------------------------------------------------------------
fid1 = fopen('psatglobs.gms','wt+');
fid2 = fopen('psatdata.gms','wt+');
fprintf(fid2,'%s\n','$onempty');

for i = 2:nargin
  if ischar(varargin{i})
    fprintf(fid1,'$setglobal %s ''%s''\n',inputname(i), ...
            varargin{i});
  elseif isnumeric(varargin{i})
    fprintf(fid2,'$kill %s\n',inputname(i));
    if length(varargin{i}) == 1
      fprintf(fid2,'scalar %s /%f/;\n',inputname(i),varargin{i});
    else
      fprintf(fid2,'parameter %s /\n',inputname(i));
      [x,y,v] = find(varargin{i});
      fprintf(fid2,'%d.%d %f\n',[x y v]');
      fprintf(fid2,'/;\n');
    end
  elseif isstruct(varargin{i})
    labels = varargin{i}.labels;
    fprintf(fid2,'$kill %s\n',varargin{i}.name);
    fprintf(fid2,'parameter %s /\n',varargin{i}.name);
    [x,y,v] = find(varargin{i}.val);
    if iscell(labels{1})
      %a = fm_strjoin(labels{1}(x),'.',labels{2}(y)',[blanks(length(x))',num2str(v)]);
      %fprintf(fid2,'%s\n',a{:});
      for j = 1:length(x)
	fprintf(fid2,'%s.%s %f\n',labels{1}{x(j)},labels{2}{y(j)},v(j));
      end
    else
      for j = 1:length(x)
	fprintf(fid2,'%s %f\n',labels{y(j)},v(j));
      end
    end
    fprintf(fid2,'/;\n');
  end
end
fprintf(fid2,'%s\n','$offempty');
fclose(fid1);
fclose(fid2);

% Lauching GAMS
%---------------------------------------------------------------------
status = 0;
t0 = clock;
disp(['gams ',varargin{1},' -error=PSAT'])
[status,result] = system(['gams ',varargin{1}]);
fm_disp([' GAMS routine completed in ',num2str(etime(clock,t0)),' s'])
if status
  fm_disp(result)
  return
end

% Reading GAMS output
%---------------------------------------------------------------------
nout = 0;
EPS = eps;
clear psatsol
psatsol
if nout < nargout
  for i = nout+1:nargout
    varargout{i} = [];
  end
end

if nout > nargout
  varargout(nargout+1:nout) = [];
end

%---------------------------------------------------------------------
function gams_mstat(status)

if isempty(status), return, end

switch status
 case 0, fm_disp(' GAMS model status: not available')
 case 1, fm_disp(' GAMS model status: optimal')
 case 2, fm_disp(' GAMS model status: locally optimal')
 case 3, fm_disp(' GAMS model status: unbounded')
 case 4, fm_disp(' GAMS model status: infeasible')
 case 5, fm_disp(' GAMS model status: locally infeasible')
 case 6, fm_disp(' GAMS model status: intermediate infeasible')
 case 7, fm_disp(' GAMS model status: intermediate non-optimal')
 case 8, fm_disp(' GAMS model status: integer solution')
 case 9, fm_disp(' GAMS model status: intermediate non-integer')
 case 10, fm_disp(' GAMS model status: integer infeasible')
 case 11, fm_disp(' GAMS model status: ???')
 case 12, fm_disp(' GAMS model status: error unknown')
 case 13, fm_disp(' GAMS model status: error no solution')
 otherwise, fm_disp(' GAMS model status: unknown model status')
end

%---------------------------------------------------------------------
function gams_sstat(status)

if isempty(status), return, end

switch status
 case 0, fm_disp(' GAMS solver status: not available')
 case 1, fm_disp(' GAMS solver status: normal completion')
 case 2, fm_disp(' GAMS solver status: iteration interrupt')
 case 3, fm_disp(' GAMS solver status: resource interrupt')
 case 4, fm_disp(' GAMS solver status: terminated by solver')
 case 5, fm_disp(' GAMS solver status: evaluation error limit')
 case 6, fm_disp(' GAMS solver status: unknown')
 case 7, fm_disp(' GAMS solver status: ???')
 case 8, fm_disp(' GAMS solver status: error preprocessor error')
 case 9, fm_disp(' GAMS solver status: error setup failure')
 case 10, fm_disp(' GAMS solver status: error solver failure')
 case 11, fm_disp(' GAMS solver status: error internal solver error')
 case 12, fm_disp(' GAMS solver status: error post-processor error')
 case 13, fm_disp(' GAMS solver status: error system failure')
 otherwise, fm_disp(' GAMS solver status: unknown solver status')
end