A Stochastic Programming Approach for the Planning and Operation of a Power to Gas Energy Hub with Multiple Energy Recovery Pathways

There is a need for energy storage to improve the efficiency and effectiveness of energy distribution with the increasing penetration of renewable energy sources. Among the various energy storage technologies being developed, ‘power-to-gas’ is one such concept which has gained interest due to its ability to provide long term energy storage and recover the energy stored through different energy recovery pathways. Incorporation of such systems within the energy infrastructure requires analysis of the key factors influencing the operation of electrolyzers and hydrogen storage. This study focusses on assessing the benefits power-to-gas energy storage while accounting for uncertainty in the following three key parameters that could influence the operation of the energy system: (1) hourly electricity price; (2) the number of fuel cell vehicles serviced; and (3) the amount of hydrogen refueled. An hourly time index is adopted to analyze how the energy hub should operate under uncertainty. The results show that there is a potential economic benefit for the power-to-gas system if it is modeled using the two-stage stochastic programming approach in comparison to a deterministic optimization study. The power-to-gas system also offers environmental benefits both from the perspective of the producer and end user of hydrogen.


Introduction
The increase in the share of intermittent renewable energy generation systems within the power grid without adequate energy storage systems leads to increased events of power supply and demand mismatch.To address this issue, research on short and long term energy storage systems have increased over the past few years.The adoption of the power to gas energy storage concept in countries with an existing natural gas system can help in an easier transition to a greener energy economy.The power to gas concept proposes to convert excess electricity to hydrogen through the water electrolysis process.The produced hydrogen will then be injected into and distributed via the existing natural gas infrastructure.In other words, the natural gas pipeline and underground storage infrastructure will act as an energy storage and distribution system.Based on a comparative analysis between compressed air, pumped hydro, and power to gas energy storage systems carried out by de Boer et al. [1], it was determined that countries with an established natural gas infrastructure can benefit from utilizing it as an energy storage sink in order to deal with increased penetration of intermittent renewable generation systems.An additional benefit of power to gas energy storage system is their ability to offer different energy recovery pathways.Victorsson et al. [2] highlight the advantages of coupling the renewable energy technologies producing cleaner power with the transportation sector.Their work specifically looks at making an economic assessment of an independent refueling station that produces its own hydrogen.They highlight the importance of off-grid solutions that are cost-competitive and benchmark realistic costs of a small scale hydrogen production system.Walker et al., Götz et al.,  highlight the various energy recovery pathways of the power to gas concept.Walker et al.'s work [3], shows the potential use of hydrogen produced from surplus electricity in the transportation, industrial sector, and the natural gas sector as a fuel.Götz et al. and Rönsch et al. [4,6] show the potential of combining the hydrogen with carbon dioxide to produce renewable natural gas.Otto et al.'s [7] work also focusses on utilizing the energy recovery pathways that produce hydrogen and methane and put them to use in the steel industry.The case study is based on the German steel industry and highlights the versatility and CO 2 emission reduction potential of using the power to gas concept to couple the industry with renewable energy generators.
In this study, a stochastic programming approach has been taken to plan and operate a power to gas energy hub that can provide hydrogen as a fuel to fuel cell vehicles, and to natural gas end users.This approach is novel and more appropriate for the planning of such energy systems as the solution is more robust and realistic.The energy hub is also modeled to participate in the demand response market for auxiliary energy services as regulated by the Independent Electricity System Operator, as it has been shown in an earlier work by this research group [8] that the addition of provision of such high value services can enhance the economic viability of the power-to-gas energy storage concept.

Literature Review
As access of renewable energy sources increases, it will eventually have a trickledown effect when the power produced by the energy sources is utilized in different sectors of the energy economy.Of late, maximizing the effective utilization of renewables has led to them being utilized to produce hydrogen for the transportation sector.Juan et al. [9] review the challenges (both strategic and operational) and potential environmental benefits of implementing the use of hydrogen and battery powered electric vehicles in urban communities.Franzitta et al., and Blanco-Fernández et al. [10,11] highlight the benefits of harnessing wind and wave energy for the production of hydrogen.This hydrogen can then be utilized for fueling fuel cell vehicles.Liu et al. [12] carry out an economic analysis on the hydrogen economy infrastructure required to sustain the introduction of hydrogen fuel cell vehicles in Ontario.They consider three scenarios each for the years 2015, 2025 and 2050 that represents the percentage share of fuel cell vehicles among all new vehicles sold in Ontario.Each of the three scenarios have different rates at which fuel cell vehicles are introduced in the market.In order to account for uncertainty in their economic analysis, they carry out a sensitivity analysis on the electricity price, water price, efficiency of the electrolysis process and the lifetime of the hydrogen production plant.Their uncertainty/sensitivity analysis shows that variability in electricity price can significantly impact the hydrogen system cost.Change in efficiency of electrolyzers although does effect the system economics, however it is seen that its impact is lower in comparison to variability in electricity price.The sensitivity analysis on price of water and the project lifetime are deemed to have insignificant impact on system economics, according to Liu et al. [12].
Kim et al. [13] carry out an optimization study to determine the optimal hydrogen supply chain technologies for the year of 2044 in different regions of South Korea.The estimated hydrogen demands for each region is used as a base (average) value around which two more scenarios are developed to account for uncertainty in hydrogen demand.The two scenarios considered assume that the demand varies by −20% of the average values for a 'below average' scenario, and +20% of the average values for an 'above average' scenario.The third scenario comprises of the average values itself.The problem is solved as a two-stage stochastic programming problem.However it should be noted that one of the assumptions of their study is that the hydrogen demand scenarios for each region is time invariant.
The conclusions of their study show that steam methane reforming is the cost optimal way to produce hydrogen.However, they also mention that with future cost improvements using water electrolysis to produce hydrogen is a cleaner alternative.
Almansoori et al. [14] develop a long-term hydrogen supply chain model for Great Britain that accounts for uncertainty in hydrogen demand through a scenario tree approach.They develop a multi-stage stochastic MILP problem.The time horizon considered in their work includes three time periods (each of 6 years) ranging from 2005-2022.Each time period has 9 scenarios.Their modeling study analyzes the effect of uncertainty on three different liquid hydrogen supply chain configuration cases where the hydrogen transportation mode has been varied.In each of the three configurations the transition from 'distributed forecourt plants' in the initial stages to 'centralized plants' at the end of the time period is observed.The primary hydrogen production source is seen to be from steam methane reforming plants.Dayhim et al. [15] develop a hydrogen supply chain for the state of New Jersey by implementing the two-stage stochastic programming approach.The uncertainty in their study is associated with the market penetration rate of fuel cell vehicles in the state.A total of 10 scenarios have been considered.The market penetration of fuel cell vehicles has been varied between 5-100%.The computational intensity is reduced as only 4 time periods have been considered.Their results show that steam methane reforming will initially be adopted for the production of hydrogen and in future time periods water electrolysis in combination with steam methane reforming will be used as the market for fuel cell vehicles grow.
Nunes et al. develop a two-stage stochastic optimization problem to design and plan a liquid hydrogen supply chain for Great Britain [16].They account for the uncertainty in hydrogen demand by generating scenarios via a first order autoregressive model (which includes a random error element following a normal distribution).The modeling time horizon spans a period of 18 years which is split into three periods comprising of six years each.They determine that 15 scenarios is a good limit to account for the uncertainty in hydrogen demand for a supply-chain planning study, without compromising computational time of solving the optimization problem.In addition to this, another conclusion of their study shows that starting with small steam methane reforming units and moving towards units having higher production capacities seems to be the most economical solution.
Taljan et al. [17] develop a model to optimally operate and assess the economic viability of a hydrogen production, storage and fuel cell subsystem of fixed capacities in Ontario.The subsystem is modeled to provide transportation fuel and provide electricity back to the power grid when selling prices are favorable.Their work focusses on assessing the viability of running the hydrogen production subsystem solely on power available from a nuclear power plant and a wind farm in Ontario.The uncertain parameters include electricity price and wind farm output.Each time point in their study has 30 electricity price scenarios (generated from a normal distribution), and only one wind scenario.Their analysis shows that electricity from the power grid when utilized for providing hydrogen to the transportation sector rather than providing energy storage, is an economically viable option.They also conclude that the utilization of the oxygen produced via water electrolysis and the heat produced from fuel cells improves the economics of the energy system.
Hydrogen supply chain designs have been proposed by work done by Han et al. and de-León Almaraz et al. for South Korea and a Mid-Pyrenees region, respectively [18,19].Both studies take a multi-objective optimization approach to carry out a comparative analysis between a centralized and a de-centralized production system based on their cost effectiveness, associated safety risks, and their environmental impact.In Han et al.'s work, steam methane reforming technologies are the most cost effective for large scale centralized system.In the centralized system, storing and transporting liquid hydrogen in tankers is the cost optimal.However, they deem this design to have higher risk and environmental costs.For the de-centralized design, a multiple number of small scale water electrolysis units, combined with compressed hydrogen tank storage and using pipelines for distributing them is deemed to be cost optimal.In this case, the system has lower risk and environmental cost, however it has a lower cost effectiveness.The hydrogen demand in Han et al.'s work is considered to be constant on a daily basis throughout the modeling horizon of one year (2044).De-León Almaraz et al.'s [19] system accounts for the growth of the demand side over the course of 30 years (2020-2050) with a 10 year time step.Their work gives preference to the use of water electrolysis for hydrogen production due to the growth of the renewable energy generation sector in France.The results presented in their study show the trade-offs for both running mono-and multi-period multi-objective optimization problems.
The stochastic analysis presented in this study tries to address the uncertainty in dynamic (hourly) electricity pricing by developing 5 possible electricity pricing scenarios for each hour over the course of a year.This is different from Liu et al.'s assumption of three different static electricity pricing scenarios for the entire period of 2015-2050.The uncertainty in hydrogen demand accounted by Kim et al. [13] is also static whereas in this study uncertainty in hourly (dynamic) hydrogen demand has been conside with respect to time in their red providing for a unique result.
There are two major points of difference when the work presented in this study is compared to work done by Almansoori et al., Dayhim et al.,.Firstly, this study focusses more on the effect of uncertainty in electricity pricing, and hydrogen demand on the sizing and dynamic (hourly) operation of an energy hub.Secondly, the two-stage stochastic optimization problem developed in this study is more focused towards developing an energy hub system that can satisfy gaseous hydrogen demand for a local region within the province of Ontario, and not the entire country or an entire state.Due to a finer time resolution of the optimization study used in this work the problem is more computationally intensive.Therefore, the number of scenarios considered here are lower in comparison to Almansoori  In comparison to Taljan et al.'s work where electricity pricing scenarios were obtained from a normal distribution, this study carries out a more rigorous probability distribution fitting analysis for electricity pricing.Different probability distribution functions have been fitted to historical electricity pricing for a period of 10 years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) in Ontario (Section 2.1).Taljan et al. and Han et al. [17,18] do not consider uncertainty and variability in hydrogen demand with respect to time in their respective studies.The uncertainty in hydrogen demand on an hourly basis has been accounted for through the development of five potential demand scenarios each hour (Section 2.2).In addition to providing hydrogen to the transportation sector, this modeling study also looks into the potential of hydrogen to be supplied to natural gas end users as hydrogen enriched natural gas and for the electrolyzers to provide ancillary services to the power grid.

