fm_gams.m

基于PSAT 软件的多目标最优潮流计算用于中小型电力系统的分析和管理

  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(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(' ---------------------------------------------------------------')
  [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([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(' ---------------------------------------------------------------')
  [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([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))])
    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.from, Line.to,Pijc, ...
               Pijmax,abs((-abs(Pijc)+Pijmax)), ...
               Pjic,Pijmax,abs((-abs(Pjic)+Pijmax))])
      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(DAE.glfp)*1e5)*MVA;
  bid_loss = 1e-5*round((sum(DAE.glfp)-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;
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;
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

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})
    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 = 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{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;
[status,result] = system(['gams ',varargin{1}]);
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