#--------------------------------------------------------------------------------------------#
#THE SCENARIO TREE
#--------------------------------------------------------------------------------------------#
#Number of years between scenario tree levels
param years;
#Calculation horizon in years = economic lifetime of investments
param lifetime;
#Number of scenario tree levels (stages)
param levels:=lifetime/years;
#Number of levels at which investments can be made
param inv_levels;
#Number of scenario tree building blocks
param blocks;
#Discount rate. (Not directly part of the scenario tree definition, but should be adjusted when scenario data is changed)
param rate;
#Set of all scenario tree building blocks
set BLOCKS:=1..blocks;
#Set of all scenarios (or paths through the scenario tree)
set SCEN;
#Set of all nodes in a scenario
set SCEN_NODES{SCEN} ordered;
#Set of all nodes
set TOT_NODES:=union{s in SCEN} SCEN_NODES[s] ordered;
#Defines the building block of each node
param block_no{n in TOT_NODES};
#The probability of each scenario to occur
param scen_prob{SCEN};
#The probability of each node
param prob{n in TOT_NODES}:=sum{s in SCEN: n in SCEN_NODES[s]} scen_prob[s];
#Auxiliary set used to express the parameter "parent"
set NODE_SCENS{n in TOT_NODES}:=setof{s in SCEN: n in SCEN_NODES[s]} s ordered;
#Shows which node is the parent of a node, that is, which branch it orginates from
param parent{n in TOT_NODES}:=if n == 0 then -1 else prev(n,SCEN_NODES[first(NODE_SCENS[n])]) symbolic;
#The level of a node
param level{n in TOT_NODES}:=if parent[n]==-1 then 0 else 1+ level[parent[n]];
#The root node
set ROOT within TOT_NODES := setof{n in TOT_NODES :parent[n]==-1} n;
#Nodes where investments can be made
set INV_NODES within TOT_NODES:=setof{n in TOT_NODES: level[n]<=inv_levels-1 && prob[n]>0} n;
#All nodes with probability >0
set NODES within TOT_NODES:=setof{n in TOT_NODES:level[n]<=levels && prob[n]>0} n ;
#--------------------------------------------------------------------------------------------#
#THE MODEL OF THE MILL
#--------------------------------------------------------------------------------------------#
#Definitions of sets
set USE ordered; #Technologies using a heat surplus at the mill to enable energy exports
set PRESSURE ordered; #Steam pressure levels (Ordered from high to low)
set MEASURE; #Steam saving measures (fixed-cost)
set REDUC; #Technology options causing pressure reductions (e.g. turbines, pressure-drop valves)
set PRODUCTS; #Energy "products" for export
set VARIED within USE; #Export options for which a finer time resolution is needed, e.g. all options affected by a varying
#district heating demand. (Includes, in this case study, all options using low-pressure steam)
set UNVARIED within USE:=USE diff VARIED; #Elements of USE not in VARIED
#--------------------------------------------------------------------------------------------#
#PARAMETERS
#...........Economic parameters..........#
#Prices, etc...
param product{USE} symbolic in PRODUCTS; #Export product of each "USE" option
param price{PRODUCTS,BLOCKS}; #Price of product [M€/GWh]
param operating_cost{USE,BLOCKS}; #Operating/running cost [M€/GWh]
#Preset value factor, used to calculated the present value of an investment cost
param pvf{l in 0..levels}:= (1+rate)^-(l*years);
#Present value sum, used to discount the value of a repeated cash flow to the beginning of a period.
param pvf_sum:=sum{i in 1..years} (1+rate)^-i;
#Residual value factor for investments with remaing economic lifetime at the end of the scenario tree.
param residual{l in 0..levels}:=((1+rate)^(l*years)-1)/((1+rate)^lifetime-1);
#Investment cost for fixed size/cost measures
param measure_cost{MEASURE}; #Investment cost [M€]
#Parameters for calculations of linearized investment cost data for energy export options
param ints{u in USE}; #Number of linearization intervals
param int_size{u in USE,1..ints[u]}; #Linearization interval sizes [MW],
param base_cost{u in USE}; #Base cost of investment [M€]
#Parameters a, b and c used to calculate the investment cost according to
# cost[M€]=aa * (size[MW] + cc)^bb
param aa{USE};
param bb{USE};
param cc{USE};
#Investment cost function
param inv_cost{u in USE,i in 0..ints[u]}:= if i=0 then base_cost[u] else aa[u]*(sum{j in 1..i} int_size[u,j]+cc[u])^bb[u];
#Slope of the investment cost function, i.e., cost per MW.
param slope{u in USE,i in 1..ints[u]}:=(inv_cost[u,i]-inv_cost[u,i-1])/int_size[u,i];
#........Steam balance parameters........#
#Steam savings at different steam pressure levels
param steam{MEASURE,PRESSURE}; #Steam savings [tonnes/h]
#States the conversion relation from steam flow to energy export
param conv{USE,PRESSURE}; #Conversion ratio: exports/steam [MW/(tonnes/h)]
#Factor used to include quench water in the steam balances. Quench water is mixed with the steam to saturate it when the steam pressure
#is reduced in e.g. a turbine or a pressure-drop valve.
param quench{REDUC,PRESSURE diff {first(PRESSURE)}}; #Steam flow out from turbine/valve (incl. quench water) / Steam flow in to turbine/valve []
#Parameter stating the maximum output which can be achieved without any steam use.
param max_free{USE}; #Max "steam-free" output [MW]
#Parameter stating the maximum decrease of steam flow to existing equipment, i.e., the current steam flow.
param max_loss{USE union REDUC, PRESSURE};#Steam flow decrease [tonnes/h]
#The current capacity of existing equipment
param existing_size{USE}; #Capacity [MW] NOTE1!
#Parameter to convert MW to GWh/year
param MW_to_GWh; #Unit conversion [GWh/year/MW]
#.....Yearly variations parameters........#
param per; #Number of periods per year for parameter/variables affected by yearly variations
param hours{1..per}; #Length of period [h].
param demand{1..per}; #Max output for options affected by yearly variations (typically district heating demand) [MW].
#.............CO2 parameters.............#
param co2{USE,BLOCKS}; #Reduction of CO2 emissions associated with energy export [ktonnes/GWh]
param co2_limit; #Req on expected "discounted" reduction (relative year 0) during the calculation horizon [ktonnes]
param co2_discount; #"Discount rate" for co2 emissions.
param co2_pvf{l in 0..levels}:= (1+co2_discount)^-(l*years); #Present value factor for emissions
param co2_pvf_sum:=sum{i in 1..years} (1+co2_discount)^-i; #Present value summation for emissions
#--------------------------------------------------------------------------------------------#
#VARIABLES
#Variables related to investment costs
var activate{MEASURE union USE, n in NODES} binary; #If =1, investment is made in the node
var deactivate{MEASURE union USE, n in NODES} binary; #Deactivation of investment
var newsize{u in USE,1..ints[u], n in NODES} >=0; #New capacity in linearization interval - investment is made in the node [MW]
var z{u in USE,0..ints[u]-1, n in NODES} binary; #Auxiliary variable for the linearization
#Variables used to define which investments have been made so far
var active{MEASURE union USE,n in NODES} >=0 <=1; #If =1, investment has been made in previous node
var size{USE,n in NODES} >=0; #Total available capacity [MW]
#Variables for the total output.
var output{USE, n in NODES diff ROOT}; #Yearly average energy output/export from technology option [MW]
var output_var{1..per, VARIED, NODES diff ROOT}>=0; #Output varying over the year [MW]
#Variables describing the steam surplus utilization for different export options.
var flow{USE union REDUC diff VARIED,PRESSURE,NODES diff ROOT}; #Steam flow used for technology option (yearly average) [tonnes/h]
var flow_var{VARIED,NODES diff ROOT,1..per}; #Varying steam flow (LP) for technology option [tonnes/h]
#Variable for the output achieved without any need for steam
var free{USE,NODES diff ROOT}>=0; #"Free" output [MW]
#--------------------------------------------------------------------------------------------#
#OBJECTIVE FUNCTIONS
#Primary objective function - Maximization of the expected value of the net present value of the investments [M€]
maximize fNPV: -sum{n in NODES} prob[n]*pvf[level[n]]*(1-residual[level[n]])*(sum{m in MEASURE} measure_cost[m]*
activate[m,n]+sum{u in USE} (activate[u,n]*base_cost[u]+sum{i in 1..ints[u]} slope[u,i]*newsize[u,i,n]))+
sum{n in NODES diff ROOT} prob[n]*pvf[level[n]-1]*pvf_sum*(sum{u in USE} (price[product[u],block_no[n]]-
operating_cost[u,block_no[n]])*MW_to_GWh*output[u,n]);
#--------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------#
#Alternative objective function - Maximization of the expected value of the CO2 emissions reductions [ktonnes] with a small adjustment for investment cost to eliminate "crazy" solutions
#~ maximize fCO2: sum{n in NODES diff ROOT} prob[n]*co2_pvf[level[n]-1]*co2_pvf_sum*(sum{u in USE} co2[u,block_no[n]]*
#~ MW_to_GWh*output[u,n])+0.1*(-sum{n in NODES} prob[n]*pvf[level[n]]*(1-residual[level[n]])*(sum{m in MEASURE}
#~ measure_cost[m]*activate[m,n]+sum{u in USE} (activate[u,n]*base_cost[u]+ sum{i in 1..ints[u]}
#~ slope[u,i]*newsize[u,i,n])));
#--------------------------------------------------------------------------------------------#
#GENERAL CONSTRAINTS
#Relates availability of measures and technologies to previously made investments.
subject to Activation{m in MEASURE union USE,n in NODES}: active[m,n]=if n in ROOT then 0 else
(active[m,parent[n]]+activate[m,parent[n]]-deactivate[m,parent[n]]);
subject to ActivationStop{m in MEASURE union USE, n in NODES diff INV_NODES}: activate[m,n]==0;
subject to DeactivationStop{m in MEASURE union USE, n in NODES diff INV_NODES}: deactivate[m,n]==0;
subject to Deactivation{u in USE, n in NODES}: active[u,n]>=deactivate[u,n];
subject to Reactivation{u in USE, n in NODES}: activate[u,n]>=deactivate[u,n];
#Linearization of investment costs
subject to Initialize{u in USE,n in NODES}: z[u,0,n] = activate[u,n];
subject to Linearization_Low{u in USE,i in 1..ints[u]-1,n in NODES}: z[u,i,n]*int_size[u,i] <=newsize[u,i,n];
subject to Linearization_High{u in USE,i in 1..ints[u],n in NODES}: newsize[u,i,n] <=z[u,i-1,n]*int_size[u,i];
#Total capacity
subject to InstalledSize{u in USE,n in NODES}:size[u,n]=if n in ROOT then 0 else
(size[u,parent[n]]+sum{i in 1..ints[u]} newsize[u,i,parent[n]] );
#Output limited by the steam flow
subject to SteamToOutput {u in UNVARIED, n in NODES diff ROOT}: output[u,n]<=sum{p in PRESSURE} flow[u,p,n]*conv[u,p]+free[u,n];
subject to SteamToOutput_var{u in VARIED,n in NODES diff ROOT, d in 1..per}: output_var[d,u,n]<=flow_var[u,n,d]*conv[u,'LP']+free[u,n];#NOTE1!
#"Free" output is limited
subject to FreeCapacity{u in USE, n in NODES diff ROOT}: free[u,n]<=max_free[u]*active[u,n];
#Steam balances on steam pressure levels (NOTE1)
subject to SteamBalance_HP {n in NODES diff ROOT}: sum{u in UNVARIED union REDUC} flow[u,'HP',n]<= sum{m in MEASURE} active[m,n]*steam[m,'HP'];
subject to SteamBalance_MP {n in NODES diff ROOT}: sum{u in UNVARIED diff REDUC} flow[u,'MP',n]<=
sum{m in MEASURE} active[m,n]*steam[m,'MP']+sum{u in REDUC} (flow[u,'HP',n]-flow[u,'MP',n])*quench[u,'MP'];
subject to SteamBalance_LP{n in NODES diff ROOT,d in 1..per}: sum{u in VARIED} flow_var[u,n,d]<=
sum{m in MEASURE} active[m,n]*steam[m,'LP']+(sum{u in REDUC} flow[u,'MP',n]*quench[u,'LP']);
#Steam flows has to be positive (except for net changes in existing equipment)
subject to PositiveFlow{u in USE union REDUC diff {VARIED}, p in PRESSURE,n in NODES diff ROOT}: flow[u,p,n]>=-max_loss[u,p];
subject to PositiveFlow_var{u in VARIED, n in NODES diff ROOT,d in 1..per}: flow_var[u,n,d]>=-max_loss[u,'LP']; #NOTE1
#Avoiding negative steam flows where the steam flow has no effect on the output.
subject to AvoidingNegative{u in UNVARIED, p in PRESSURE, n in NODES diff ROOT}: if conv[u,p]==0 then flow[u,p,n]=0;
#Calculates the average output for options with yearly variations
subject to AverageOutput{u in VARIED, n in NODES diff ROOT}: output[u,n]*sum{d in 1..per} hours[d]<=
sum{d in 1..per} output_var[d,u,n]*hours[d];
#Output is limited by invested capacity
subject to Capacity{u in UNVARIED, n in NODES diff ROOT}: output[u,n]<= size[u,n]- existing_size[u]*active[u,n];
subject to Capacity_var{u in VARIED, n in NODES diff ROOT,d in 1..per}: output_var[d,u,n] <=size[u,n]; #NOTE1!
#--------------------------------------------------------------------------------------------#
#CO2 CONSTRAINTS
#Reduced Co2 emissions [ktonnes]
subject to CO2Decrease: sum{n in NODES diff ROOT} prob[n]*co2_pvf[level[n]-1]*co2_pvf_sum*
(sum{u in USE} co2[u,block_no[n]]*MW_to_GWh*output[u,n])>=co2_limit;
#Reduced Co2 emissions in each scenario [ktonnes]
subject to CO2Scenario{s in SCEN}: sum{n in NODES diff ROOT: n in SCEN_NODES[s]}co2_pvf[level[n]-1]*co2_pvf_sum*
(sum{u in USE} co2[u,block_no[n]]*MW_to_GWh*output[u,n])>=0;
#--------------------------------------------------------------------------------------------#
#CASE-SPECIFIC CONSTRAINTS (and sets and variables)
set PIVAPS; #All process-integrated evaporation plants
set CONVAPS; #All modern conventional evaporation plants
set LIGVAPS; #Evaporation plants adapted for ligning extraction
set EVAPS; #All evaporations plants
set ANTI_PIVAP; #All measures that cannot be combined with process-integrated evaporation
set DH; #Alla district heating options
#Definition of specific variables
var ligvapcap{NODES}; #The evaporation plant capacity for ligning extraction [MW]
var newligvapcap{NODES}; #New capacity for lignin extraction (compare size/newsize) [MW]
#Specific requirements on combinations of measures
#Measures that cannot be combined with process-integrated evaporation
subject to Req1{m1 in ANTI_PIVAP, n in NODES diff ROOT}: active[m1,n]+sum{m2 in PIVAPS} active[m2,n]<=1;
#Measures that are required in order to have process-integrated evaporation
subject to Req2{n in NODES diff ROOT}:sum{m in PIVAPS} active[m,n]-active['HWWS',n]<=0;
#Only one evaporation plant can be "active"
subject to Req3{n in NODES diff ROOT}:sum{m1 in PIVAPS} active[m1,n]+sum{m2 in CONVAPS}active[m2,n]+active['ConVap',n]<=1;
#Investments in lignin extraction require an evaporation plant adapted for lignin extraction
subject to Req4{n in NODES diff ROOT}: active['LIG',n]-sum{m in LIGVAPS} active[m,n]<=0;
#The flash cannot be combined with DH100.
subject to Req5{n in NODES diff ROOT}: active['Flash',n]+active['DH100',n]<=1;
#The evaporation plant is designed for a certain level (or capacity) of lignin extraction.
#The capacity is the same as in the previous node with the possible addition of new capacity (which can be negative)
subject to Req6A{n in NODES}: ligvapcap[n]= if n in ROOT then 0 else ligvapcap[parent[n]]+newligvapcap[parent[n]];
#New capacity must equal zero if not any of the "LIGVAPS" are activated in which case newligvapcap is unconstrained
subject to Req6B{n in NODES}: newligvapcap[n]>=0 - 125*sum{m in LIGVAPS} activate[m,n];
subject to Req6C{n in NODES}: newligvapcap[n]<=0 + 125*sum{m in LIGVAPS} activate[m,n];
#The lignin output must equal the capacity of the evaporation plant unless the evaporation plant has been deactivated, in which case the output must equal zero. The latter case is
#handled by the general constraints above.
subject to Req6D{n in NODES diff ROOT}: output['LIG',n]>=ligvapcap[n]-125*(1-sum{m in LIGVAPS} active[m,n]);
subject to Req6E{n in NODES diff ROOT}: output['LIG',n]<=ligvapcap[n];
#A reduction of the turbine steam flow must be equally distributed on high-pressure and medium-pressure steam. This is not a requirement at the time of a new investment
subject to Req7A{n in NODES diff ROOT}: (flow['BP','HP',n]-flow['BP','MP',n])+
(if parent[n]<>0 then -(flow['BP','HP',parent[n]]-flow['BP','MP',parent[n]]))-activate['BP',parent[n]]*5000<=0;
subject to Req7B{n in NODES diff ROOT}: (flow['BP','HP',n]-flow['BP','MP',n])+
(if parent[n]<>0 then -(flow['BP','HP',parent[n]]-flow['BP','MP',parent[n]]))+activate['BP',parent[n]]*5000>=0;
#Capacity can only increase for existing turbine (NOTE2!)
subject to Req7C{n in NODES}: sum{i in 1..ints['BP']} newsize['BP',i,n]>=existing_size['BP']*(activate['BP',n]-deactivate['BP',n]);
#The following constraints concern the production increase
#Either the recovery boiler must be upgraded (RBU) or lignin has to be extracted
subject to Req8: activate['RBU',0]+activate['LIG',0]>=1;
#RBU can only be chosen in node 0
subject to Req9{n in NODES diff ROOT}: activate['RBU',n]==0;
#If RBU is not chosen, at least 53.6 MW of lignin has to be extracted
subject to Req10{n in NODES diff ROOT} :output['LIG',n]>=53.6*(1-active['RBU',n]);
#A new evaporation plant is required for the production increase
subject to Req11: sum{m in EVAPS} activate[m,0]==1;
#District heating deliveries are not allowed to decrease
subject to Req12A{n in NODES: n<>0 && parent[n]<>0}: sum{u in DH} output[u,n]>=
sum{u in DH} output[u,parent[n]];
subject to Req12B{n1 in NODES diff ROOT, n2 in NODES diff ROOT:parent[n1]==parent[n2]
&& n2<>n1}: sum{u in DH} output[u,n1]=sum{u in DH} output[u,n2];
#District heating piping is needed if any of the district heating options are chosen
subject to Req13{u in DH,n in NODES diff ROOT}: active['Piping',n]>=active[u,n];
#The district heating demand is limited by a demand curve
subject to Req14{n in NODES diff ROOT,d in 1..per}: sum{u in DH} output_var[d,u,n]<=demand[d];
#NOTE1: Not sufficiently generalized. Have to be adapted for other case studies
#NOTE2: Need to be generalized further