fm_gams.m

电力系统的psat

  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);
  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.*tgphi,Bus.n,1)+Bus.Ql)*MVA);
PayS = -ro(1:Bus.n).*PG;
PayD = ro(1:Bus.n).*PL;
ISOPay = -sum(ro(1:Bus.n).*DAE.glfp*MVA);

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

  fm_disp
  fm_disp(' Power Supplies')
  fm_disp(' ---------------------------------------------------------------')
  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([Bus.con(Supply.bus,1),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([Bus.con(Supply.bus,1),Ps*MVA,Psmax*MVA,Psmin*MVA])
  end
  fm_disp
  fm_disp(' Power Demands')
  fm_disp(' ---------------------------------------------------------------')
  if GAMS.method == 7
    fm_disp({'Bus','Pd','Pd max','Pd min','dPd_up','dPd_dw'})
    fm_disp({'<i>','[MW]','[MW]','[MW]'})
    fm_disp([Bus.con(Demand.bus,1),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([Bus.con(Demand.bus,1),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([Bus.con(iBQg,1),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([Bus.con(iBQg,1),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([Bus.con(:,1),DAE.V, DAE.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([Bus.con(:,1),ro(1:Bus.n), NCP, PayS, PayD])
  else
    fm_disp({'Bus','LMP','Pay S','Pay D'})
    fm_disp({'<i>','[$/MWh]','[$/h]','[$/h]'})
    fm_disp([Bus.con(:,1),ro(1:Bus.n), 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([Bus.con(:,1),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.from, Line.to,Pij, ...
             Pijmax,abs((-abs(Pij)+Pijmax)), ...
             Pji,Pijmax,abs((-abs(Pji)+Pijmax))],1)
    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'},1)
     case 2,
      fm_disp({'From Bus','To Bus','Pijc','Pijcmax', ...
               'Pijc margin','Pjic','Pjicmax','Pjic margin'},1)
     case 3,
      fm_disp({'From Bus','To Bus','Sijc','Sijcmax', ...
               'Sijc margin','Sjic','Sjicmax',['Sjic margin']},1)
    end
    fm_disp({'<i>','<j>','[p.u.]','[p.u.]', ...
	     '[p.u.]','[p.u.]','[p.u.]','[p.u.]'})
    if GAMS.flow
      fm_disp([Line.from, Line.to,Pijc, ...
               Pijmax,abs((-abs(Pijc)+Pijmax)), ...
               Pjic,Pijmax,abs((-abs(Pjic)+Pijmax))],1)
      fm_disp
    end
  end

  fm_disp
  fm_disp(' Totals')
  fm_disp(' ---------------------------------------------------------------')
  if GAMS.method >= 4,
    fm_disp([' omega = ',num2str(omega(1))],1)
    fm_disp([' lambda_c = ',num2str(lambdac),' [p.u.]'],1)
    fm_disp([' kg = ',num2str(kg),' [p.u.]'],1)
  end
  if GAMS.method == 1
    fm_disp([' Market Clearing Price = ',num2str(MCP),' [$/MWh]'])
  end
  fm_disp([' Total Losses = ',num2str(1e-5*round((sum(DAE.glfp))* ...
                                                 1e5)*MVA),' [MW]'],1)
  fm_disp([' Bid Losses = ',num2str(1e-5*round((sum(DAE.glfp)- ...
                                                Snapshot(1).Ploss)* ...
                                               1e5)*MVA),' [MW]'],1)
  fm_disp([' Total demand = ',num2str(sum(Pd)*MVA),' [MW]'],1)
  fm_disp([' Total Transaction Level = ', ...
           fvar(sum(Pd)*MVA+sum(PQ.con(:,4))*MVA,8),' [MW]']);
  if GAMS.method == 4 | GAMS.method == 6
    fm_disp([' Maximum Loading Condition = ', ...
             fvar((1+lambdac)*(sum(Pd)*MVA+sum(PQ.con(:,4))*MVA),8),' [MW]']);
    fm_disp([' Available Loading Capability = ', ...
             fvar(lambdac*(sum(Pd)*MVA+sum(PQ.con(:,4))*MVA),8),' [MW]']);
  end
  if GAMS.method == 5
    fm_disp([' Maximum Loading Condition = ', ...
             fvar(lambdac*(sum(Pd)*MVA+sum(PQ.con(:,4))*MVA),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.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;

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

global Settings Bus Line PQ PV Demand Supply GAMS

PQold = PQ;
PVold = PV;
Busold = Bus;
Demand.con(:,7) = Pd;
Supply.con(:,6) = Ps;
for i = 1:Demand.n
  a = find(PQ.bus == Demand.bus(i));
  PQ.con(a,4) = PQ.con(a,4) + Demand.con(i,7);
  PQ.con(a,5) = PQ.con(a,5) + Demand.con(i,4)* ...
      Demand.con(i,7)/Demand.con(i,3);
end
for i = 1:Supply.n
  a = find(PV.bus == Supply.bus(i));
  if ~isempty(a)
    PV.con(a,4) = PV.con(a,4) + Supply.con(i,6);
  end
end
show_old = Settings.show;
Settings.show = 0;
Settings.locksnap = 1;
fm_spf
Settings.locksnap = 0;
Settings.show = show_old;
gams_flow = GAMS.flow;
if ~gams_flow
  gams_flow = 1;
end
[Pij, Jij, Hij, Pji, Jji, Hji] = ...
    fm_flows(gams_flow,ones(Line.n,1),ones(Line.n,1));
Pij = sqrt(Pij);
Pji = sqrt(Pji);
Qg = Bus.Qg(iBQg);
Bus = Busold;
PQ = PQold;
PV = PVold;


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

global DAE SW GAMS Line Bus

Vold = DAE.V;
aold = DAE.a;
J11old = DAE.J11;
J12old = DAE.J12;
J21old = DAE.J21;
J22old = DAE.J22;

DAE.V = V;
DAE.a = a;
fm_lf(2)
Jlfv = [DAE.J11,DAE.J12;DAE.J21,DAE.J22];
Jlfv(SW.bus,:) = [];
Jlfv(:,SW.bus) = [];
gams_flow = GAMS.flow;
if ~gams_flow
  gams_flow = 1;
end
[Fij, Jij, Hij, Fji, Jji, Hji] = ...
    fm_flows(gams_flow,ones(Line.n,1),ones(Line.n,1));
dH_dtV = Jij'*mFij + Jji'*mFji + [zeros(Bus.n,1); ...
                    mV];
dH_dtV(SW.bus,:) = [];
NCP = Jlfv'\dH_dtV;
NCP = [NCP(1:SW.bus-1);0;NCP(SW.bus:end)];
NCP = NCP([1:Bus.n]);

DAE.V = Vold;
DAE.a = aold;
DAE.J11 = J11old;
DAE.J12 = J12old;
DAE.J21 = J21old;
DAE.J22 = J22old;

% ==========================================================================
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})
    if Settings.octave
      fprintf(fid1,'$setglobal %s ''%s''\n',argn(i,:), ...
              varargin{i});
    else
      fprintf(fid1,'$setglobal %s ''%s''\n',inputname(i), ...
              varargin{i});
    end
  elseif isnumeric(varargin{i})
    if Settings.octave
      fprintf(fid2,'$kill %s\n',argn(i,:));
      if length(varargin{i}) == 1
        fprintf(fid2,'scalar %s /%f/;\n',argn(i,:),varargin{i});
      else
        fprintf(fid2,'parameter %s /\n',argn(i,:));
        [x,y,v] = find(varargin{i});
        fprintf(fid2,'%d.%d %f\n',[x y v]');
        fprintf(fid2,'/;\n');
      end
    else
      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
    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 = strcat(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{x(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;
if Settings.octave
  [result,status] = system(['gams \"',Path.psat,varargin{1},'\"']);
elseif Settings.hostver < 6.1
  if isunix
    [status,result] = unix(['gams ',varargin{1}]);  
  else
    [status,result] = dos(['gams ',varargin{1}]);  
  end
else  
  [status,result] = system(['gams ',varargin{1}]);
end
fm_disp([' GAMS routine completed in ',num2str(etime(clock,t0)),' s'])
if status == 1
  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