Multi-objective optimization of a power-to-hydrogen system for mobility via two-stage stochastic programming

A systematic way for the optimal design of renewable-based hydrogen refuelling stations in the presence of uncertainty in the hydrogen demand is presented. A two-stage stochastic programming approach is used to simultaneously minimize the total annual cost and the CO2 footprint due to the electricity generation sources. The first-stage (design) variables correspond to the sizing of the devices, while the second-stage (operation) variables correspond to the scheduling of the installed system that is affected by uncertainties. The demand of a fleet of fuel cell vehicles is synthesized by means of a Poisson distribution and different scenarios are generated by random sampling. We formulate our problem as a large-scale mixed-integer linear program and we rely on a two-level approximation scheme to keep the problem computationally tractable. A solely deterministic setting which does not take into account uncertainties leads to underestimated device sizes, resulting in a significant fraction of demand remaining unserved with a consequent loss in revenue. The multi-objective optimization produces a convex Pareto front, showing that a reduction in carbon footprint comes with increasing costs and thus diminishing profit.


Introduction
The transportation sector is reported to contribute for around 40% of CO 2 emissions in Switzerland, making it a good target for decarbonization strategies. Indeed, the road-based private transport sector in Switzerland has an estimated CO 2 reduction potential of up to 40% by 2040 [1]. This is due to the rise of alternatives to fossil-fuelled vehicles, such as Electric Vehicles (EVs) and Fuel Cell Vehicles (FCVs). While EVs represent a more mature technology, the proliferation of FCVs is highly correlated to the development of new infrastructure, namely a network of hydrogen refuelling stations. To support the decarbonization strategy (e.g. as presented in [1]) it is crucial that the hydrogen required to run FCVs is produced via a carbonneutral process, with water electrolysis being such a possibility. However, a proper evaluation of the CO 2 impact of electrolysis has to be traced back to the sources producing the electricity required for running the electrolyzer: it is evident that if the electricity used derives from fossilbased sources, the electrolysis cannot be considered carbon-neutral.
Previous works have mainly studied the optimal design of charging stations for electric vehicles (see [2] and references therein), while corresponding studies on hydrogen refuelling  [3], the authors investigate the optimal sizing of the electrolyzer in a refuelling station, while [4] focuses on a feasibility study of renewable hydrogen stations. To the best of our knowledge, no previous studies on optimal design of hydrogen refuelling stations have considered the presence of uncertainties affecting the system.
In this work, we propose a systematic way for approaching the optimal design of renewable electricity-based hydrogen refuelling stations by considering uncertainties in the hydrogen demand. We cast our problem in the framework of two-stage stochastic programming. The first-stage ("here-and-now") decisions represent the sizing of the devices considered. The secondstage ("wait-and-see") decisions represent the scheduling of the installed system which is affected by uncertainties and therefore the corresponding cost is evaluated as expectation over multiple demand scenarios. We claim two advantages with respect to similar works in literature. First, we include the presence of uncertainties for the hydrogen demand in our formulation. Second, we propose a multi-objective optimization framework to simultaneously minimize the economic costs and the CO 2 footprint of the system evaluated in terms of Global Warming Potential (GWP) per kilowatt of electricity consumed.
The paper is structured as follows. In Section 2, we introduce the models for the devices constituting the hydrogen refuelling station. In Section 3, we formulate the optimization framework and we propose a two-level approximation scheme to solve it efficiently. In Section 4, we present the results of the numerical study. Finally, in Section 5, we conclude with ideas for refinement of our approach.

System modelling
We consider a hydrogen refuelling station composed of a PEM electrolyzer, a storage tank and an array of photovoltaic panels. We model the dynamics of the system as, where the vectors x t , u t and ξ t comprise respectively all the states, inputs and uncertainties affecting the system at time t ∈ T . Throughout we use the notation q (with q being a generic quantity) to denote the collection of q t with t ∈ T . The states considered here model the energy stored in the hydrogen storage tank and the electrolyzer stack temperatures, while the efficiency of the electrolyzer is modelled as static piecewise affine map. The uncertainty ξ t corresponds to the hydrogen demand from a fleet of fuel cell vehicles. The matrices A, B, C are assumed to be known and of appropriate dimension. We consider bilinear operational constraints for the system devices (e.g. capacity limits) of the form where the vector s indicates the design variables (i.e. the capacity of the installed devices) with appropriate upper and lower bounds ϕ(s) ≤ 0. b indicates a vector of binary decision variables describing the status on/off of the corresponding device during the whole horizon. Matrices E, F and functions h and ϕ are assumed to be known. The electricity carrier gives rise to the following energy balancing constraint, where matrices H p , H r , H u , H n are used to model the power flows affecting the balancing node.
Here p t is the power purchased from the grid, r is the power produced by renewable sources (i.e. the local photovoltaic array) which is assumed to be perfectly known for the whole horizon and n is associated to the curtailment of excess renewable energy (alternatively, one could also sell back to the grid the excess in solar power, adding an additional revenue). The dependence of the economic costs, namely investment cost J inv , maintenance cost J maint and operational cost J oper , from design and operational variables is modelled as, where C I , W I , C M are matrices of appropriate dimension, c t is the electricity price at time t in T and a is the annuity factor calculated according to the standard definition considering an investment rate of 6% and a system lifetime of 10 years. T indicates the number of time steps considered in the optimization: we consider one year with hourly resolution, in other words 365 · 24 = 8760 time steps. All parameters in this section are taken from [7]. Revenues R from selling hydrogen are evaluated as, where p H 2 ,t is the selling price of hydrogen at time t, assumed constant at 15 CHF/kg H 2 .

Optimization framework
Based on the model derived in Section 2, a two-stage stochastic program [9] is formulated to simultaneously maximize the annual profit and minimize the carbon footprint by relying on the -constraint method.
The objective function comprises a deterministic first-stage cost J inv (s) + J maint (s) which only depends on the "here-and-now" design variables s and an expected second-stage cost E[J oper (u, b, x, p, n, ξ) − R(ξ)] depending on the operation variables u, b, x, p, n and on the disturbances ξ. The "sustainable" objective J CO 2 (u, b, x, p, n, ξ) = t∈T g T t p t (where g t indicates the "dirtiness" of electricity at time t) on the second stage is enforced via a constraint that should hold in expectation; a more "aggressive" emission policy could be enforced by constraining all scenarios to have CO 2 footprint lower than [6].
The stochastic program (8) is solved by randomly sampling a finite number of scenarios (namely, the realizations of the uncertain parameter ξ over the whole horizon) from an underlying probability distribution P ξ to approximate the expectation in the second-stage cost and in the -constraint. Note that in (8) (8) where s represents a scenario in the set of all considered scenarios S with probability of occurrence π s , leading to the demand ξ s .
Program (9) is typically intractable as the number of second-stage variables grows with T ·|S| where T = 8760 and |S| can be in the order of hundreds (see [8]). To reduce the computational burden, a two-level approximation scheme is introduced. The first approximation is based on representative days (see [7] for a complete explanation of the representative days approach): the entire year is represented by means of R representative days (with R N ) which are selected according to a k-mean clustering approach based on the deterministic input data (in this work, electricity price and solar irradiance).
The second approximation is based on the selection of a subset S r of S with |S r | |S|. We model the daily hydrogen demand according to a Poisson distribution under the assumption that each day is independent from the others. We generate M samples (i.e. daily patterns) of the hydrogen demand by randomly sampling from the underlying distribution and reduce them to m representative samples with m M , using reduction techniques (see [6] and references therein). The probability of each reduced sample m i is calculated as the number of samples in the cluster m i divided by total number of samples generated. The resulting daily samples m are then concatenated to form |S r | = m R scenarios (resulting in yearly patterns) by considering all sequences deriving from all possible combinations of the reduced samples along the representative days. The corresponding probability of each scenario is given by the product of the probabilities of all samples composing the scenario.
In this way, the number of second-stage variables is reduced to K · R · m R , which for an appropriate choice of the daily time interval K, the representative days R and the reduced scenarios m R represents a significant reduction. The logic behind the two-level approximation scheme of (9) is represented in Figure 1.