Method
Techno-economic optimization models of power to gas energy hubs involve the use of data signals or data sets that can be categorized under the term 'random'.From wind speed, to electricity demand, there exist many such data sets that have an associated randomness with them.Such data sets/signals could have 'n' number of potential values at any given instant.In other words, there can be 'n' number of scenarios.Deterministic techno-economic optimization approaches look at solving the problem by considering expected values of such data signals (or sets).However, by taking the expected value, the accuracy of the physical system can be compromised.Another approach could be to solve numerous deterministic optimization problems for each of the possible scenarios of the random data sets [20].The results obtained could then be used to draw inferences from the results from each of the deterministic runs.However, this approach can be time consuming.Considering the randomness associated with these data sets is important as they can influence the operation of the energy hubs.Therefore, developing a stochastic optimization problem could be an alternative.
The study focusses on developing a power to gas energy hub that is able to provide three primary services, namely: (1) providing hydrogen as a clean alternative fuel for an initial influx of fuel cell vehicles; (2) injection of hydrogen in to the natural gas grid for potential energy arbitrage service; and (3) providing hourly demand response service to the power grid in Ontario.
The randomness associated with hourly electricity price, and hourly hydrogen demand data has been considered for the two stage stochastic optimization problem.The development of probability distribution functions with respect to these two data sets was chosen due to the availability of historical data for both of them.

Stochastic Hourly Ontario Electricity Price Data
In this study, the historical electricity price data in Ontario has been used in order to develop potential electricity prices that can occur in an hour.Information for this has been gathered from the Independent Electricity System Operator's (IESO) historical data archives [21].
Historical Hourly Ontario Electricity Price for the period of 2003-2013 has been considered for this study.The gathered data for each of the years have been first reoriented from an 8760 (h) × 1 (year) or 8784 (h) × 1 (year) matrix to a 365 (days) × 24 (h) or 366 (days) × 24 (h) matrix.This is done so that one can better analyze the trends in how the electricity price for each hour of the 365 days in a year looks like.However, before starting to carry out the fit test, the entire year's data is segregated into 4 seasons, namely winter (91 days), spring (91 days), summer (91 days), and fall (92 days).
After segregating each of the years (2003-2013) in the order described above, the data for winter, spring, summer and fall seasons were aggregated and arranged into 4 separate data sets with each having values from the period of 2003-2013 for the respective season.
In order to describe how the fitting tests are carried out, the fitting test on the winter season has been described.The fitting tests were carried out using the EasyFIT XL software (MathWave Technologies, Dnepropetrovsk, Ukraine) which is an add-on in Microsoft Excel.The software is used for each of the 24 h.To begin with, we select the column having hour 1's data and fit the distribution and check the fit rank (based on the chi squared test).The fit rank helps to determine the best probability distribution that fits the histogram plot of the column data.Upon fitting the different distributions available in the EasyFIT XL software to each of the 24 h in the accumulated data for every the winter season from 2003-2013, it is seen that the log-logistic (3P) distribution is the best fit.Similarly, after carrying out the same steps for the other three seasons, it is seen that log-logistic (3P) is the best fit.96 (24 h × 4 seasons) fit tests have been carried out to develop the stochastic data for the hourly Ontario electricity price.In order to give the reader an idea of the values that the fit parameters can take on, Table 1 lists values obtained for hour '1' from each of the four seasons.The symbols α, β, and γ are called the shape, scale and location parameters, respectively.Once the best fit distribution is decided, we use the gathered fit parameter for each of the 24 h (in each season) to generate random numbers that are basically the potential values that the electricity price can take for that specific hour.Now, once the potential values for each of the 24 h for every season have been produced, the new generated numbers can be combined and rearranged to develop a year worth of potential electricity prices.
Energies 2017, 10, 868 6 of 27 Note that this analysis could have been carried out with just one year worth of data, however having more data makes our case of suggesting the chosen best fit probability distribution more justifiable as it captures a greater number of possible values that the electricity price can take in a particular hour.

Stochastic Hourly Hydrogen Demand Data
In this study the power to gas energy hub proposes to meet hourly demand at a hydrogen refueling station located in the greater Toronto area (this is a large urban area with about 6 million people).The first step in preparing the stochastic hydrogen demand data includes setting an estimate on potential fuel cell vehicle market penetration.Oak Ridge National Laboratory carry out a cost and benefit analysis of ramping up fuel cell vehicle market penetration in the US [22].They propose 3 market penetration scenarios, where the number of hydrogen vehicles deployed starting from 2012 to 2025 have been reported.According to the report, scenarios 1 and 2 have been developed such that they follow a market penetration trend similar to what has been observed for the increase in the share of hybrid electric vehicles in North America in the past.The only difference between scenario 1 and 2 is that, scenario 2 assumed an earlier ramp up in market penetration beginning in 2015, whereas scenario 1 assumes a later ramp up in the year 2018.Scenario 3 is defined as an optimistic scenario, and it assumed that availability of refueling infrastructure will not be a constraining factor for the deployment of the fuel cell vehicles in the market.This study uses the market penetration scenario 2 projection for the year 2016.As the Oak Ridge National Laboratory Report presents scenarios for the US market, this study scales the figures down to a value for the greater Toronto area.This is done by initially converting the per capita fuel cell vehicle car ownership in the US to per capita fuel cell vehicle ownership in Canada (based on population data obtained via [23,24]).The per capita fuel cell vehicle ownership in Canada is then further scaled down to the per capita fuel cell vehicle ownership for the greater Toronto area [25].The total fuel cell vehicles in the greater Toronto area at the end of 2016 is calculated to be 1766.This value was used in this study to have an estimate of the size of the power to gas energy hub required to meet demands from a large fleet.Figure 1 above shows the planning region within which the 1766 fuel cell vehicles are going to be serviced.people).The first step in preparing the stochastic hydrogen demand data includes setting an estimate on potential fuel cell vehicle market penetration.Oak Ridge National Laboratory carry out a cost and benefit analysis of ramping up fuel cell vehicle market penetration in the US [22].They propose 3 market penetration scenarios, where the number of hydrogen vehicles deployed starting from 2012 to 2025 have been reported.According to the report, scenarios 1 and 2 have been developed such that they follow a market penetration trend similar to what has been observed for the increase in the share of hybrid electric vehicles in North America in the past.The only difference between scenario 1 and 2 is that, scenario 2 assumed an earlier ramp up in market penetration beginning in 2015, whereas scenario 1 assumes a later ramp up in the year 2018.Scenario 3 is defined as an optimistic scenario, and it assumed that availability of refueling infrastructure will not be a constraining factor for the deployment of the fuel cell vehicles in the market.This study uses the market penetration scenario 2 projection for the year 2016.As the Oak Ridge National Laboratory Report presents scenarios for the US market, this study scales the figures down to a value for the greater Toronto area.This is done by initially converting the per capita fuel cell vehicle car ownership in the US to per capita fuel cell vehicle ownership in Canada (based on population data obtained via [23,24]).The per capita fuel cell vehicle ownership in Canada is then further scaled down to the per capita fuel cell vehicle ownership for the greater Toronto area [25].The total fuel cell vehicles in the greater Toronto area at the end of 2016 is calculated to be 1766.This value was used in this study to have an estimate of the size of the power to gas energy hub required to meet demands from a large fleet.Figure 1 above shows the planning region within which the 1766 fuel cell vehicles are going to be serviced.It is important for a supply chain system to be able to account for the uncertainty in the demand that it is catering to.This study tries to account for the uncertainty in the hydrogen demand on an hourly basis over the course of a year by developing probability distribution functions of two essential factors to account for at a refueling station, namely: (1) The percentage of total number of refueling events over the course of a 24 h period at a refueling station; and (2) The different fueling amounts a station can encounter for a single refueling cycle.It is important for a supply chain system to be able to account for the uncertainty in the demand that it is catering to.This study tries to account for the uncertainty in the hydrogen demand on an hourly basis over the course of a year by developing probability distribution functions of two essential factors to account for at a refueling station, namely: (1) The percentage of total number of refueling events over the course of a 24 h period at a refueling station; and (2) The different fueling amounts a station can encounter for a single refueling cycle.
The data for the percentage of total number of refueling events for a 24 h period at a refueling station has been obtained from a presentation prepared by the National Renewable Energy Laboratory (NREL) on the "Performance Status of Hydrogen Stations and Fuel Cell Vehicles" [26]. Figure 2 represents the data gathered from the presentation.The data taken from the presentation has been segregated into two time frames: (1) 9 AM to 9 PM; and (2) 10 PM to 8 AM.In order to have a less complicated distribution function, the data was fit to two normal distribution functions.The distribution for the time period of 9 AM to 9 PM has a mean of 6.26% and standard deviation of 1.12, respectively.Similarly, for the time period of 10 PM to 8 AM, the distribution has a mean of 1.7% and a standard deviation of 0.9, respectively.Upon generating data from the two normal distributions, they have been combined to form data for a 24 h period for each of the days in a single year.Five different scenarios have been developed for the percentage of total number of refueling events occurring every hour.Note that the total number of refueling events is dependent on the total number of cars that has been estimated above to be 1766.
As mentioned earlier, the second essential factor considered in the study is the different fueling amounts in kg of H2 that a refueling station can encounter.The range of the data has been taken from NREL's data archive [27].To maintain the simplicity in the distribution considered, a normal distribution has been fitted to the data having a mean and standard deviation of 3.45 and 1.9, respectively.The amount of hydrogen refueled for a single refueling event ranges from 0.7-6.95kg.The final data was converted to kmol by dividing the scenarios developed by the molecular weight of H2 (2 kg per kmol).Figure 3 shows the normal distribution plot used for the uncertain parameter 'amount of hydrogen refueled'.The data taken from the presentation has been segregated into two time frames: (1) 9 AM to 9 PM; and (2) 10 PM to 8 AM.In order to have a less complicated distribution function, the data was fit to two normal distribution functions.The distribution for the time period of 9 AM to 9 PM has a mean of 6.26% and standard deviation of 1.12, respectively.Similarly, for the time period of 10 PM to 8 AM, the distribution has a mean of 1.7% and a standard deviation of 0.9, respectively.Upon generating data from the two normal distributions, they have been combined to form data for a 24 h period for each of the days in a single year.Five different scenarios have been developed for the percentage of total number of refueling events occurring every hour.Note that the total number of refueling events is dependent on the total number of cars that has been estimated above to be 1766.
As mentioned earlier, the second essential factor considered in the study is the different fueling amounts in kg of H 2 that a refueling station can encounter.The range of the data has been taken from NREL's data archive [27].To maintain the simplicity in the distribution considered, a normal distribution has been fitted to the data having a mean and standard deviation of 3.45 and 1.9, respectively.The amount of hydrogen refueled for a single refueling event ranges from 0.7-6.95kg.The final data was converted to kmol by dividing the scenarios developed by the molecular weight of H 2 (2 kg per kmol).Figure 3 shows the normal distribution plot used for the uncertain parameter 'amount of hydrogen refueled'.
As mentioned earlier, the second essential factor considered in the study is the different fueling amounts in kg of H2 that a refueling station can encounter.The range of the data has been taken from NREL's data archive [27].To maintain the simplicity in the distribution considered, a normal distribution has been fitted to the data having a mean and standard deviation of 3.45 and 1.9, respectively.The amount of hydrogen refueled for a single refueling event ranges from 0.7-6.95kg.The final data was converted to kmol by dividing the scenarios developed by the molecular weight of H2 (2 kg per kmol).Figure 3 shows the normal distribution plot used for the uncertain parameter 'amount of hydrogen refueled'.Note that the number of refueling events data set is time dependent.However, since the data for the amount refueled was only a range, it is assumed for the sake of simplicity to be time independent for each of its five scenarios developed for every hour in a year.

Two-Stage Stochastic Optimization Formulation
This sub-section describes the constraints and the objective function of the mixed integer stochastic linear programming problem.The optimization problem has been formulated in the general algebraic modeling system software (GAMS, Version 22.6 (GAMS Software GmbH, Frechen, Germany)).
The power to gas energy hub is modeled to provide primarily three services, namely: (1) meet pure hydrogen demand from fuel cell vehicles in the greater Toronto area; (2) inject hydrogen in to the natural gas distribution system to provide a cleaner burning fuel (hydrogen enriched natural gas) and at the same time also utilize the low electricity prices to achieve energy arbitrage between the power and natural gas grids; and (3) provide demand response service to the power grid as an ancillary service.Figure 4 gives an overview of the power to gas energy hub. Figure 4 below shows a block diagram layout of the power to gas energy hub modeled in this study.
Energies 2017, 10, 868 8 of 26 Note that the number of refueling events data set is time dependent.However, since the data for the amount refueled was only a range, it is assumed for the sake of simplicity to be time independent for each of its five scenarios developed for every hour in a year.

Two-Stage Stochastic Optimization Formulation
This sub-section describes the constraints and the objective function of the mixed integer stochastic linear programming problem.The optimization problem has been formulated in the general algebraic modeling system software (GAMS, Version 22.6 (GAMS Software GmbH, Frechen, Germany)).
The power to gas energy hub is modeled to provide primarily three services, namely: (1) meet pure hydrogen demand from fuel cell vehicles in the greater Toronto area; (2) inject hydrogen in to the natural gas distribution system to provide a cleaner burning fuel (hydrogen enriched natural gas) and at the same time also utilize the low electricity prices to achieve energy arbitrage between the power and natural gas grids; and (3) provide demand response service to the power grid as an ancillary service.Figure 4 gives an overview of the power to gas energy hub. Figure 4 below shows a block diagram layout of the power to gas energy hub modeled in this study.As mentioned earlier in Section 2.2, the fuel cell vehicle projection used in the study was for the year 2016.However, the data pertaining to the demand response service requirement from the power grid and natural gas flow in the distribution pipelines is only available for the period of November 2012-October 2013.Therefore, the modeling study assumes the November 2012-October 2013 timeframe and looks at the economics of a power to gas energy hub that can meet fuel demand from a high penetration of fuel cell vehicles (ahead of time in the market).
The two stage stochastic programming problem employs two types of decision variables, namely: (1) the first stage decision variables; and ( 2  As mentioned earlier in Section 2.2, the fuel cell vehicle projection used in the study was for the year 2016.However, the data pertaining to the demand response service requirement from the power grid and natural gas flow in the distribution pipelines is only available for the period of November 2012-October 2013.Therefore, the modeling study assumes the November 2012-October Energies 2017, 10, 868 9 of 27 2013 timeframe and looks at the economics of a power to gas energy hub that can meet fuel demand from a high penetration of fuel cell vehicles (ahead of time in the market).
The two stage stochastic programming problem employs two types of decision variables, namely: (1) the first stage decision variables; and (2) the second stage decision variables.The first stage decision variables are assigned values before having known the actual values of the random parameters.The uncertainty associated with the values of the random parameters requires at times a recourse action.The second stage variables fall into the recourse action category and have a particular cost associated with each of them.
The first stage decision variables in this study determine the size of both the electrolyzer system, and the compressed hydrogen storage system.The second stage decision variables determine how the power to gas energy hub operates.Some examples of the second stage variables include: amount of hydrogen produced, stored, withdrawn or purchased, amount of hydrogen injected in to the pipeline, and amount of load curtailment (demand response) offered.Appendix A lists the description of the symbols used for denoting the first and second stage decision variables in the stochastic optimization formulation.Appendix B lists all the parameters used in the study.
Equation ( 1) is used to estimate the hydrogen produced (G H 2 ,h,s ) by polymer electrolyte membrane (PEM) electrolyzers which are the chosen technology for this study.The electrolyzer efficiency is denoted by η El in kmol per kWh.The amount of energy consumed by the electrolyzers is denoted by E h,s , kWh.The subscripts 'h' and 's' in Equation ( 1) denote the hour and the subsequent scenario for the variables.These two subscripts hold the same meaning for the associated variables and parameters that will be described in the remainder of the section: The hydrogen produced by the electrolyzer can be broken down in to three streams (Equation ( 2), namely: (1) hydrogen injected into the on-site compressed storage system (I H 2 ,h,s , kmol); (2) hydrogen sent to the pressure reduction station for injection in to the natural gas distribution pipelines (PRS H 2 ,h,s , kmol); and (3) the stream that bypasses storage in the on-site tanks and is sent directly to the fuel cell vehicle refueling station dispenser (B H 2,h,s , kmol): The product of number of cars (N Car ), with each of the five scenarios for both stochastic parameters, namely, the percentage of total number of refueling events (RE H 2 ,h,s ), and the amount fueled (RA H 2 ,h,s ) gives five scenarios for the hourly hydrogen demand: Equation ( 3) shows the flow balance at the hydrogen dispensing point.A combination of hydrogen coming directly from: (1) the electrolyzers (storage bypass stream,B H 2,h,s , kmol); (2) the storage tank (O H 2 ,h,s , kmol); and (3) industrial producer (hydrogen purchased,P H 2 ,h,s , kmol) can be used to meet hydrogen demand.
The rated energy consumption capacity of the PEM electrolyzers used in the study is 1000 kWh (E Max ).Due to the fast ramping capabilities of the electrolyzers it is assumed that they can be turned off and restarted fairly quickly.Therefore, the lower operating limit (E Min , kWh) is set at 0 kWh [28].The constraint shown in Equation ( 4) is used to determine the number of PEM electrolyzers to be installed at the energy hub (N Ele ): The energy consumed by the electrolyzers is one of the second stage variables of the stochastic optimization problem.It can be regulated according the change in electricity price, the demand response required and the hydrogen demand placed on the energy hub.Equation (5) shows two variables that are used in order to denote the reduction in energy consumption by the electrolyzers.During hours when demand response is to be offered, the variable DR h,s (kWh) is used to denote the reduction in energy consumption.Consequently, the variable ER h,s (kWh) is used to denote the amount by which energy consumption is reduced when there is no demand response service that needs to be offered to the power grid.These two variables are subtracted from the maximum energy consumption limit of the electrolyzer system, denoted by the product of the number of electrolyzers (N Ele ) in the energy hub and the maximum energy consumption rating (E Max , kWh) of a single electrolyzer unit: The usage of the two variables ER h,s and DR h,s to denote energy consumption reduction will also require one to create a clear distinction between the hours in which demand response needs to be offered and the hours in which it does not have to be offered.The binary parameter α h (Equation ( 6)) is used to denote hours when demand response is required and it takes a value of 1 during these hours.There is no subscript 's' associated with α as demand response is not considered to be a stochastic random parameter in this study.The process of determining the hours in which demand response needs to be offered has been described in detail in work done by Mukherjee et al. [8].Equation ( 6) is used to set an upper and lower bound on the demand response capacity offered by the energy hub for a particular hour and scenario.According to the IESO, when participating in the demand response market, a load has to provide a minimum of 1000 kWh when it has been contracted to provide the service [29].The maximum reduction in load that the system can offer is its rated capacity.In other words the electrolyzers can be completely turned off.This upper bound is shown by the right hand side of the inequality constraint in Equation ( 6): Similarly, during the hours when no demand response is required (denoted by the term 1 − α h in Equation ( 7)), the variable ER h,s (kWh) can have a maximum value equivalent to the rated capacity of the electrolyzer system: The demand response service provider is issued a clawback cost by the IESO in hours when the entire contracted demand response capacity is not offered.The unit clawback cost ($ per kWh) is the same as the unit price offered to the load to provide demand response.The notation Y DR ($ per kWh) is defined to denote this unit price [30].The contracted capacity of the power to gas energy hub is taken to be equal to the rated capacity (N Ele × E Max , kWh) of the electrolyzer system as shown in Equation (8).The clawback (X h,s , $) is estimated by multiplying the difference between the rated capacity and the demand response offered with the demand response incentive (Y DR , $ per kWh): Equation ( 9) is used in the model to constrain the amount of hydrogen injected into the natural gas pressure reduction station (PRS H 2 ,h,s , kmol per hour).It is a rearrangement of the equation used to calculate the concentration of hydrogen in the natural gas system, where the denominator (PRS NG,h,s + PRS H 2 ,h,s , kmol) is taken to the right hand side of the inequality.At a particular instant, the concentration of hydrogen within the natural gas system is constrained to a maximum of 5 mole%.This limit is set by the notation θ (=0.05) in Equation ( 9).Melaina et al. define 5 mole% as a safe injection limit into natural gas distribution system in North America [31]: Energies 2017, 10, 868 11 of 27 Equation ( 10) is a simple energy balance equation that equates the energy content of the hydrogen and natural gas flow through the distribution pipeline to the hourly energy demand from the natural gas end user (ED h , MMBtu per hour).The energy demand from the natural gas end user is not a stochastic random parameter and therefore does not have the subscript 's' associated with it.The terms HHV H 2 and HHV NG are the high heating values of hydrogen and natural gas (in MMBtu per kmol): The power to gas energy hub sizes an on-site tank storage system used for the fuel cell vehicle refueling station.Equation (11) has been used to constrain the maximum hydrogen storage system capacity denoted by the variable St max (kmol).Since the goal of the optimization is to minimize net cost (Equation (20)), one can set a large enough arbitrary number as the upper bound for St max : The optimization solver then decides the value that St max should take such that it satisfies the inequality constraints shown in Equation ( 12).Since the problem optimizes net cost, the value of St max will most likely going to be a multiple of the maximum capacity of a single tank unit (C Max ).In this study, tanks with a maximum unit capacity of 45.4 kmol (or ~90 kg) at a maximum storage pressure of 172 bar have been used [32]: A minimum inventory level inside each tank should be maintained so that the tank pressure does not fall below 70 bar.The lower inventory limit was calculated using the ideal gas equation with compressibility factor of hydrogen used to account for real gas behavior.The lower limit on the tank storage comes from the knowledge that hydrogen withdrawn from storage is sent to a booster compressor unit that can compress gas from 70 bar to 825 bar to complete a refueling cycle for a fuel cell vehicle that has maximum on-board tank storage pressure of 700 bar [32,33].
The hydrogen inventory (Inv H 2 ,h,s , kmol) in the tank storage system at the end of every hour for each scenario cannot exceed St max .This relation is portrayed by Equation (13) below: Equation ( 14) shows the tank storage system's inventory balance.The inventory is updated at the end of every hour based on the amount of hydrogen injected (I H 2 ,h,s , kmol per hour) and/or withdrawn (O H 2 ,h,s , kmol per hour) from storage: The hydrogen goes through a compressor before being injected into the tank storage system.In this study a compressor module with a maximum flow handling capacity (F Max,C ) of 21 kmol per hour (or 42 kg per hour) has been used.The reciprocating compressor can take an inlet pressure of 42 psig and has an outlet pressure of 4500 psig.Therefore the optimization problem is constrained such that the product of number of compressors (N C ) and its unit flow capacity (F Max,C , kmol per hour) is greater than or equal to the hydrogen flow directed to storage for injection (I H 2 ,h,s , kmol per hour) at any given hour and scenario: The hydrogen enriched natural gas used to satisfy end user energy demand offsets a certain fraction of the natural gas that would have been sent if no hydrogen is injected in to the distribution pipelines.Equation ( 16) calculates this offset (NG O f f set,h,s , kmol per hour) by subtracting the chosen values of natural gas flow through the distribution pipelines (PRS NG,h,s , kmol per hour) by the optimization solver from the ratio of natural gas energy demand and high heating value of natural gas (HHV NG , MMBtu per kmol): There are two source points of emissions from the energy hub: (1) emissions associated with electricity bought for hydrogen production; and (2) emissions associated with hydrogen purchased from the market.
The emissions associated with the electricity bought is accounted for by multiplying the gas produced by the electrolyzers (G H 2 ,h,s , kmol per hour) with the ratio of power grid emission factor in Ontario (EF Grid,h , kg CO 2 per kWh) and electrolyzer efficiency (η El , kmol per kWh).
The hydrogen purchased from the market is assumed to be produced from the steam methane reforming process.Therefore, the product of hydrogen purchased (P H 2 ,h,s , kmol per hour) and the emission factor of the steam methane reforming process (EF SMR , kg CO 2 per kmol H 2 ) [34,35] summed over all scenarios and hours gives in technicality, the CO 2 emissions purchased by the energy hub.Equation ( 17) adds these two source points of CO 2 emission to give the total annual emission incurred (A CO 2 , kg of CO 2 emission incurred): In Equation ( 18), the CO 2 emission offset associated with the use of hydrogen enriched natural gas (as described above) has been calculated by summing the product of natural gas usage offset at the end user (NG O f f set,h,s , kmol per hour) and the well-to-wheel emission factor of natural gas (EF NG , kg CO 2 per kmol NG) [35], over all hours and scenarios.
The second part of Equation ( 18) calculates the emissions deferred by not using the steam methane reforming process to produce hydrogen.This is estimated by summing the product of total hydrogen flow coming from the electrolyzers (G H 2 ,h,s , kmol per hour) and the emission factor of hydrogen produced by the steam methane reforming process (EF SMR , kg CO 2 per kmol of H 2 ).The CO 2 emissions offset from selling hydrogen enriched natural gas to the natural gas end user and deferring the use of the steam methane reforming process for producing hydrogen are summed to calculate the total CO 2 emission offsets possible on annual basis (S CO 2 , kg of CO 2 emissions offset): The difference between S CO 2 and A CO 2 give the net CO 2 emissions offset by the power to gas energy hub (Net CO 2 ,O f f set , kg of CO 2 ): The objective of the optimization study is to minimize the annual net cost (Net Cost, $) of operating the energy hub, and is calculated from the difference between the annual cost (C, $) and annual revenue (R, $) (Equation (20)): The different contributors to the annual cost variable 'C' have been shown in Equation ( 21) below: Some of the new symbols in Equation ( 21) that have not been described in earlier equations include:

•
C Ele : Annual amortized cost of a PEM electrolyzer unit [28].This cost includes the cost of replacement of stacks which are assumed to have a stack life of 5 years.The replacement cost is taken to be equal to 15% of the capital cost [36].This includes the cost of running compressors along the natural gas pipeline [38].
Note that all costs were updated to 2012 Canadian dollars and amortized for a project lifetime of 20 years and an interest rate of 8%.The parameters listed above have been listed in Appendix C.
The first three terms in Equation ( 21) denote the cost of installing electrolyzer, tank and compressor units.The following terms in sequential order denote:

•
The annual cost of buying and transmitting electricity to electrolyzers and compressors.

•
Annual cost of buying water for H 2 production • Annual cost of distributing H 2 through the natural gas distribution system.

•
Annual clawback cost associated with participating and failing to provide demand response in the ancillary service market.

•
Annual cost associated with the hydrogen purchased from a third party vendor.
Equation (22) shows the terms that are included the annual revenue: The terms shown in Equation ( 22) that have not been described earlier include: • SP H 2 : The unit selling price of hydrogen when hydrogen is sold to the refueling station ($ per kmol) [39].

•
MP NG,h : Market price of natural gas.In this case it is assumed to be the Henry Hub Spot Price ($ per MMBtu).Hydrogen injected into the distribution line is sold to the natural gas utility that distributes it to its end users on an energy basis at this price.
In sequential order (as they appear), Equation ( 22) includes: • Annual revenue from selling hydrogen to fuel cell vehicles.

•
Annual revenue from selling hydrogen to the natural gas utility on an energy value basis.

•
Annual revenue from providing the demand response service.

•
Annual revenue earned from offsetting CO 2 emissions at the end user and using a cleaner method in comparison to SMR for producing H 2 .
The hydrogen sold to the natural gas utility in this study is undervalued.Therefore the optimization problem tries to minimize the price differential between its production price and selling price to the natural gas grid.Equation ( 23) compares revenues earned from selling hydrogen to the natural gas utility at the Henry hub spot price with what can be earned if the gas is sold at the same value as SP H 2 , which is the selling price of hydrogen to fuel cell vehicles.This comparison yields the annual revenue loss (RL PRS,H 2 , $): The energy hub uses the natural gas system as an energy sink when utilizing low electricity prices (that are an indication of surplus generation) to produce hydrogen.Even though the cost of production based on electricity price is low in periods of low electricity price, the energy hub currently does not receive any additional incentive to account for the costs associated with its capital investment when it injects hydrogen into the natural gas system.Therefore there is a loss incurred in this pathway.In order to constrain this loss, Equations (24)-( 27) have been incorporated in to the optimization problem.Equation ( 24) is used to estimate the annual revenue earned from meeting fuel cell vehicle hydrogen demand (R FCV,H 2 ) has been estimated in Equation (24).In this study, the selling price (SP H 2 ) is set at $16 per kmol.The National Hydrogen Association estimates market price of hydrogen as a transportation fuel between $5-$10 per kg (or $10-$20 per kmol) [39]: The selling price of hydrogen set in this study can be marked up by $4 per kmol.Equation ( 25) is used to estimate the additional revenue (∆ low ) that can be earned when the current selling price (SP H 2 ) is increased by $1 per kmol: Similarly, Equation (24) has been used to estimate the additional revenue that can be earned when (SP H 2 ) is increased by $4 per kmol: The variables ∆ low and ∆ up are then used to constrain the variable RL PRS,H 2 in Equation ( 27).In other words, Equation ( 25) can be looked at as a way to transfer the losses from selling hydrogen to the natural gas utility at the Henry hub spot price to what the energy hub system pays by not selling hydrogen as a transportation fuel at a higher price: Equation ( 28) could be used to estimate the capacity factor (CF Ele ) of the electrolyzer system.The term µ is used to denote the product of number of hours in a year (8760 h) and the number of scenarios chosen for the stochastic parameters (five scenarios): Equation ( 29) can then be used to set a lower and upper bound on the capacity factor.In this study, the lower bound (CF Ele,min ) has been set as a modest 0.65.This value is assumed to be the minimum value that the electrolyzer system must have to make a business case, according to an industry partner's feedback.CF Ele,max is set as 1: However, implementing Equation ( 28) makes the optimization problem a mixed integer non-linear programming (MINLP) problem.A non-linear stochastic programming problem may prove to be highly computationally intensive.Therefore, this study reformulates Equation ( 28) such that the problem can be made a mixed integer linear programming problem.
The first part of the reformulation involves the use of a geometric series with six terms (with a = 1, r = 2, and n = 0, 1, 2,..., 5).It should be noted that the ratio of sum of energy consumed by the electrolyzer system for each of the five scenarios and the number of electrolyzers creates the non-linearity in Equation (28).Since the number of electrolyzers is an integer variable, it is easier to relate it with a geometric series.This geometric series is shown in Equation (30).Each of the six geometric series terms have been multiplied by binary variables (β 0 , .., β 5 ) so that the optimization problem can decide on what terms of the geometric series to pick such that their sum adds up to the value of N Ele : As mentioned earlier, the right hand side of Equation ( 28) is a nonlinear term.The inequality constraint shown in Equation (31) has been used as a reformulation of Equation (28).The variables N Ele and CF Ele have been combined to form a new continuous variable denoted by γ n , on the left hand side.As N Ele = ∑ 5 n=0 (β n × a × r n ) (Equation ( 30)), one could say that γ n is used to represent the expression CF Ele × ∑ 5 n=0 (β n × a × r n ), where n = 0, 1, 2, .., 5. 5 In order for the optimization solver to understand what γ n represents, Equations ( 32)-( 34) have been included in the problem formulation.Equation (32) constrains γ n to have an upper limit denoted by the product of capacity factor of the electrolyzer system and the nth term of the geometric series (where: n = 0,1,2,..,5).The lower bound on γ n is set as the product of the binary variable, β n , the nth term of the geometric series (parameter) and the parameter CF Ele,min .In this way, for the β n 's (on the left hand side of Equation ( 30)) that are greater than zero, the corresponding γ n is forced to be greater than zero: By adding an additional upper bound on γ n denoted by Equation ( 33), the capacity factor of the electrolyzer (CF Ele ) gets bound to not exceed the maximum value of 1, denoted by the parameter CF Ele,max in Equation ( 33): Although Equation ( 34) is redundant, it is used as an additional aid on top of Equation ( 33), for the optimization solver to determine the value for the capacity factor of the electrolyzer system: The power to gas energy hub acts as a supplier of hydrogen to a refueling station collocated with it.It is assumed that the power to gas energy hub is the sole provider of hydrogen to this refueling station, even if the energy hub has to purchase hydrogen to meet demand from the refueling station.Modeling of the components associated with dispenser, booster compressor and cascaded tank storage unit within the refueling station have been ignored to reduce the computational intensity required to solve the problem.
The goal of this study is to assess the benefit of accounting for the uncertainty associated in developing of power to gas systems supplying multiple services.The study does not focus on the different ways of addressing uncertainty and therefore does not include a robust optimization analysis.

Stochastic Programming Concepts: EVPI and VSS
The value of the solution obtained from a stochastic optimization problem can be determined with the help of two existing theories, namely: (1) the expected value of perfect information (EVPI); and (2) the value of stochastic solution (VSS) [20].
The EVPI is a measure of the cost one is willing to pay in order to attain accurate information of the random parameters (or data sets/signals).This can also be interpreted as the cost one incurs for using 'prediction techniques' [41].In order to estimate the EVPI, one needs to first define the terms RP solution and WSS.The solution from the two stage stochastic programming problem is called the recourse problem solution (RP solution).WSS is defined as the wait-and-see solution, and can be calculated by initially solving the optimization problem for each of the possible scenarios of the random parameters, one by one.In others, one needs to solve as many deterministic optimization problems as the number of scenarios for the random parameter one considers in their work.Once the deterministic solutions are obtained, the mean of solutions to the objective function obtained from every deterministic optimization problem gives the WSS.The difference between the RP solution and WSS is called the EVPI (Equation ( 35)): In order to estimate the VSS, one needs to first estimate the solution to the expected value (EV) problem.The solution to the optimization problem when it is solved by taking the expected value of the random parameter gives one the EV.The values of the first stage decision variables obtained for the expected value problem (EV) need to be then used and fixed in the two stage stochastic optimization problem.Upon fixing the first stage decision variables (based on the values obtained from the EV problem solution) and solving the RP, one can estimate the EEV (expected result of using the EV solution).The EEV helps in determining the solution to the second stage decision variables when the first stage decision variables have been fixed.Subsequently, the difference between EEV and the RP solution gives the VSS (Equation ( 36)): According to Birge and Louveaux [20], the VSS helps in determining whether it is beneficial to set the first stage decision variables in the stochastic optimization problem (based on the solution obtained from the EV problem).In other words, if the difference between EEV and RP is negative, it is less beneficial to account for uncertainty in the parametric data sets/signals and the expected value solution is good enough.

Results
This section highlights the results of the stochastic optimization problem.The section presents analyses on the effects of the uncertain parameters on the first and second stage decision variables.The values of EVPI and VSS have also been estimated.The section ends with highlighting the importance of accounting uncertainty rather than taking a deterministic approach.
Note that Figure 5 is used to highlight how the energy hub operates in a particular hour for the five probabilistic scenarios of both electricity price (denoted by the line) and hydrogen demand (denoted by the height of the bar graph plots):

•
Probabilistic scenario 1: High electricity price ($0.11 per kWh), and moderately high hydrogen demand (184.9 kmol).The values of EVPI and VSS have also been estimated.The section ends with highlighting the importance of accounting uncertainty rather than taking a deterministic approach.Note that Figure 5 is used to highlight how the energy hub operates in a particular hour for the five probabilistic scenarios of both electricity price (denoted by the line) and hydrogen demand (denoted by the height of the bar graph plots):
The hydrogen demand in a particular hour can be met by: (1) withdrawing from storage; (2) producing hydrogen from electrolyzer and directly sending it to a dispenser; and (3) purchasing hydrogen from the market.Figure 5 shows how the energy hub operates for a particular hour when the HOEP (Hourly Ontario Electricity Price) and the hydrogen demand are stochastic inputs to the optimization problem.Both the HOEP and the hourly hydrogen demand comprise five different scenarios.In the above figure, it is observed that for scenarios 1, 3, 4 and 5, the HOEP is higher in comparison to scenario 2. Therefore, the hydrogen demand is met by withdrawing hydrogen from storage in scenarios 1, 3, 4 and 5.When the HOEP decreases for scenario 2, the share hydrogen coming directly from the electrolyzer is seen to increase.Despite the HOEP being the lowest in scenario 2, the high hydrogen demand in scenario 2 leads the optimization problem to satisfy 25.9% of the demand by purchasing hydrogen from the market.The electrolyzers cover 57.8% of the total hydrogen demand in scenario 2. The remaining share (16.3%) of hydrogen demand in scenario 2 is satisfied by withdrawing hydrogen from on-site hydrogen storage tanks.Therefore, even if HOEP is low and the hydrogen demand is high for a particular scenario, the energy hub may have to purchase hydrogen from the market to satisfy demand.Through Figure 6, it is shown how the energy hub varies its operating regime to provide demand response while still meeting the stochastic hydrogen demand for a particular hour in the optimization time frame.The stochastic parameter HOEP is also plotted in Figure 6 to show its influence in the The hydrogen demand in a particular hour can be met by: (1) withdrawing from storage; (2) producing hydrogen from electrolyzer and directly sending it to a dispenser; and (3) purchasing hydrogen from the market.Figure 5 shows how the energy hub operates for a particular hour when the HOEP (Hourly Ontario Electricity Price) and the hydrogen demand are stochastic inputs to the optimization problem.Both the HOEP and the hourly hydrogen demand comprise five different scenarios.In the above figure, it is observed that for scenarios 1, 3, 4 and 5, the HOEP is higher in comparison to scenario 2. Therefore, the hydrogen demand is met by withdrawing hydrogen from storage in scenarios 1, 3, 4 and 5.When the HOEP decreases for scenario 2, the share hydrogen coming directly from the electrolyzer is seen to increase.Despite the HOEP being the lowest in scenario 2, the high hydrogen demand in scenario 2 leads the optimization problem to satisfy 25.9% of the demand by purchasing hydrogen from the market.The electrolyzers cover 57.8% of the total hydrogen demand in scenario 2. The remaining share (16.3%) of hydrogen demand in scenario 2 is satisfied by withdrawing hydrogen from on-site hydrogen storage tanks.Therefore, even if HOEP is low and the hydrogen demand is high for a particular scenario, the energy hub may have to purchase hydrogen from the market to satisfy demand.
Through Figure 6, it is shown how the energy hub varies its operating regime to provide demand response while still meeting the stochastic hydrogen demand for a particular hour in the optimization time frame.The stochastic parameter HOEP is also plotted in Figure 6 to show its influence in the energy hub operating regime.Figure 6a shows two similar trends to Figure 5: (1) Hydrogen is withdrawn from storage to satisfy FCV demand in case of high HOEP; and (2) The share of hydrogen demand satisfied directly from electrolyzer increases when HOEP decreases.
Energies 2017, 10, 868 18 of 26 energy hub operating regime.Figure 6a shows two similar trends to Figure 5: (1) Hydrogen is withdrawn from storage to satisfy FCV demand in case of high HOEP; and (2) The share of hydrogen demand satisfied directly from electrolyzer increases when HOEP decreases.At a particular hour, the maximum and minimum demand response capacity that can be offered by the electrolyzer system within the energy hub is 17,000 kWh and 1000 kWh.In Figure 6b, it is seen that the maximum demand response capacity is offered for scenarios 1, and 4.This trend can be easier to interpret when one looks at the high value of HOEP and low value of hydrogen flow from storage for scenario 1 in Figure 6a.Therefore when Figure 6 are compared for the case of scenario 1, it is seen that the electrolyzer is able to provide its full demand response capacity (Figure 6b and satisfy the low hydrogen demand by only needing to withdraw hydrogen from on-site storage (Figure 6a).
Another interesting inference can be made when one looks at a scenario where electricity price (HOEP) is low and hydrogen flow in the energy hub is high (implying that hydrogen demand is high, see scenario 3 in Figure 6a).In this case the lower HOEP and higher hydrogen flow in comparison to scenario 1 (in Figure 6a), leads the optimization problem to decide that running electrolyzers at a At a particular hour, the maximum and minimum demand response capacity that can be offered by the electrolyzer system within the energy hub is 17,000 kWh and 1000 kWh.In Figure 6b, it is seen that the maximum demand response capacity is offered for scenarios 1, and 4.This trend can be easier to interpret when one looks at the high value of HOEP and low value of hydrogen flow from storage for scenario 1 in Figure 6a.Therefore when Figure 6 are compared for the case of scenario 1, it is seen that the electrolyzer is able to provide its full demand response capacity (Figure 6b and satisfy the low hydrogen demand by only needing to withdraw hydrogen from on-site storage (Figure 6a).
Another interesting inference can be made when one looks at a scenario where electricity price (HOEP) is low and hydrogen flow in the energy hub is high (implying that hydrogen demand is high, see scenario 3 in Figure 6a).In this case the lower HOEP and higher hydrogen flow in comparison to scenario 1 (in Figure 6a), leads the optimization problem to decide that running electrolyzers at a higher operating level and paying the clawback cost of not providing the entire demand response capacity (see Figure 6b) is a more profitable decision.In other words, it can be concluded that the demand response offered by the electrolyzers is directly proportional to the value of HOEP and inversely proportional to the value of hydrogen demand at any particular hour and scenario.
The power to gas energy hub in this study injects hydrogen in to the natural gas distribution pipeline system at a pressure reduction station.Figure 7 have been plotted to observe how the two stochastic parameters (i.e., Hourly H 2 demand and Hourly Ontario Electricity Price) affect hydrogen concentration within the natural gas pipeline.Both figures represent values for the five scenarios of the stochastic parameters from the same particular hour in the optimization run's result.Figure 7a compares the second stage decision variable i.e., H 2 concentration with the stochastic hydrogen demand, while Figure 7b compares H 2 concentration with stochastic electricity price.The hydrogen concentration is seen to have an inverse relationship with both the hydrogen demand and electricity price.In Figure 7a, it can be seen that when hydrogen demand is high, especially for scenarios 3 and 4, hydrogen concentration is zero.Similarly, the high electricity prices in scenarios 3 and 4 also contribute to the hydrogen concentration being zero.It can also be seen that the electricity price (in Figure 7b) is low in scenarios 1 and 4, but the hydrogen concentration is higher in scenario 4 in comparison to scenario 1.This can be explained by looking at Figure 7a where it is seen that scenario 1 has a higher hydrogen demand in comparison to scenario 4. The higher hydrogen demand leads the system to inject less hydrogen into the pipeline for the case of scenario 1.
Energies 2017, 10, 868 19 of 26 concentration within the natural gas pipeline.Both figures represent values for the five scenarios of the stochastic parameters from the same particular hour in the optimization run's result.Figure 7a compares the second stage decision variable i.e., H2 concentration with the stochastic hydrogen demand, while Figure 7b compares H2 concentration with stochastic electricity price.The hydrogen concentration is seen to have an inverse relationship with both the hydrogen demand and electricity price.In Figure 7a, it can be seen that when hydrogen demand is high, especially for scenarios 3 and 4, hydrogen concentration is zero.Similarly, the high electricity prices in scenarios 3 and 4 also contribute to the hydrogen concentration being zero.It can also be seen that the electricity price (in Figure 7b) is low in scenarios 1 and 4, but the hydrogen concentration is higher in scenario 4 in comparison to scenario 1.This can be explained by looking at Figure 7a where it is seen that scenario 1 has a higher hydrogen demand in comparison to scenario 4. The higher hydrogen demand leads the system to inject less hydrogen into the pipeline for the case of scenario 1.
Electricity prices influence the second stage decision variables of hydrogen injection and withdrawal from the storage tank.On an annual average basis, it is seen that hydrogen injection in to tank is higher than hydrogen withdrawal from it when electricity prices in each of the five scenarios are around the $0.019 per kWh (σ: 0.02).Hydrogen withdrawal is greater than hydrogen injection in to the tank when the electricity price is seen to go above the $0.055 per kWh (on annual average basis, σ: 0.034) mark for all of the five scenarios of electricity price.
As shown by Equation ( 3), hydrogen demand at a particular hour can be met by the power to gas energy hub via a combination of hydrogen from tank storage and a direct stream from the electrolyzers to the end user.Therefore, the stochastic electricity price also influences the hydrogen directly sent to meet hydrogen demand.On annual average basis, when electricity prices in each scenarios are low and near the $0.021 per kWh (σ: 0.0198) mark, producing and sending more hydrogen directly to meet demand is more optimal than withdrawing from the tank storage system.Electricity prices influence the second stage decision variables of hydrogen injection and withdrawal from the storage tank.On an annual average basis, it is seen that hydrogen injection in to tank is higher than hydrogen withdrawal from it when electricity prices in each of the five scenarios are around the $0.019 per kWh (σ: 0.02).Hydrogen withdrawal is greater than hydrogen injection in to the tank when the electricity price is seen to go above the $0.055 per kWh (on annual average basis, σ: 0.034) mark for all of the five scenarios of electricity price.
As shown by Equation (3), hydrogen demand at a particular hour can be met by the power to gas energy hub via a combination of hydrogen from tank storage and a direct stream from the electrolyzers to the end user.Therefore, the stochastic electricity price also influences the hydrogen directly sent to meet hydrogen demand.On annual average basis, when electricity prices in each scenarios are low and near the $0.021 per kWh (σ: 0.0198) mark, producing and sending more hydrogen directly to meet demand is more optimal than withdrawing from the tank storage system.
Implementing the power to gas energy hub concept within an energy economy transitioning towards increased share of renewable energy generation reduces the emissions associated with producing hydrogen.The net CO 2,e emissions offsets achievable by the power to gas energy hub (See Equations ( 16)-( 19)) accounts for : (1) the CO 2,e emissions offset from utilizing the water electrolysis process to produce hydrogen in comparison to using the steam methane reforming technology; (2) The CO 2,e emissions incurred while purchasing hydrogen from the market (based on the knowledge that the hydrogen bought is produced via the steam methane reforming process); and (3) The lifecycle CO 2,e emissions (associated with using 100% natural gas) offset achieved by selling a cleaner blend of fuel to the natural gas end user (hydrogen-enriched natural gas ).It is estimated that on an annual basis the energy hub can offset 9470.6 tonnes of CO 2,e emissions.
From the solutions to the recourse problem (RP) and the EEV problem presented in Table 2 above, the VSS is estimated to be 104,276.673(VSS = EEV−RP).The positive VSS denotes that the accounting for uncertainty in the modeling of the energy hub is beneficial.The WSS for the objective function (Net Cost) has been calculated to be −9241480.626.Based on this, the EVPI has been calculated to be 57,211.014(EVPI = RP−WSS).The high value of EVPI shows that one must pay a high cost to obtain exact information on the values of the random parametric data sets/signals.Therefore adopting the usage of stochastic programming is beneficial and one can save 104,276.673 in cost.Based on this trend it is likely to see the saving increase for greater variances in the data set.However, the study does not go further in to validating this.The levelized cost of producing hydrogen for the fuel cell vehicle fleet is estimated to be $3.145 per kg.Therefore, selling hydrogen to a fuel cell vehicle fleet at a premium price of $8 per kg [39] is a profitable energy recovery pathway for the power to gas energy hub.A kg of H 2 is equal to ~1 gallon gasoline equivalent.Therefore, the $8 per kg H 2 premium price and the $3.145 per kg production price if converted to gasoline equivalents, translate to a gasoline market price of $2.12 per liter and $0.83 per liter.The gasoline price in the region of Toronto was near $1 per liter of gasoline in 2016 [42].The gasoline equivalent price of selling hydrogen ($2.12 per liter, above) is roughly twice at what gasoline is being sold in the present day.This price paid by the end user of hydrogen seems to be high and can be lowered further by either providing the refueling station/power to gas energy hub with tax rebates that allows them to maintain their profit and at the same time sell the hydrogen closer to their production price of $3.145 per kg or $0.83 per liter gasoline equivalent.Another option could be that the government provides price rebates to the end user for using a cleaner hydrogen fuel.
The cost associated with the purchase of the Toyota Mirai (a fuel cell vehicle) in the US market is $58,385 (CAD $ 78,105) [43].Ontario, has created incentives for early adopters of the electric vehicle technology, where consumers receive an incentive of $6000-$10,000, depending on the size of battery pack.The consumers receive an additional incentive of $3000, if the price of the electric vehicle ranges from $ 75,000-$150,000.Creating such incentives for fuel cell vehicles can definitely create a huge incentive for the use of this technology.The government can implement further non-fiscal incentives (similar to what has been observed in Europe for the electric vehicles [44]).Some of these non-fiscal incentives could include, no charges on express toll routes, free parking in municipal parking lots, access to bus lanes, and reduced ferry rates.
Table 3 represents the comparison of the solutions obtained from solving the optimization problem for each of the five different scenarios of the uncertain parameters (HOEP, and Hourly H 2 Demand) individually.Each of the five scenarios are treated as five individual problems and the solutions obtained are treated as deterministic solutions.From the figures presented in Table 3, it can be concluded that the recourse problem gives a better solution in comparison to most of the deterministic solutions.It is only for scenario 4, where the deterministic solution is better.In this case the system has had to produce less hydrogen.This indicates the importance of the stochastic programming approach taken to primarily use the energy system to produce and satisfy high hydrogen demands under uncertainty.The deterministic solution will be better by a small margin only in cases where the hydrogen demand is lower than the hydrogen demand met by the recourse problem.Favorably lower deterministic hourly electricity prices in comparison to uncertain hourly electricity prices that may be higher for the recourse problems will also influence the difference in deterministic and stochastic solutions.In other words, the recourse problem that accounts for different uncertainty in hydrogen demand and electricity has more reliability in comparison to a deterministic problem with deterministic parameter inputs.

Conclusions
The study proposes the use of the two-stage stochastic programming approach to plan and operate a power to gas energy hub located in an urban area.The goal is to highlight the benefit of accounting for uncertainty in parameters that influence the operation of the power to gas energy hub.The energy hub system is modeled to be able to provide three services, namely: (1) meeting hydrogen demand from a fuel cell vehicle refueling station; (2) selling hydrogen to the natural gas end users as a mixed hydrogen enriched natural gas (HENG) fuel (not more the 5% hydrogen); and (3) provide hourly demand response service by participating in the demand response market (in this case as operated in Ontario by the IESO. The hourly electricity price, amount refueled at the refueling station, and the number of refueling events have been taken as the random (uncertain) parameters.The two different stages of decisions taken by the stochastic optimization problem include: (1) First Stage: The size of the individual components of the power to gas energy hub (electrolyzer, tank, and compressor); and (2) Second Stage: Decisions taken to operate the components of the energy hub.
Due to the finer time index (hourly) adopted in the modeling study, the two stage stochastic programming problem only includes five scenarios of each of the three stochastic parameters (electricity price, refueling amount, number of refueling events).The refueling amount and the number of refueling amounts at a particular time point have been multiplied with an optimistic projection of fuel cell vehicle market penetration to estimate the hydrogen demand.
The stochastic parameters are seen to affect the operation of the power to gas energy hub.Higher electricity price scenarios for a particular hour leads to the hydrogen storage being utilized to meet the hydrogen demand rather than direct from the electrolysor.Consequently, when the electricity price for a particular scenario in an hour is low, the share of hydrogen coming directly from the electrolyzer increases.So in jurisdiction where electricity price is often low it is thus likely that the storage requirements would be lower, and the electrolyser availability higher.However, the hydrogen demand is also stochastic and varies for the five scenarios in an hour.Thus when a scenario has low electricity price but very high hydrogen demand, the energy hub will need to purchase hydrogen from the market as shown in scenario 2 of Figure 3.
The demand response service offered is also seen to be affected by the stochastic nature of electricity price and hydrogen demand.From the analysis presented in Section 3, it is observed that the demand response offered by the electrolyzers is directly proportional to the value of electricity price and inversely proportional to the value of hydrogen demand at any particular hour or scenario.The power to gas energy hub is also able to offset approximately 9470.6 tonnes of CO 2,e emissions by forgoing the use of the steam methane reforming technology for producing hydrogen and using hydrogen enriched natural gas (HENG) that offsets lifecycle CO 2,e emissions incurred in using 100% natural gas at the natural gas end user.
The benefit of accounting for uncertainty is highlighted through the two stochastic programming concepts called Value of Stochastic Solution (VSS) and Expected Value of Perfect Information (EVPI).A positive VSS shows that by accounting for uncertainty one can have an economic benefit of $104,276 over a deterministic solution.A positive EVPI of $57,211 highlights the cost one must pay to obtain exact information about the random parameters.Future work will focus on looking at different ways for accounting for uncertainty, like robust optimization and also assessing further how the incentive program can be developed for the hydrogen economy by taking the example of the existing electric vehicle incentive program in Ontario.
et al., Dayhim et al., and Nunes et al.'s work.Considering more scenarios leads to an increase in the number of variables and that in turn increases the convergence time to the solution.Also in this stochastic study it is assumed that the energy hub is co-located at the refueling station, therefore there is no exchange of hydrogen between production centers as shown in each of Almansoori et al., Dayhim et al., and Nunes et al.'s work, and distributed hydrogen generation is thus examined in this work.

Figure 1 .
Figure 1.Ontario Planning Region Map denoting the fuel cell vehicle service area covered by the power to gas energy system.

Figure 1 .
Figure 1.Ontario Planning Region Map denoting the fuel cell vehicle service area covered by the power to gas energy system.

Figure 2 .
Figure 2. Variation in number of refueling events over the course of a day.Data is presented as a percentage of total number of refueling events occurring in a day.

Figure 2 .
Figure 2. Variation in number of refueling events over the course of a day.Data is presented as a percentage of total number of refueling events occurring in a day.

Figure 3 .
Figure 3. Normal distribution plot for the uncertain parameter 'Amount Refueled'.

Figure 3 .
Figure 3. Normal distribution plot for the uncertain parameter 'Amount Refueled'.

Figure 4 .
Figure 4. Power to gas energy hub system block diagram.
) the second stage decision variables.The first stage decision variables are assigned values before having known the actual values of the random parameters.The uncertainty associated with the values of the random parameters requires at times a recourse action.The second stage variables fall into the recourse action category and have a particular

Figure 4 .
Figure 4. Power to gas energy hub system block diagram.

Figure 5 .
Figure 5. Variation in hydrogen stream flows used to meet hydrogen demand for different electricity price and hydrogen demand scenarios for a particular hour.The line denotes the value of the electricity price for each of its five probabilistic scenarios in a particular hour.The height of the bars in the plot resemble the hydrogen demand for each of its five scenarios in a particular hour.

Figure 5 .
Figure 5. Variation hydrogen stream flows used to meet hydrogen demand for different electricity price and hydrogen demand scenarios for a particular hour.The line denotes the value of the electricity price for each of its five probabilistic scenarios in a particular hour.The height of the bars in the plot resemble the hydrogen demand for each of its five scenarios in a particular hour.

Figure 6 .
Figure 6.The bars in (a) show the different options utilized to meet the five probabilistic hydrogen demand scenario in an hour; (b) shows the change in demand response capacity offered for different scenarios of electricity price and hydrogen demand in an hour.(Note: hydrogen demand is only shown in (a)).(a,b) denote the values of the five probabilistic electricity price scenarios through the lines with markers.

Figure 6 .
Figure 6.The bars in (a) show the different options utilized to meet the five probabilistic hydrogen demand scenario in an hour; (b) shows the change in demand response capacity offered for different scenarios of electricity price and hydrogen demand in an hour.(Note: hydrogen demand is only shown in (a)).(a,b) denote the values of the five probabilistic electricity price scenarios through the lines with markers.

Figure 7 .
Figure 7. Effect of stochastic hydrogen demand (a) and stochastic electricity price (b); for a particular hour) on hydrogen concentration within the natural gas pipeline.

Figure 7 .
Figure 7. Effect of stochastic hydrogen demand (a) and stochastic electricity price (b); for a particular hour) on hydrogen concentration within the natural gas pipeline.

Table 1 .
Values of fit parameters for hour '1' in each of the four seasons.

Table 2 .
Values of Objective Function and Decision Variables for the Recourse Problem (RP), Expected Value Problem and EEV Problem.

Table 3 .
Comparison of Objective Function and Decision Variables for the Recourse Problem (RP), and Deterministic Problems.
H 2 Hydrogen inflow to the tank (kmol) PRS H 2 ,h,s Hydrogen injected into pressure reduction station (kmol) PRS NG,h,s Natural gas flowing through pressure reduction station (kmol) NG O f f set Amount of natural gas offset at the pressure reduction station (kmol) A CO 2 Total CO 2,e emissions associated with production and purchase of hydrogen (kg) S CO 2 Total CO 2,e emissions curbed while substituting natural gas with hydrogen and from not using steam methane reforming to produce on-site hydrogen.(kg)NetCO 2 ,O f f set Net CO 2,e emissions offset (kg) X Clawback cost for not offering scheduled demand response ($ per kWh) O&M C 2 Operating and maintenance cost of booster compressor modules that includes electricity consumption and transmission charges ($) RL PRS,H 2 Annual revenue loss in selling hydrogen at natural gas spot price ($) R FCV,H 2 Annual revenue earned from meeting hydrogen demand ($)∆ lowThe additional annual revenue that can be earned when hydrogen as a transportation fuel is sold at $17 per kmol∆ upThe additional annual revenue that can be earned when hydrogen as a transportation fuel is sold at $20 per kmolCF EleAnnual average capacity factor of electrolyzers β 0 , . . ., β 3 , . . ., β 5 Binary variables γ 0 , . . ., γ 3 , . . ., γ 5 Not Sure how to define it: Product of geometric series of constant ratio 2 and capacity factor variable