Optimized Configuration and Operating Plan for Hydrogen Refueling Station with On-Site Electrolytic Production

Abstract

Hydrogen refueling stations (HRSs) are critical for the popularity of hydrogen vehicles (fuel cell electric vehicles—FCEVs). However, due to high installation investment and operating costs, the proliferation of HRSs is difficult. This paper studies HRSs with on-site electrolytic production and hydrogen storage devices and proposes an optimization method to minimize the total costs including both installation investment and operating costs (OPT-ISL method). Moreover, to acquire the optimization constraints of hydrogen demand, this paper creatively develops a refueling behavior simulation method for different kinds of FCEVs and proposes a hydrogen-demand estimation model to forecast the demand with hourly intervals for HRS. The Jensen–Shannon divergence is applied to verify the accuracy of the hydrogen-demand estimation. The result: 0.029 is much smaller than that of the estimation method in reference. Based on the estimation results and peak-valley prices of electricity from the grid, a daily hydrogen generation plan is obtained, as well as the optimal capacities of electrolyzers and storage devices. As for the whole costs, compared with previous configuration methods that only consider investment costs or operating costs, the proposed OPT-ISL method has the least, 8.1 and 10.5% less, respectively. Moreover, the proposed OPT-ISL method shortens the break-even time for HRS from 11.1 years to 7.8 years, a decrease of 29.7%, so that the HRS could recover its costs in less time.

1. Introduction

Air pollution and the lack of fossil fuels have become severe in recent years. To solve the two problems, conventional internal combustion engine vehicles will be replaced by plug-in electric vehicles (EVs) [1] or hydrogen-powered vehicles (FCEVs) [2]. Compared with EVs, FCEVs can be more promising because of their superior features of fast refueling rate, high mileage range, and zero pollution [3].

Generally, hydrogen can be produced through electrolytic, thermal, and biochemical processes [4]. Among these technologies, electrolysis is expected to be the main generation method [5]. It is a mature technology of applying direct electric current to water to dissociate it; moreover, the green hydrogen obtained has a high purity reaching 99.999 vol.% [6].

Normally the producing process happens in electrolyzers and the produced hydrogen is stored in storage-equipped fueling stations to fulfill the demand for FCEVs [7]. For hydrogen refueling stations (HRSs), hydrogen can be generated on-site or at large central production plants and be transported to the stations [8]. The latter approach costs less for hydrogen production; however, the spend on transporting hydrogen can be excessive [9]. Ref. [10] develops a hydrogen supply chain planning model that determines the least-cost mix of H₂ generation, storage, transmission, and compression facilities to meet H₂ demand. A wide range of hydrogen-related technology options is studied.

By contrast, the on-site hydrogen stations eliminate the transportation cost and enhance the safety assurance for the hydrogen industry.

Even though the on-site hydrogen stations with electrolyzers are proving viable [8], their proliferation is limited by exorbitant construction costs and the high price of electricity for operation [6]. Therefore, researchers are studying various methods to realize cost reduction.

The efficiency of commercial electrolyzers is less than 60% [5]; thus a massive amount of electricity is needed to satisfy hydrogen demand [8]. Some researchers utilize cheap renewable energy, replacing traditional electricity from the grid to supply power to electrolyzers [11]. Most of them study the scheduling method during an HRS’s operation and concentrate on the stability and economy of the whole system. Ref. [12] analyzes the modeling and control of a hybrid-drive wind turbine with hydrogen energy storage and studies the regulation function of hydrogen in a wind power generation system. Ref. [13] proposes a cooperative operation model for the wind turbines and HRSs considering the individual benefit. A planning model for an electricity-hydrogen integrated energy system (EH-IES) is proposed in [14], with considerations of hydrogen production and storage technologies. Similarly, they both consider the risks deriving from electricity price uncertainties [13] and generation-load uncertainties [14]. Moreover, seasonal storage and interregional hydrogen supply chains (HSCs) are employed by some researchers to eliminate the imbalance between renewable energy and hydrogen demand [15].

Expect for cheap electricity from renewable energy, some researchers focus on the peak-valley difference of electricity price from the power grid, decreasing the operation cost by giving the electrolyzer different production scheduling hourly. In [7], the dramatically changeable electricity prices provide access for the station to participate in the power market with a hydrogen-storage tank. The operating reserve provision model is proposed in [8] to intensify the economic feasibility of the investment. Ref. [16] proposes a model for optimal scheduling of privately owned hydrogen storage stations by exploiting the lower electricity market prices to reduce the power purchase cost.

Even though investment costs have been considered in the optimization of some integrated energy systems [17], hydrogenation station scenarios with hydrogen energy demand estimations are an original study point. Ref. [18] proposes a model for optimal day-ahead scheduling of power-to-grid storage and studies the gas load management in electricity and gas markets to minimize the cost of gas consumption for the gas load. A gas demand forecasting algorithm is integrated into the optimal scheduling model using soft constraints, slack variables, and penalizing mechanisms. However, this forecasting algorithm only considers the information from data but ignores the behavior of hydrogen vehicles.

Installation cost accounts for over 30 percent of the total investment; therefore, it is supposed to be considered for cost reduction. To reduce operating costs, it is feasible to make proper use of the peak-valley difference in electricity prices [16].

To minimize the total cost during the HRS life cycle, this paper proposes a configuration method for installation and hourly generation scheduling for electrolyzers. As for installation investment, the optimal capacity for electrolyzers, as well as for storage devices, is investigated. The peak-valley difference of electricity prices is in full use to reduce electricity costs. The operation constraints of HRS have been mentioned in previous literature. However, precise demand constraints are needed to solve the optimization problem, which has not been studied yet. In this paper, a refueling behavior simulation method for FCEVs is proposed, contributing to hydrogen demand estimation for HRS, which is significant to solve the optimization problem.

The main contributions of this paper are summarized below:

• Mimicking the research into electric vehicles, the refueling behavior of FCEVs is simulated and a hydrogen demand estimation method is proposed. The estimation results provide relatively accurate demand constraints for the optimization model to minimize the total cost of HRS.
• An OPT-ISL method is developed for both the installation investment and operation cost. The optimal configuration of electrolyzers and storage devices is obtained with this method. The peak-valley difference of electricity prices from the power grid is fully utilized with hydrogen generation plan from the OPT-ISL method.
• The economical efficiency to reduce the total cost of HRS with the proposed OPT-ISL method is proven to be optimal.

The remainder of this paper is organized as follows. Section 2 proposes a simulation method for FCEVs’ refueling behavior and a hydrogen demand estimation method for HRS with a certain amount of FCEVs. An optimization formulation with objective function and constraints is presented in Section 3. Case studies and profit analysis are conducted in Section 4. Finally, conclusions are given in Section 5.

2. Hydrogen Demand Estimation

In the micro-grid with HRS, researchers generally use the previous hydrogen demand curve to predict the subsequent curve, since hydrogen demand is difficult to estimate in complex micro-grid scenarios [9]. In this study, as the HRS takes power directly from the grid and the fluctuation of electricity price and accurate estimation of hydrogen demand are used to minimize the cost, the research scenario is relatively simple. Therefore, it is feasible to realize an accurate hydrogen demand estimation.

As the output constraint of HRS, hydrogen demand is a critical factor to optimize the configuration of hydrogen refueling stations and calculate the benefits. This depends on the time, the location of the station, and other issues [17]. In previous studies, a typical hydrogen demand is usually assumed, with random fluctuations to make it more reliable [10]. However, the daily hydrogen demand follows certain rules, which is not a typical value and can be estimated in some ways.

This chapter proposes a more accurate method to estimate hydrogen demand if the vehicles served are decided.

2.1. Research Background

Compared with HRSs, the demand for gas stations and charging stations [19] are studied more often. The engineering analogy method performs well for gas stations to estimate hydrogen demand. In this method, the total amount of hydrogen refueling stations is estimated based on the number of urban motor vehicles; that is, the station serves a certain number of motor vehicles, which can match the supply and demand. Since HRSs and gas stations share many similarities, the engineering analogy method can possibly be applied to estimate the hydrogen demand.

It is assumed that the HRS operates the whole 24 h in a day and each FCEV refuels itself only during its driving time on the road. Fuel quantity for FCEVs depends on their fuel tank volumes. The refueling rate is related to the minimum amount of residual hydrogen acceptable to an FCEV [12].

Thus, the hydrogen demand per hour from each FCEV can be expressed as the product of possibility on the road, refueling rate, and fuel quantity, which can be estimated respectively to obtain hydrogen demand, since for individual FCEVs there are only two states: on the road or not.

There are three different types of FCEVs: taxis, private cars, and buses classified by their travel habits. We derive three driving patterns for each of them. On the basis of the Fermi Estimation Principle [20], this will produce a more reasonable estimation.

According to studies of plug-in electric vehicle (PEV) driving patterns by using statistics from the National Household Travel Survey in [21], the home arrival (office departure) times and the office arrival (home departure) times of plug-in electric vehicles can be simulated following the normal distribution. Their daily driving distances can be modeled by a logarithmic normal distribution. Since the lifestyles of the owners are hardly affected by whether the vehicle is FCEV or PEV, the driving pattern of private FCEVs can be considered similar.

2.2. Private FCEVs

Private FCEVs seem to leave home in the morning and arrive back home at night. Similarly, the pattern of PEVs can be applied to private FCEVs.

The home departure times and home arrival times both satisfy the normal distribution $N(\mu ,{\sigma }^2)$ as Equation (1), while the daily driving distances d can be modeled by a logarithmic normal distribution shown in Equation (2), assuming the only time for private hydrogen cars to refuel is when they are on road, from home to workplace, and backtracing. For individual private FCEVs, the one-hot encoding shown in Equation (3) is employed to show whether it is on-road or not:

$$f_{pri}^k(T)=\frac{1}{{\sigma }_{pri}^k\sqrt{2\pi }}\exp \left[-\frac{{(T-{\mu }_{pri}^k)}^2}{2{\sigma }_{pri}^k{}^2}\right]$$

(1)

where k = l, a to represent departure times or home arrival times.

$$f_{pri}(d)=\frac{1}{d{\sigma }_{pri}^d\sqrt{2\pi }}\exp \left[-\frac{{(\ln d-{\mu }_{pri}^d)}^2}{2{\sigma }_{pri}^d{}^2}\right]$$

(2)

$${\kappa }_{pri}^i(t)=\left\{\begin{array}{l}1,\hspace{1em}t=t_{pri}^{il}\\ 1,\hspace{1em}t=t_{pri}^{ia}\\ 0,\hspace{1em}t\in (0,24]\hspace{1em}and\hspace{1em}t\ne t_{pri}^{il}|t_{pri}^{ia}\end{array}\right.$$

(3)

Since there is no direct relationship between daily driving distance and departure time, home arrival time, the average daily distance $e^{{\mu }_{pri}^i}$ is applied to replace $d_{pri}^i$.

For private FCEVs, whether they would refuel themselves depends on the state of charge (SOC), which is normally 0 to 1. Each FCEV should ensure its SOC is above $SO{C}_{\min }^i$ when on the road, which is determined by the driving distance between HRS and home (office) $l_{pri}^i$, hydrogen consumption per kilometer $q_{pri}^i$, and nominal hydrogen capacity $C_{pri}^i$. However, obtaining $l_{pri}^i$ for each private FCEV is impractical; hence, $d_{pri}^i$ replaces $l_{pri}^i$ to describe $SO{C}_{pri\min }^i$ as Equation (4) shows:

$$SO{C}_{pri\min }^i={\alpha }^i+\frac{d_{pri}^i\cdot q_{pri}^i}{2\cdot C_{pri}^i}\times 100\%$$

(4)

where ${\alpha }^i$ is assumed as a constant value between 15 and 25%. Therefore, for private FCEV i, the refueling rate can be expressed as Equation (5):

$${\lambda }_{pri}^i=\frac{d_{pri}^i\cdot q_{pri}^i}{C_{pri}^i(1-SO{C}_{pri\min }^i)}$$

(5)

2.3. Non-Private FCEVs

There are roughly two kinds of non-private FCEVs: taxis and buses. Taxi drivers also have starting and leaving time of work, which satisfy the normal distribution. For taxis with hydrogen, the possibility on roads can be expressed as Equation (3), while the parameters for $f_{taxi}^k(T)$ are different. For individual taxis, the time on the road is also decided by home departure time and home arrival time, as shown in Equation (6):

$${\kappa }_{taxi}^i(t)=\left\{\begin{array}{l}1,\hspace{1em}t=t_{taxi}^{il}\\ 1,\hspace{1em}t=t_{taxi}^{ia}\\ 0,\hspace{1em}t\in (0,24]\hspace{1em}and\hspace{1em}t\ne t_{taxi}^{il}|t_{taxi}^{ia}\end{array}\right.$$

(6)

Taxis are assumed to leave home earlier and arrive back home later. Whereas, according to the study of taxis, the refueling time of them is mostly in the period of 11:00–14:00 and 20:00–24:00, which is different from private cars [22].

Therefore, for hydrogen taxis, the refueling times also obey normal distribution in these two periods, and the daily driving distances obey logarithmic normal distribution. This is because the paths of taxis commonly overlap on the road when taking customers from home to workplace. The SOC also decides whether the taxis refuel themselves, which satisfies normal distribution with a large mean value and little square deviation compared with private cars. This is because the expected travel distance for the next customer is uncertain and drivers may comply with specifications set by the company. $SO{C}_{taxi\min }^i$ and ${\lambda }_{\mathrm{taxi}}^i$ are shown in Equations (7) and (8). Their daily driving distances d can be modeled by a logarithmic normal distribution shown in Equation (9):

$$\begin{array}{l}f(SO{C}_{taxi\min }^i)=\frac{1}{{\sigma }_{taxi}^{soc}\sqrt{2\pi }}\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\times \exp \left[-\frac{{(SO{C}_{taxi\min }^i-{\overline{SOC}}_{taxi\min })}^2}{2{({\sigma }_{taxi}^{soc})}^2}\right]\end{array}$$

(7)

$${\lambda }_{taxi}^i=\frac{d_{taxi}^i\cdot q_{taxi}^i}{C_{taxi}^i(1-SO{C}_{taxi\min }^i)}$$

(8)

$$\begin{array}{l}{f}_{taxi}(d_{taxi}^i)=\frac{1}{d{\sigma }_{taxi}\sqrt{2\pi }}\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\times \exp \left[-\frac{{(\ln d_{taxi}^i-{\mu }_{taxi}^d)}^2}{2{\sigma }_{taxi}{}^2}\right]\end{array}$$

(9)

Hydrogen buses can be scheduled in terms of passenger flow and road congestion [23], which is not the point of this paper. Here we assume HRS refuels the buses just before they leave for work or arrive after work.

Specifically, buses go on duty (off duty) at regular intervals; their hydrogen demand is distributed evenly during the period between the first-on-duty (first-off-duty) bus and the last-on-duty (last-off-duty) bus. Generally, half of the buses refuel before duty and half after duty. Therefore, assuming the two periods of time are $T_1=t_{lod}-t_{fod}$ and $T_2=t_{lfd}-t_{ffd}$, the number of buses refueling in HRS can be expressed as (10):

$${\kappa }_{bus}^i(t)=\frac{1}{T_1+T_2}\hspace{1em}{t}_{lod}\le t<t_{fod}\hspace{1em}or\hspace{1em}{t}_{lfd}\le t<t_{ffd}$$

(10)

2.4. Fuel Quantity

The fuel quantity is the difference value of SOC between arriving time and leaving time. We assume that FCEVs are always refueled when their SOC reaches a minimum value and leave HRS with 100% SOC. For private FCEVs and taxis, the fuel quantity depends on the vehicle’s hydrogen capacity and minimum SOC, while for buses, it is relative to hydrogen consumption per kilometer and daily mileage. Fuel quantity for three types of FCEVs is described in Equation (11) as follows:

$$\left\{\begin{array}{l}{\chi }_{pri}^i=C_{pri}^i(1-SO{C}_{pri\min }^i)\\ {\chi }_{taxi}^i=C_{taxi}^i(1-SO{C}_{taxi\min }^i)\\ {\chi }_{bus}^i=q_{bus}^i\cdot d_{bus}^i\end{array}\right.$$

(11)

2.5. Summary

Above all, the hydrogen demand for HRSs on weekdays can be expressed as Equation (12), while the hydrogen demand for private FCEVs and taxis on weekends is assumed to be 70 and 120%, respectively.

$$\begin{array}{l}{f}_{hydrogen}(t,N_1,N_2,N_3)={\displaystyle \sum _{i=1}^{i=N_1}{\kappa }_{pri}^i(t){\lambda }_{pri}^i{\chi }_{pri}^i}\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}+{\displaystyle \sum _{i=1}^{i=N_2}{\kappa }_{taxi}^i(t){\lambda }_{taxi}^i{\chi }_{taxi}^i}\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}+{\displaystyle \sum _{i=1}^{i=N_3}{\kappa }_{bus}^i(t){\lambda }_{bus}^i{\chi }_{bus}^i}\hspace{1em}t\in (0,24]\end{array}$$

(12)

The estimation algorithm for daily hydrogen demand with hourly intervals is shown in Algorithm 1.

Algorithm 1. Hydrogen demand estimation algorithm.

$\mathrm{Step} 1: \mathrm{Input}{N}_1,N_2,N_3,D$. $\mathrm{Step} 2: \mathrm{For} \mathrm{integer}{i}_1\in (0,N_1]$$, \mathrm{generate} {\kappa }_{pri}^i(t)$$\mathrm{for} {\mathrm{i}}_{\mathrm{th}} \mathrm{private} \mathrm{FCEV} \mathrm{with} \mathrm{home} \mathrm{departure} \mathrm{time} t_{pri}^{il}$$\mathrm{and} \mathrm{home} \mathrm{arrival} \mathrm{time} t_{pri}^{ia}$$; \mathrm{generate} d_{pri}^i$$; \mathrm{calculate} SO{C}_{pri\min }^i$$, {\lambda }_{pri}^i$$\mathrm{and} {\chi }_{pri}^i$.

$\mathrm{Step} 3: \mathrm{For} \mathrm{integer} \hspace{1em}{i}_2\in (0,N_2]$$, \mathrm{generate} {\kappa }_{taxi}^i(t)$$\mathrm{for} {\mathrm{i}}_{\mathrm{th}} \mathrm{hydrogen} \mathrm{taxis} \mathrm{with} \mathrm{home} \mathrm{departure} \mathrm{time} t_{\mathrm{taxi}}^{il}$$\mathrm{and} \mathrm{home} \mathrm{arrival} \mathrm{time} t_{taxi}^{ia}$$; \mathrm{generate} d_{taxi}^i$$\mathrm{and} SO{C}_{taxi\min }^i$$; \mathrm{calculate} {\lambda }_{taxi}^i$$\mathrm{and} {\chi }_{taxi}^i$.

$\mathrm{Step} 4: \mathrm{For} \mathrm{integer} \hspace{1em}{i}_3\in (0,N_3]$$, \mathrm{calculate} {\kappa }_{bus}^i(t)$$, {\lambda }_{bus}^i$$\mathrm{and} {\chi }_{bus}^i$.

$\mathrm{Step} 5: \mathrm{Calculate} \mathrm{hydrogen} \mathrm{demand} \mathrm{with} \left(12\right). \mathrm{If} D=6|7$, multiply private FCEVs’ demand by 70%, taxis’ demand by 120%.

3. Optimization Formulation

The investment and operating costs of HRSs are decided by the size of electrolyzers and storage, as well as the hydrogen demand. The costs should be minimized, subject to constraints on the operation of the system and capacities of the electrolyzers and storage [17]. Based on the fluctuation of electricity prices, extra hydrogen can be generated and be stored during the low-price hours, while less hydrogen can be generated and the stored hydrogen makes up for inadequate supplies during the high-price hours. Meanwhile, the capacities of electrolyzers and storages should be considered, since the larger the capacities are, the higher the installation costs are. Therefore, for the optimization formulation, the variables are the capacities of electrolyzers and storage, as well as the hydrogen generation scheduling. Figure 1 shows the structure of the proposed optimization model.

3.1. Objective Function

There are no transportation costs for on-site HRSs, thus the installation costs mainly include storage devices and electrolyzers. Electricity costs account for a high proportion of operating costs [7], which is the critical factor in the objective function shown in Equation (13):

$$\underset{P_{hydrogen},P_{elec}^{\max },Q_{st}^{\max }}{\min }{C}_{elec}+C_{st}+C_{e,operation}+C_{s,operation}$$

(13)

where $P_{hydrogen}$ represents the hydrogen generation plan, $C_{elec}$ and $C_{st}$ represent the installation investment of electrolyzers and storage devices, and $C_{e,operation}$ and $C_{s,operation}$ represent the operating costs related to the electrolyzers and storage devices. $P_{hydrogen}$, $P_{elec}^{\max }$, and $Q_{st}^{\max }$ are the variables. The optimization objective is the total cost of HRS.

It is worth mentioning that we indirectly optimize the cost of the whole life cycle of the HRS by optimizing the cost for one year. The present value interest factors of annuity $(A/P,i,n)$ are applied to calculate the annual installation cost. Moreover, since the operation schedule is regulated with hourly intervals to make use of the peak-valley difference of electricity prices, daily operating costs should be calculated and summed up annually. Thus, uniform hierarchical time discretization [24] is applied to decompose the time domain into different sub-domains: hourly intervals $h\in H$, daily intervals $d\in D$, weekly intervals $w\in W$, and yearly intervals $y\in Y$.

3.1.1. Installation Costs

Usually, the installation cost for electrolyzers is approximate to the first-order function with the capacity of the electrolyzers. Similarly, for the storage devices, it is proportional to the capacity of the hydrogen storage. Therefore, the installation costs are decided by the two variables $P_{elec}^{\max }$ and $Q_{st}^{\max }$ as Equation (14).

$$\left\{\begin{array}{l}{C}_{elec}+C_{st}=n(A/P,r,n)({\beta }_{P\max }{P}_{elec}^{\max }+{\beta }_{Q\max }{Q}_{st}^{\max })\\ (A/P,r,n)=\frac{r\times {(1+r)}^n}{{(1+r)}^n-1}\end{array}\right.$$

(14)

3.1.2. Operating Costs

Electricity from the power grid is applied for both hydrogen generation and hydrogen compression with unit compression cost ${\omega }_c{c}_e$. The operating costs of hydrogen storage devices depend on the hydrogen flow and unit daily storage cost $c_s$.

$$C_{e,operation}={\displaystyle \sum _{w=1}^{\left|W\right|}{\displaystyle \sum _{d=1}^{\left|D\right|}{\displaystyle \sum _{h=1}^{\left|H\right|}(c_e{P}_{EL,w,d,h}}}}+{\omega }_c{c}_e{Q}_{in,w,d,h})$$

(15)

$$C_{s,operation}={\displaystyle \sum _{w=1}^{\left|W\right|}{\displaystyle \sum _{d=1}^{\left|D\right|}(c_s(Q_{in,w,d}+Q_{out,w,d})}})$$

(16)

Above all, the annual cost for an HRS is presented as the objective function in this chapter, concerning hydrogen generation scheduling, the capacity of the electrolyzers, and storage. $P_{EL,w,d,h}$ is related to variable $P_{hydrogen}$, while they are both influenced by the hydrogen demand of HRS.

If hydrogen storage capacity is very small, $P_{EL,w,d,h}$ should realize real-time tracking to hydrogen demand, which may lead to high electricity costs. In the opposite situation, if hydrogen storage and electrolyzer capacity are large enough, hydrogen can only be produced during the time of the lowest electricity price per day. However, this may result in an expensive installation investment.

Therefore, the objective of the proposed method is to find a balance between installation investment and operating costs, contributing to minimize total costs.

3.2. Operation Constraints

The relationship between variables should be modeled to solve the optimization formula. There are various operation constraints for HRS including electrolyzer constraints, hydrogen demand constraints, storage constraints, and grid constraints [17], of which the first three are critical constraints for the scenario in this paper. Hundreds of intermediate variables like daily hydrogen storage values are considered in the constraints and are forecast in the conclusion part.

3.2.1. Electrolyzer Constraints

The operating function of hydrogen generation is shown in Equation (17), where $LH{V}_{H_2}$ is the low heat value of hydrogen. For electrolyzers, the capacity constraint is expressed as Equation (18):

$$Q_{in,w,d,h}=\frac{{\eta }_{EL}{P}_{EL,w,d,h}}{LH{V}_{H_2}},\forall w\in W,d\in D,h\in H$$

(17)

$$0\le P_{EL,w,d,h}\le P_{elec}^{\max },\forall w\in W,d\in D,h\in H$$

(18)

3.2.2. Hydrogen Storage Constraints

For hydrogen storage devices, their quantity changes hourly related to hydrogen-in from electrolyzers and hydrogen-out to FCEVs, which are shown in Equation (19):

$$\left\{\begin{array}{l}{Q}_{st,w,d,h}=Q_{st,w,d,h-1}+[Q_{in,w,d,h}{\eta }_{H_2,in}-Q_{out,w,d,h}/{\eta }_{H_2,out}]\hspace{1em}\\ \forall w\in W,d\in D,2\le h\le \left|H\right|\\ Q_{st,w,d}=Q_{st,w,d-1}+[Q_{in,w,d}{\eta }_{H_2,in}-Q_{out,w,d}/{\eta }_{H_2,out}]\hspace{1em}\\ \forall w\in W,2\le d\le \left|D\right|\end{array}\right.$$

(19)

To balance the hydrogen storage of each year, year-end hydrogen storage will be emptied by transporting or any other approach. Therefore, at the beginning of each year, $Q_{st,1,1,1}=0$.

Besides operating function, Equation (20) presents maximum constraints for hydrogen flow and quantity in storage devices:

$$\left\{\begin{array}{l}0\le Q_{st,w,d,h}\le Q_{st}^{\max }\\ 0\le Q_{in,w,d,h}\le Q_{in}^{\max }\\ 0\le Q_{out,w,d,h}\le Q_{out}^{\max }\end{array}\right.\hspace{1em}\forall w\in W,d\in D,h\in H$$

(20)

3.2.3. Hydrogen Demand Constraints

HRS should satisfy hydrogen demand all the time, which acts as the balance between hydrogen-out from storage and hydrogen demand from FCEVs.

$$Q_{out,w,d,h}=f_{hydrogen}(w,d,h)$$

(21)

3.3. Summary

To sum up, with the objective of minimizing the total costs of HRS, the whole optimization formulation can be concluded as follows:

$$\mathrm{s}.\mathrm{t} \left(17\right)–\left(21\right) \underset{P_{hydrogen},P_{elec}^{\max },Q_{st}^{\max }}{\min }{C}_{elec}+C_{st}+C_{e,operation}+C_{s,operation},$$

According to the hydrogen demand estimation results, the optimization formulation can be solved to obtain the optimal configuration of devices and the hydrogen generation plan. The proposed optimization model takes installation investment and operating costs into consideration.

4. Case Study

The hydrogen estimation and optimization algorithms are developed in Matlab 7.12 and the computer programs are executed in a computer with the following specifications: Core i5-8265U, 3.40GHz CPU, 8GB RAM, and 64-b system. Figure 2 shows the scenario of the case study. The hydrogen energy output of HRS is determined by the driving and refueling behavior of the three types of hydrogen vehicles, from which HRS obtains income. There are three important parts in HRS: compressor, electrolyzer, and storage, which contribute to most installation costs and operating costs, except electricity costs. Electricity for HRS comes directly from the power grid.

4.1. Hydrogen Demand Estimation

The estimation model for hydrogen demand is performed in Section 2.

Table 1. Parameters for hydrogen demand estimation.

Private FCEVs	Value	Taxis	Value	Buses	Value
$N_1$	60	$N_2$	20	$N_3$	50
${\mu }_{pri}^d$	3.2	${\mu }_{\mathrm{taxi}}^d$	5	$t_{lod}$	5:00 a.m.
${\overline{v}}_{pri}^i$	55	$q_{taxi}^i$	0.01	$t_{fod}$	8:00 a.m.
${\mu }_{pri}^k$	8/17.6	${\mu }_{taxi}^k$	14/22	$t_{lfd}$	20:00 p.m.
${\alpha }^i$	0.2	$C_{\mathrm{taxi}}^i$	5.1	$t_{ffd}$	23:00 p.m.
$q_{pri}^i$	0.01	${\overline{SOC}}_{taxi\min }$	0.3	$q_{bus}^i$	0.034
$C_{pri}^i$	5.1	${\sigma }_{taxi}^k$	3.6	$C_{bus}^i$	20
${\sigma }_{pri}^k$	3.6	${\sigma }_{taxi}^d$	3	${\overline{v}}_{bus}^i$	45
${\sigma }_{pri}^d$	0.88	${\sigma }_{taxi}^{SOC}$	0.001

shows the parameters applied for the case study [21]. For hydrogen buses, their working hours are set from 6:00 to 23:00 [23], and the three-hour period is considered for duty intervals between first-leaving and last-leaving buses. It is worth mentioning that even though the buses are on duty for about 15 h, there may be extra time for them to stop or wait at the departure station. Thus, their driving times are assumed to be 10 h.

4.1.1. Single Private FCEV Refueling Behavior

The driving behavior and refueling behavior of each FCEV are simulated with the proposed estimation model. Figure 3 shows the state of hydrogen of six representative hydrogen cars among them and the map presents the serial number of each vehicle. The curve goes down when the corresponding car is running on the road, while the curve rises steeply when the corresponding car is being refueled. Different colors represent different cars. Applying the simulation model proposed in Section 2, each car has its own driving and refueling behavior. Daily consumption and refueling behavior are presented.

During the week, hydrogen consumption happens every day, while the refueling behaviors only take place when the hydrogen storage is lower than the expectant value of drivers. For example, the state of hydrogen for number 17 car is at the lowest at about 8:00 p.m. on Tuesday, and the driver refuels the car at about 9:00 a.m. on Friday, so that the state of hydrogen rises to the highest.

4.1.2. Taxis’ Hydrogen Demand Distribution

According to the simulation model proposed, the hydrogen demand distribution for taxis is shown in Figure 4. Previous studies indicate that hydrogen vehicles have a similar refueling behavior to traditional vehicles [25]. Researchers collected 1,342,957 refueling events of taxis to obtain the petrol consumption distribution during the day [26].

Compared with the actual distribution of petrol consumption in [26], they have envelopes of similar shape, which is the verification of the proposed estimation method for taxis.

4.1.3. Jensen–Shannon Divergence

The National Renewable Energy Laboratory recorded 18,568 fueling events of FCEVs [24] and provided actual refueling time during the day. The Jensen–Shannon divergence in Equation (22) can be applied to measure the degree of difference between two probability distributions. Thus, it is used to measure the error between the hydrogen demand estimation results and the actual results.

$$\left\{\begin{array}{l}KL(P||Q)={\displaystyle \sum _{i=1}^n{P}_i\log (\frac{P_i}{Q_i})}\\ JS(P||Q)=\frac{1}{2}KL(P||\frac{P+Q}{2})+\frac{1}{2}KL(Q||\frac{P+Q}{2})\end{array}\right.$$

(22)

where P and Q are two discrete distributions. The range of JS is from 0 to 1, and the closer it gets to 0, the smaller the error between P and Q is.

Considering the average value of the prediction results of taxis and cars in 1 week as P and the actual refueling distribution as Q, the result of JS(P||Q) is 0.029. Figure 5 shows the two hydrogen demand values, in which the shapes of the two distributions are similar.

The previous study also profiles the daily hydrogen demand of HRS [7]. JS divergence between its profile and the actual value is 0.0989, which is obviously larger than that of our prediction model. Therefore, the proposed hydrogen demand estimation model is proved to be accurate. It is worth mentioning that the total daily hydrogen demand is set as the same for the above distribution to make a comparison. Since the hourly hydrogen demand of buses mostly depends on the dispatching method, Figure 5 shows the demand distribution except for buses. The demand distribution of buses is shown in

Table 2. Hydrogen demand of buses with presented dispatching method.

Time of Day (Hour)	0:00–5:00	5:00–8:00	8:00–20:00	20:00–23:00	23:00–24:00
Hydrogen Demand of Buses (kg)	0	191.25	0	191.25	0

, according to the dispatching method presented in Section 2. Including the demand from buses, the assumed HRS’s daily hydrogen demand is about 8000 kg.

4.2. Capacity Configuration Optimization

According to the obtained hydrogen demand curve, the capacity of electrolyzers and storage can be optimized with the optimization formulations in Section 3. The optimization model is solved as linearly constrained optimization to determine the optimized capacity of electrolyzers and storage, and hydrogen generation scheduling to minimize the costs of both investment and operation.

4.2.1. Parameters Design

Table 3. Necessary parameters for capacity configuration.

Parameters	Value
${\beta }_{P\max }$	USD 454/kW [28]
${\eta }_{EL}$	0.6 [2]
$LH{V}_{H_2}$	39.72 kWh/kg [9]
${\eta }_{H_2,in}$${\eta }_{H_2,out}$	95%
$c_s$	USD 0.0746/kg/d
$Q_{in}^{\max }$	$0.2Q_{st}^{\max }$
$Q_{out}^{\max }$	$0.2Q_{st}^{\max }$
${\omega }_c$	1 kWh/kg
$n$	10
$r$	5% [17]
${\beta }_{Q\max }$	37.31 dollar/kg [17]

shows the necessary parameters for capacity configuration. Since the service lives of electrolyzers and hydrogen storage tanks are 10 years for each, we set |Y| = 10, |W| = 52, |D| = 7, and |H| = 24. The real-time electricity prices shown in Figure 6 are from the Illinois Power Company [27], which indicates the fluctuation of electricity prices during the day.

4.2.2. Mixed-Integer Linear Programming

It seems that the optimization formulas for HRS are linear optimization problems, which are effortless to solve. However, since this paper considers electricity prices with hourly hydrogen demand, it results in yearly optimization with hourly intervals, resulting in a large number of variables, and hourly hydrogen storage (24 × 364). Thus, some tips are employed to reduce the number of variables for optimization.

• Since daily hydrogen demand follows certain rules, the generation scheduling for electrolyzers can be 1-day scheduling, repeating every day. Regular generation scheduling can reduce labor costs and increase algorithm stability. Thus, there are just 24 variables for hydrogen generation scheduling.
• Equation (23) expresses hydrogen storage at D = d, H = h as hydrogen storage at D = d, H = 1 combined with hydrogen demand and hydrogen generation between the two times. Therefore, 364 variables are applied for daily hydrogen storage.

$$\begin{array}{l}{Q}_{st,w,d,h}=Q_{st,w,d,1}+{\displaystyle \sum _{s=1}^{s=h-1}(Q_{in,w,d,s}{\eta }_{H_2,in}-Q_{out,w,d,s}/{\eta }_{H_2,out})},\\ \forall w\in W,d\in D,2\le h\le \left|H\right|\end{array}$$

(23)

Above all, with two variables for capacities of electrolyzers and storage, there is a total of 390 variables in the input vector.

According to the optimization formulas mentioned above, the constraints and goals are described with a matrix. Real-time hydrogen storage is calculated to make sure it is less than the capacity, contributing to the [9124*390] inequality constraints matrix with [364*390] equality constraints matrix. Figure 7 shows the mixed-integer linear programming algorithm for the optimization problem.

4.3. Optimization Results and Comparison

As mentioned above, the optimization objective for the proposed configuration method is to minimize the annual cost for HRS, considering both installation cost and operation cost, while previous studies usually consider just one of them [8]. This chapter analyzes the optimization results and compares them with the configuration methods, just considering investment or operation.

OPT Method: Optimize the operation cost only, with prescribed capacity value.

According to the proposed estimation algorithm, the daily hydrogen demand is up to 11,000 kg; therefore, the configuration for electrolyzers and storage of HRS is set at 26,000 kg/d and 22,000 kg [7]. With this fixed configuration, the optimization formula is applied to minimize the operation cost.

ISL Method: Optimize the installation cost without regard to operating cost.

Configuration of HRS is optimized to minimize the installation cost, but the operation rules are irrelevant to electricity price. Thus, the operation cost is not optimized.

OPT-ISL Method: Our proposed method optimizes the installation cost and operation cost.

Since the daily operation rules depend on the fluctuation of electricity prices, the operation cost can be minimized. Meanwhile, the installation cost is taken into consideration in the optimization objective.

Table 4. Optimization results and capacity configuration.

	OPT Method	ISL Method	OPT-ISL Method
Electrolyzer capacity (kW)	71,720	26,145	41,770
Storage capacity (kg)	22,000	12,907	13,901
Electrolyzer investment (USD/year)	4,216,700	1,537,100	2,455,800
Storage investment (USD/year)	106,290	62,366	67,168
Electricity cost (USD/year)	5,851,700	8,060,100	6,318,700
Other operation cost (USD/year)	227,640	467,260	467,260
Total cost (USD/year)	10,402,000	10,127,000	9,309,100
End-year hydrogen rest (kg)	8284

shows the optimization results of the above three methods. To ensure the accuracy of comparison, the three cases end with the same amount of hydrogen remaining in storage annually and the same scenario of hydrogen demand. Even though the ISL method performs better in electricity costs, for the whole cost the proposed OPT-ISL method is the least, 8.1% less than the ISL method and 10.5% less than the OPT method. Moreover, placing a 71 MW electrolyzer in a fueling station is impossible, which means it should be located at the service station. Therefore, it seems the traditional configuration method cannot satisfy the demand for hydrogen refueling in the future.

The different hydrogen generation schedulings of the three methods are shown in Figure 8. For the OPT method, the electrolyzer works in more centralized periods to use the cheapest electricity to generate hydrogen, which may lead to a large capacity of electrolyzer and storage. For the ISL method, hydrogen is generated evenly hourly to minimize the installation costs but increase the electricity costs. Compared with previous methods, the proposed OPT-ISL method balances the installation cost and operation cost, making full use of the peak-valley difference of electricity prices and considering the appropriate capacity of hydrogen generation electrolyzers and storage.

As shown in Figure 9, the hydrogen storage of HRS can be forecast with optimal configuration for devices and a hydrogen generation plan for electrolyzers. To be precise, the hydrogen storage of HRS changes all the time and is difficult to forecast; however, applied with the proposed OPT-ISL method, hourly hydrogen generation, and demand are obtained, the daily storage change is thus acquired. It is also important to keep hydrogen storage within the permissible limits of capacity.

The percentage of HRS expenses in a year with the three configuration methods is shown in Figure 10, which illustrates that whichever method is applied to plan an HRS, the installation cost is a non-negligible part of the total cost and the electricity cost is always the most expensive part. Therefore, as shown in Table 4, the proposed OPT-ISL method takes both of them into consideration, contributing to more than a 40% decrease in electrolyzer investment compared with the OPT method and more than a 20% decrease in electricity cost with the ISL method.

4.4. Profit Analysis

It is noteworthy that hydrogen compressors play an important role in installation cost, which is not the optimized goal in this paper. According to the control variable method, the installation cost of the hydrogen compressor is assumed to be 30% of the investment for electrolyzers and storage devices with the proposed OPT-ISL method, so that the profits can be calculated. The operating costs of compression electricity consumption are included in Equation (15).

Annual hydrogen demand is 3,051,500 kg in the prescribed scenario. Hydrogen price per kilogram is adjustable from USD 2.45 to USD 6.61 [29,30], which is influenced by the electrolyzer’s efficiency, electricity price, and so on. The revenue and break-even years of the three configuration methods can be obtained with the hydrogen price decided. Present-value interest factors of annuity (PVIFA) are considered when we analyze profitability.

As shown in Figure 11, the proposed OPT-ISL method shortens the break-even time for HRS from 11.1 years to 7.8 years, decreasing 29.7%. Even though the ISL method shows much less initial investment, it results in 10.7 years for HRS to achieve break-even. Therefore, the OPT method and ISL method cannot earn their costs within 10 years with hydrogen at USD 3.5/kg.

It is important to decide on a reasonable hydrogen price for the sake of a sustainable market [31]; normally the price of hydrogen depends on investment during construction and operation. Break-even time indicates income from investment, which can probably help with hydrogen price setting.

The relationship between hydrogen price and break-even time of the three methods is shown in

Table 5. Break-even Years with Different Hydrogen Price.

Hydrogen Price (USD/kg)	Break-Even Time (Years)
	OPT Method	ISL Method	OPT-ISL Method
3	20.7	>100	15.6
3.5	11.1	10.7	7.8
4	7.9	5.8	5.4
4.5	6.0	3.9	4.1
5	4.9	3.0	3.3
5.5	4.1	2.4	2.8
6	3.6	2	2.4
6.5	3.1	1.7	2.1

.

As shown in Table 4, a higher hydrogen price contributes to a shorter break-even time. However, if the hydrogen price is less than USD 3/kg, the ISL method will be impractical for configuration. Figure 12 indicates that the proposed OPT-ISL method has the largest range of hydrogen prices when the break-even time is set less than 10 years. With hydrogen prices between USD 3.2/kg and USD 4.3/kg, the profits with the OPT-ISL method are optimal, while with hydrogen prices higher than USD 4.3/kg, the ISL method appears a little better. However, from the perspective of market economics, high hydrogen prices go against the popularity of FCEVs.

The long-term cost target for dispensed hydrogen is USD 2–3/kg [29]. Therefore, the proposed OPT-ISL method is beneficial to realize the target, helping to promote the hydrogen market. Meanwhile, it contributes to considerable profit as long as the price is set higher than USD 3.5/kg, the profit ratio can be larger than 20%.

5. Conclusions

In the aspect of the economy of HRS, capacity configuration and generation scheduling are critical research directions. This paper proposes an optimization method (OPT-ISL method) for HRSs with on-site electrolytic production, aiming to minimize the total costs of HRS.

First, a refueling behavior simulation method for three kinds of FCEVs, private cars, taxis, and buses, is creatively developed. With the simulation results of drivers’ travel habits, three critical variables are abstracted out: possibility on the road, refueling rate, and fuel quantity. The hourly hydrogen demand for an HRS with given served vehicles can be estimated with the three critical variables.

Second, an optimization model for HRS is proposed, taking advantage of the peak-valley prices of electricity from the grid and the hydrogen demand estimation results. On the one hand, the installation costs and hourly operating costs are considered; on the other hand, a variety of constraints including electrolyzers, hydrogen storage, and hydrogen demand are analyzed.

Finally, the case study validates the accuracy of the hydrogen demand estimation method and the economical efficiency of the proposed OPT-ISL method. The result of the Jensen–Shannon divergence is applied to verify the accuracy of the hydrogen demand estimation, which is much smaller than that of the estimation method in reference. In comparison with previous configuration methods that consider only installation cost or operating cost, the OPT-ISL method achieves an obvious reduction of total costs. Simultaneously, the break-even time of HRS is evidently shortened with the proposed OPT-ISL method. In other words, the OPT-ISL method achieves a lower hydrogen price with the expectant break-even time, which is beneficial to the sustainable development of the hydrogen energy market.

To sum up, this paper creatively proposes a hydrogen demand estimation method. Based on the estimation results, the proposed OPT-ISL method provides optimal configuration of devices and hydrogen generation plans for future HRS, as well as guidance on hydrogen prices.