📄 fm_gams.m
字号:
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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -