$OnText
This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
Written by Behrang Shirizadeh, Quentin Perrier and Philippe Quirion, October 2020
$Offtext
*-------------------------------------------------------------------------------
*                                Defining the sets
*-------------------------------------------------------------------------------
sets     h                                               /0*8759/
         first(h)        'first hour'
         last(h)         'last hour'
         m               'month'                         /jan, feb, mar, apr, may, jun, jul, aug, sep, oct, nov, dec/
         tec             'technology'                    /offshore, onshore, pv, river, lake, biogas, gas, phs, battery, methanation/
         gen(tec)        'power plants'                  /offshore, onshore, pv, river, lake, biogas, gas/
         vre(tec)        'variable tecs'                 /offshore, onshore, pv/
         ncomb(tec)      'non-combustible generation'    /offshore, onshore, pv, river, lake, phs, battery/
         comb(tec)       'combustible generation techs'  /biogas, methanation/
         str(tec)        'storage technologies'          /phs, battery, methanation/
         frr(tec)        'technologies for upward FRR'   /lake, phs, battery, gas/
;
first(h) = ord(h)=1;
last(h) = ord(h)=card(h);
alias(h,hh);
*-------------------------------------------------------------------------------
*                                Inputs
*-------------------------------------------------------------------------------
parameter month(h)  /0*743 1, 744*1439 2, 1440*2183 3, 2184*2903 4
                    2904*3647 5, 3648*4367 6, 4368*5111 7, 5112*5855 8
                    5856*6575 9, 6576*7319 10, 7320*8039 11, 8040*8759 12/
$Offlisting
parameter load_factor(vre,h) 'Production profiles of VRE'
/
$ondelim
$include  inputs/vre_profiles2006.csv
$offdelim
/;
parameter demand(h) 'demand profile in each hour in GW'
/
$ondelim
$include inputs/demand2050_ademe.csv
$offdelim
/;
Parameter lake_inflows(m) 'monthly lake inflows in GWh'
/
$ondelim
$include  inputs/lake_inflows.csv
$offdelim
/ ;
parameter gene_river(h) 'hourly run of river power generation in GWh'
/
$ondelim
$include  inputs/run_of_river.csv
$offdelim
/ ;
parameter epsilon(vre) 'additional FRR requirement for variable renewable energies because of forecast errors'
/
$ondelim
$include  inputs/reserve_requirements.csv
$offdelim
/ ;
parameter capa_ex(tec) 'existing capacities of the technologies by December 2017 in GW'
/
$ondelim
$include  inputs/existing_capas.csv
$offdelim
/ ;
parameter capa_max(tec) 'maximum capacities of the technologies in GW'
/
$ondelim
$include  inputs/max_capas.csv
$offdelim
/ ;
parameter capex(tec) 'annualized power capex cost in M€/GW/year'
/
$ondelim
$include  inputs/annuities.csv
$offdelim
/ ;
parameter capex_en(str) 'annualized energy capex cost of storage technologies in M€/GWh/year'
/
$ondelim
$include  inputs/str_annuities.csv
$offdelim
/ ;
parameter fOM(tec) 'annualized fixed operation and maintenance costs M€/GW/year'
/
$ondelim
$include  inputs/fO&M.csv
$offdelim
/ ;
Parameter vOM(tec) 'Variable operation and maintenance costs in M€/GWh'
/
$ondelim
$include  inputs/vO&M.csv
$offdelim
/ ;

parameter s_capex(str) 'charging related annuity of storage in M€/GW/year' /PHS 0, battery 0, methanation 84.16086/;
parameter s_opex(str)    'charging related fOM of storage in M€/GW/year'   /PHS 0, battery 0, methanation 59.25/;
parameter eta_in(str) 'charging efifciency of storage technologies' /PHS 0.95, battery 0.9, methanation 0.59/;
parameter eta_out(str) 'discharging efficiency of storage technolgoies' /PHS 0.9, battery 0.95, methanation 0.45/;
scalar pump_capa 'pumping capacity in GW' /9.3/;
scalar max_phs 'maximum volume of energy can be stored in PHS reservoir in TWh' /0.18/;
scalar max_biogas 'maxium energy can be generated by biogas in TWh' /15/;
scalar load_uncertainty 'uncertainty coefficient for hourly demand' /0.01/;
scalar delta 'load variation factor'     /0.1/;
parameter fixed_costs(tec) 'yearly fixed cost of each tec in M€/GW/year' ;
fixed_costs(tec) = capex(tec) + fOM(tec);
*-------------------------------------------------------------------------------
*                                Model
*-------------------------------------------------------------------------------
variables        GENE(tec,h)     'hourly energy generation in TWh'
                 CAPA(tec)       'overal yearly installed capacity in GW'
                 STORAGE(str,h)  'hourly electricity input of battery storage GW'
                 S(str)          'charging power capacity of each storage technology'
                 STORED(str,h)   'energy stored in each storage technology in GWh'
                 CAPACITY(str)   'energy volume of storage technologies in GWh'
                 RSV(frr,h)      'required upward frequency restoration reserve in GW'
                 COST            'final investment cost in b€'

positive variables GENE(tec,h),CAPA(tec),STORAGE(str,h), S(str),STORED(str,h),CAPACITY(str),RSV(frr,h);

equations        gene_vre        'variables renewable profiles generation'
                 gene_capa       'capacity and genration relation for technologies'
                 combustion      'the relationship of combustible technologies'
                 capa_frr        'capacity needed for the secondary reserve requirements'
                 storing         'the definition of stored energy in the storage options'
                 storage_const   'storage in the first hour is equal to the storage in the last hour'
                 lake_res        'constraint on water for lake reservoirs'
                 stored_cap      'maximum energy that is stored in storage units'
                 storage_capa1   'the capacity with hourly charging relationship of storage'
                 storage_capa2   'storage power limit'
                 biogas_const    'maximum energy can be produced by biogas'
                 reserves        'FRR requirement'
                 adequacy        'supply/demand relation'
                 obj             'the final objective function which is COST';

gene_vre(vre,h)..                GENE(vre,h)             =e=     CAPA(vre)*load_factor(vre,h);
gene_capa(tec,h)..               CAPA(tec)               =g=     GENE(tec,h);
combustion(h)..                  GENE('gas',h)           =e=     sum(comb,GENE(comb,h));
capa_frr(frr,h)..                CAPA(frr)               =g=     GENE(frr,h) + RSV(frr,h);
storing(h,h+1,str)..             STORED(str,h+1)         =e=     STORED(str,h) + STORAGE(str,h)*eta_in(str) - GENE(str,h)/eta_out(str);
storage_const(str,first,last)..  STORED(str,first)       =e=     STORED(str,last)+STORAGE(str,last)*eta_in(str)-GENE(str,last)/eta_out(str);
lake_res(m)..                    lake_inflows(m)         =g=     sum(h$(month(h) = ord(m)),GENE('lake',h))/1000;
stored_cap(str,h)..              STORED(str,h)           =l=     CAPACITY(str);
storage_capa1(str,h)..           S(str)                  =g=     STORAGE(str,h);
storage_capa2(str)..             S(str)                  =l=     CAPA(str);
biogas_const..                   sum(h,GENE('biogas',h)) =l=     max_biogas*1000;
reserves(h)..                    sum(frr, RSV(frr,h))    =e=     sum(vre,epsilon(vre)*CAPA(vre))+ demand(h)*load_uncertainty*(1+delta);
adequacy(h)..                    sum(ncomb,GENE(ncomb,h))+GENE('gas',h)    =g=     demand(h)+sum(str,STORAGE(str,h)) ;
obj..                            COST                    =e=     (sum(tec,(CAPA(tec)-capa_ex(tec))*capex(tec))+ sum(str,CAPACITY(str)*capex_en(str))+sum(tec,(CAPA(tec)*fOM(tec)))+ sum(str,S(str)*(s_capex(str)+s_opex(str))) + sum((tec,h),GENE(tec,h)*vOM(tec)))/1000;
*-------------------------------------------------------------------------------
*                                Initial and fixed values
*-------------------------------------------------------------------------------
CAPA.up('phs') = pump_capa;
CAPA.fx('river')= capa_ex('river');
CAPA.up('lake') = 13;
CAPA.lo('lake') = capa_ex('lake');
GENE.up('river',h) = gene_river(h)*capa_ex('river');
S.up('phs') = pump_capa;
CAPACITY.up('phs') = max_phs*1000;
CAPA.up('offshore') = capa_max('offshore');
CAPA.up('onshore') = capa_max('onshore');
CAPA.up('pv') = capa_max('pv');
*-------------------------------------------------------------------------------
*                                Model options
*-------------------------------------------------------------------------------
model RES_FR /all/;
*-------------------------------------------------------------------------------
option solvelink=0;
option RESLIM = 1000000;
option lp=CPLEX;
option Savepoint=1;
option solveopt = replace;
option limcol = 0;
option limrow = 0;
option SOLPRINT = OFF;
*-------------------------------------------------------------------------------
*                                Solve statement
*-------------------------------------------------------------------------------
$If exist RES_FR_p.gdx execute_loadpoint 'RES_FR_p';
Solve RES_FR using lp minimizing COST;
*-------------------------------------------------------------------------------
*                                Display statement
*-------------------------------------------------------------------------------
parameter sumdemand      'the whole demand per year in TWh';
sumdemand =  sum(h,demand(h))/1000;
parameter gene_tec(tec) 'Overall yearly energy generated by the technology in TWh';
gene_tec(tec) = sum(h,GENE.l(tec,h))/1000;
parameter sumgene        'the whole generation per year in TWh';
sumgene = sum((gen,h),GENE.l(gen,h))/1000 - gene_tec('gas');
parameter sum_FRR 'the whole yearly energy budgeted for reserves in TWh';
sum_FRR = sum((h,frr),RSV.l(frr,h))/1000;
parameter reserve(frr) 'capacity allocated for reserve from each FRR tech in GW';
reserve(frr) = smax(h,RSV.l(frr,h));
Parameter lcoe(gen);
lcoe(gen) = (CAPA.l(gen)*(fOM(gen)+capex(gen))+ gene_tec(gen)*vOM(gen)*1000)/gene_tec(gen);
parameter lcoe_sys1;
lcoe_sys1 = cost.l*1000/sumgene;
parameter lcoe_sys2;
lcoe_sys2 = cost.l*1000/sumdemand;
parameter lcos(str);
lcos(str) = (CAPA.l(str)*(fOM(str)+capex(str))+ gene_tec(str)*vOM(str) + S.l(str)*(s_capex(str)+s_opex(str))+ CAPACITY.l(str)*capex_en(str))/gene_tec(str);
parameter cf(gen) 'load factor of generation technologies';
cf(gen) = gene_tec(gen)*1000/(8760*CAPA.l(gen));
parameter str_loss 'yearly storage related loss in % of power production';
str_loss = (sum((str,h),STORAGE.l(str,h))-sum(str,gene_tec(str)*1000))/(sumgene*10);
parameter lc 'load curtailment of the network';
lc = ((sumgene - sumdemand)*100/sumgene) - str_loss;
parameter spot_price(h) 'marginal cost'    ;
spot_price(h) = 1000000*adequacy.m(h);
parameter marginal_cost 'average value over the year of spot price in €/MWh';
marginal_cost = sum(h,spot_price(h))/8760;
parameter gas_price;
gas_price = -1000000*biogas_const.m;
*-------------------------------------------------------------------------------
display cost.l;
display capa.l;
display sumdemand;
display sumgene;
display gene_tec;
display sum_FRR;
display lcoe;
display lcoe_sys1;
display lcoe_sys2;
display lcos;
display CAPACITY.l;
display cf;
display lc;
display str_loss;
display marginal_cost;
display S.l;
display gas_price;
*-------------------------------------------------------------------------------
*                                Output
*-------------------------------------------------------------------------------
$Ontext
two main output files;
The .txt file just to have a summary and general idea of the key numbers
The .csv file to have a fine output with hourly data for final data processing and analysis
$Offtext
*-------------------------------------------------------------------------------
file results /'outputs/2006.txt'/ ;
*the .txt file
put results;
put '   the main results ' //
//
'I)Overall investment cost is' cost.l 'b€' //
//
'II)the Renewable capacity ' //
'Offshore                'CAPA.l('offshore')'    GW'//
'onsore                  'CAPA.l('onshore')'     GW' //
'PV                      'CAPA.l('PV')'  GW'//
'run of river            'CAPA.l('river') 'GW' //
'lake                    'CAPA.l('lake') 'GW' //
'biogas                  'CAPA.l('biogas')' GW'//
'gas                     'CAPA.l('gas')' GW'//
'Battery                 'CAPA.l('battery')' GW'//
'PHS                     'CAPA.l('phs')' GW'//
'methanation             'CAPA.l('methanation')' GW'//
//
//
'III)Needed storage volume and storage charging power' //
'Battery Storage         'CAPACITY.l('battery')'       GWh' //
'PHS Storage             'CAPACITY.l('phs')'       GWh'//
'methane storage         'CAPACITY.l('methanation')'   GWh'//
'Battery Storage         'S.l('battery')'       GW' //
'PHS Storage             'S.l('phs')'       GW'//
'methane storage         'S.l('methanation')'   GW'//
//
'IV)Secondary reserve requirements'//
'lake                    'reserve('lake') 'GW'//
'gas                     'reserve('gas')  'GW'//
'Pumped Storage          'reserve('phs') 'GW'//
'Battery                 'reserve('battery') 'GW'//
//
'V)Overall yearly energy generation of each technology'//
'Offshore                'gene_tec('offshore')'    TWh'//
'onsore                  'gene_tec('onshore')'     TWh' //
'PV                      'gene_tec('PV')'  TWh'//
'run of river            'gene_tec('river') 'TWh' //
'lake                    'gene_tec('lake') 'TWh' //
'biogas                  'gene_tec('biogas')' TWh'//
'gas                     'gene_tec('gas')'  TWh'//
'battery                 'gene_tec('battery')' TWh'//
'phs                     'gene_tec('phs')'  TWh'//
'Methanation             'gene_tec('methanation')'TWh'//
//
'VI)more details'//
'LCOE for Offshore               ' lcoe('offshore')' €/MWh'//
'LCOE for Onshore                ' lcoe('onshore')' €/MWh'//
'LCOE for PV                     ' lcoe('pv')' €/MWh'//
'LCOE for Run-of-river           ' lcoe('river')' €/MWh'//
'LCOE for Lake                   ' lcoe('lake')' €/MWh'//
'LCOE for Biogas                 ' lcoe('biogas')' €/MWh'//
'LCOE for Gas                    ' lcoe('gas')' €/MWh'//
'LCOS for battery                ' lcos('battery')' €/MWh'//
'LCOS for pumped storage         ' lcos('phs')' €/MWh'//
'LCOS for methanation            ' lcos('methanation')' €/MWh'//
//
'Load Curtailment' lc '% and the storage loss ' str_loss '%'//
//
'LCOE of the system is 'lcoe_sys1' €/MWh or 'lcoe_sys2' €/MWh '//
//
'Average spot price of the system is ' marginal_cost' €/MWh'//
//
'capacity factors '//
'Offshore                 'cf('offshore')//
'onsore                   'cf('onshore')//
'PV                       'cf('pv') //
'run of river             'cf('river') //
'lake                     'cf('lake') //
'biogas                   'cf('biogas')//
'gas                      'cf('gas')//
//
//
;
*-------------------------------------------------------------------------------
file hourly_generation /'outputs/2006.csv' / ;
*the .csv file
parameter nSTORAGE(str,h);
nSTORAGE(str,h) = 0 - STORAGE.l(str,h);
put hourly_generation;
hourly_generation.pc=5;
put 'hour'; loop(tec, put tec.tl;); put 'demand' , 'Pump' , 'ElecStr' , 'CH4' , 'price'/ ;
loop (h,
put h.tl; loop(tec,  put gene.l(tec,h);) put demand(h), nSTORAGE('PHS',h), nSTORAGE('battery',h), nSTORAGE('methanation',h), spot_price(h)/
;);
*-------------------------------------------------------------------------------
*                                The End :D
*-------------------------------------------------------------------------------