Numerical results
We consider solar irradiance data for the city of Dübendorf (CH) in 2019 and electricity price data from EPEX Spot Market for Switzerland in the same year. We synthesize the hydrogen demand data simulating a small-size fleet of fuel cell vehicles by means of a Poisson distribution. We select R = 4 (corresponding approximately to "winter", "spring", "summer" and "autumn" conditions) as compromise between computational complexity and model accuracy and we reduce the number of daily samples to m = 3 (corresponding to "high", "medium" and "low" hydrogen demand). The MILP in (9) is solved using Gurobi branch-and-bound algorithm to the default optimality gap on an Intel i7-10510U CPU with 16 GB RAM.
We initially solve (9) for the maximum annual profit only (effectively set = +∞). Once (9) is solved, additional scenarios are generated in order to evaluate the out-of-sample performance of the method and the corresponding expected economic indicators for the system (Figure 2). Optimization results for the single-objective program for the deterministic and the stochastic settings are compared in Table 1. Table 1 shows that the deterministic case, which assumes a "medium" demand, tends to underestimate the device size for the electrolyzer and the storage tank by 26%: this leads to an expected unserved demand of 47'469 kWh/y (≈ 10% of the overall demand) with a consequent economic loss of 21'191 CHF/y. On the contrary the stochastic case leads to an expected unserved demand of 7'510 kWh/y with an economic loss of only 3'353 CHF/y. In both settings the optimal photovoltaic panels size for minimum costs achievement is 0 m 2 , corresponding to not having this device installed: this is due to the cheap electricity price with respect to photovoltaic installation costs. The key performance indicators for the stochastic setting are annualized Return On Investment (ROI) of 6.8 %, Net Present Value (NPV) of 948'460 CHF and payback time of 2.8 years. The net cash flow over system lifetime is reported in Figure 3 under the assumption that all inputs, namely electricity price, hydrogen selling price and hydrogen demand distribution, remain constant throughout the system lifetime. A more realistic setup should, however, consider changes in these inputs (e.g. the demand of FCVs is expected to grow).
A Pareto front is constructed by solving the two single-objective optimization problems (respectively, maximize profit and minimize carbon footprint) and relying on the -constraint method for the intermediate points between the two extremes of the front.
The Pareto front ( Figure 4) presents a sharp slope at the very beginning corresponding to a change in the system design. For low GWP value, larger electrolyzer and photovoltaic panels are required with a consequent increase in the investment cost up to 556 kW, 500 m 2 and 296'346 CHF/y respectively for the extremum point corresponding to a "sustainable" single-objective optimization (point B). As we allow for more and more emissions, the device sizes diminish down to the results recovered in Table 1 for the economic single-objective optimization (point A). At around 1'600 GWP/y (point D), the system configuration and the total annual cost do not change significantly anymore. The convexity of the Pareto front favours the inclusion of "sustainable" objectives into the design of hydrogen refuelling stations: for instance, we can get a reduction of 68% in GWP/y by only increasing the total annual cost by 12% (point C).

Conclusions
This paper presents a two-stage stochastic programming approach for the optimal design of hydrogen refuelling stations considering uncertainties in the hydrogen demand. We rely on two approximation schemes, namely representative days and scenarios reduction. We show that the deterministic case tends to underestimate the sizing of the devices, leading to the incapability of serving a significant fraction of the demand in out-of-sample scenarios. On the contrary, the stochastic formulation provides the reassurance of an optimal design with better out-of-sample performances. The multi-objective optimization produces a convex Pareto front indicating that the introduction of the "sustainable" goal into the optimal design of hydrogen refuelling stations is still economically viable. As future work, we plan to extend our work in order to incorporate uncertainties also in the other inputs and to consider the optimal selection of the facility location as part of the first-stage objective.