Multi-period Planning for Hydrogen Distribution Network in the Transport Sector

(1)

University of Groningen/Newcastle University

Multi-period Planning for

Hydrogen Distribution Network in the

Transport Sector

Hang Yin

MSc Double Degree in Operations Management Master Thesis

(2)

Abstract

Insufficient refuelling infrastructure limits the adoption of hydrogen fuel cell vehicles in the transport sector. Constrained by many practical issues, a dynamic approach is required for the step-by-step development of hydrogen refuelling station networks. This research aims to develop a mathematical model which is able to maximize the staged profits of a hydrogen distribution network by optimally locating hydrogen stations in each planning period. Two customised greedy heuristic algorithms, node weight greedy algorithm and proximity greedy algorithm, which based on Forward-myopic method are developed in order to provide effective and efficient solutions for the model. Results show the proximity greedy algorithm performs well in most situations. The node weight greedy algorithm can also provide satisfying solutions within short solution time.

Keywords: Facility location-allocation model; Hydrogen refuelling stations; Hydrogen

(3)

Chapter 1. Introduction

According to the 2018 BP Energy Outlook, fossil fuels remain the dominant energy source for the transport sector and more than 70% of transport energy consumption is ascribed to fossil-fuel-based vehicles (BP p.l.c., 2018). However, the high degree of dependence on fossil fuels results in considerable environmental impacts. In 2010, 14% of the total greenhouse gas (GHG) emissions were corresponded to the transport sector, of which 72% were emitted by road transport (PICC, 2014). Taking into consideration GHG emissions, along with air pollution health hazards, fossil fuel depletion (Höök and Tang, 2013) and political issues, it is necessary to transit gradually to a more sustainable transportation system with alternative-fuels. Among many alternative options, hydrogen is increasingly gaining attention and considered as a promising solution (Demirbas, 2008), mainly due to its environment-friendliness and high energy density on a mass basis (Marbán and Valdés-Solís, 2007).

Hydrogen has been used for industrial purposes for decades. With the development of hydrogen production techniques (Mori and Saito, 2004) and the increasing demand of renewable energy, several studies have been conducted and the potential benefits of using hydrogen as a fuel have been confirmed (Demirbas, 2008; Dincer, 2007; Marbán and Valdés-Solís, 2007; Sinigaglia et al., 2017). In terms of the transport sector, the main advantage of hydrogen fuel is the absence of GHG emissions as well as other air pollutants during conversion (Cockroft and Owen, 2007). Another compelling advantage is its high gravimetric density. The amount of energy produced by the combustion of hydrogen is higher than that released by any other fuel with the same mass (Marbán and Valdés-Solís, 2007). In comparison to diesel and gasoline, hydrogen has the gravimetric density of 120 MJ/kg, which is almost three times of that of the first two energy carriers (Sinigaglia et al., 2017).

As the advantages of hydrogen fuel are identified, car manufactures’ investment in hydrogen-fuelled vehicles and hydrogen converting technology grows. There are two major approaches to convert hydrogen into energy for automotive applications. One basic approach is to convert hydrogen directly into heat energy by using internal combustion engines, and the other most commonly used method is through reaction between hydrogen and oxygen in fuel cells to obtain electric energy (Sinigaglia et al., 2017). Although it is costlier, hydrogen fuel cells gain advantages over internal combustion engines in terms of energy efficiency, near-zero tailpipe GHG emissions and low air pollutants emissions (Cockroft and Owen, 2007).

(5)

scale of HFCVs to generate demand to be economically feasible (Ogden et al., 2014; Campíñez-Romero et al., 2018). In order to tackle the ‘chicken or egg’ dilemma described above, efficient distribution networks of hydrogen fuel and more economical techniques of HFCVs production need to be developed synchronously.

The challenge of developing hydrogen fuel distribution networks lies in its high complexity and dependence on investments (Ogden et al., 2014). The design of a hydrogen fuel distribution network can be considered as a complicated systematic project, which affects the competitiveness of hydrogen as a fuel in the transport sector compared to conventional fossil fuels. The paramount consideration in the distribution network design is locating HRSs. From the perspective of an HRS investor, HRS siting should be able to address questions about where and how many HRSs needed to be built and determine the size of each HRS, in order to maximize the profits while ensure a certain level of fuel accessibility for customers.

Many researchers (Kuby et al., 2009; Kim and Kuby, 2012; Lin et al., 2008; Melaina, 2003; Campíñez-Romero et al., 2018), institutions and government departments have been dedicated to finding the best solution for locating HRSs. In comparison to using qualitative analysis to determine location selection criteria (Melaina, 2003), developing mathematical models is considered as a more effective and efficient approach for HRS siting (Kim and Kuby, 2012; Lin et al., 2008). Furthermore, hydrogen distribution infrastructure needs to be introduced gradually over time into the transport sector due to many practical constraints (Chung and Kwon, 2015). A dynamic model which is able to decide optimal locations for new HRS in each planning period could be a potential solution for step-by-step developing a hydrogen station network.

Therefore, the main objective of this research is to develop a mathematical model which is able to decide the optimal locations for new hydrogen stations in each planning period in order to maximize the staged profits of a hydrogen distribution network and investigate possible extensions of the base model by considering more details in the hydrogen refuelling network. Meanwhile, a customised solution method is expected to be provided to make sure the models can be solved effectively and efficiently. The research question and sub-questions are formulated as follows:

Research questions

How to optimally locate new hydrogen stations in each planning period in order to maximize the staged profits of a hydrogen distribution network?

RQ1. What are the costs that need to be considered in a hydrogen distribution network? RQ2. How can the revenue of a hydrogen distribution network be maximized?

(6)

Chapter 2. Literature review

The facility location problem has been a well-established research area in the discipline of Operations Research (Melo, Nickel and Saldanha-da-Gama, 2009). With the number of Supply Chain Management related research increased, Operations Research as a tool gradually entered into the Supply Chain Management domain and helped to solve many practical issues (Chopra and Meindl, 2007). Within the scope of Supply Chain Management, facility location decisions are considered as strategic level decisions for the design of supply networks and need to be addressed effectively and efficiently (Farahani and Hekmatfar, 2009, chapter 20). Therefore, facility location models developed using Operations Research strategies have been introduced in the supply chain context.

From a supply-demand perspective, each facility that provides supply needs to be located at the optimal location with assigned capacity to fulfil a certain level of customer demand, meanwhile, each customer with demand should be allocated to the available facilities (Farahani and Hekmatfar, 2009, chapter 20). These are two objectives that location-allocation models are designed to achieve. As a vital factor for formulating a facility location model, different understanding of demand results in different thinking logic that lies in the models. Classical facility location models assume customer demands appear on nodes and these nodes are the potential locations for facilities. Typical examples are the fixed charge facility location model (Balinski. 1965; Daskin, Snyder and Berger, 2005), the p-median model (Hakimi, 1964; ReVelle and Swain, 1970), the set covering location model (Toregas et al., 1971) and the maximal covering location model (MCLM) (Church and ReVelle, 1974).

Another type of models, which distinguishes themselves significantly from classical facility location models, view demand as it appears in the form of traffic flows that pass by the facilities. This strategy has been advanced by realizing many facilities, such as gas stations and ATMs, are not always the destination of people’s trips (Hodgson, 1990). One of the earliest models based on the idea of demand flow, the flow-capturing location-allocation model (FCLM), was initiated by Hodgson (1990). FCLM, as known as the flow intercepting location model (FILM), has the structure of MCLM and tries to locate a certain number of facilities to serve demand flows as many as possible. It considers a demand flow is satisfied (captured/intercepted) as long as one facility is located on the shortest path between or on the origin-destination nodes of this flow. FCLM was later adopted to locate vehicle inspection stations (Hodgson, Rosing and Zhang, 1996) and billboards (Hodgson and Berman, 1997).

(7)

station location problem, but it failed to interpret the realistic problem in vehicle refuelling, which is, sometimes a vehicle needs to be refuelled more than once to travel a longer trip due to its limited driving range. That means more than one facility is needed to be located on one path to serve a demand flow.

As the need for sufficient alternative-fuel infrastructure grows (Campíñez-Romero et al., 2018), more and more studies for optimally deploy alternative-fuel stations have been conducted, some of them are developed specifically for locating HRSs (Kim and Kuby, 2012; Kuby et al., 2009; Lin et al., 2008; Melaina, 2003; Nicholas, Handy and Sperling, 2004). Based on the practice in conventional gas station siting, Melaina (2003) analysed this problem qualitatively and suggested four criteria for selecting hydrogen station locations. Melaina (2003) also proposed three approaches to estimate the number of hydrogen stations that is sufficient for initiation. Post et al. (2017) proposed a simulation-based approach to address the uncertainty of customer demand and derive minimum infrastructure requirements in early market phases of an alternative-fuel, in terms of the number of alternative-fuel stations.

One of the flow-based models for locating alternative-fuel stations developed by Kuby and Lim (2005), which tries to address the problem in FCLM by considering the driving range of the vehicle, is known as flow refuelling location model (FRLM). FRLM has the same objective function with FCLM and inputs, but it introduces a situation where more than one refuelling station is needed to work in combination with each other to fulfil the demand for longer round trips. Kuby and Lim (2005) also provides an algorithm to determine feasible combinations for each path. This model demonstrates its fitness for practical vehicle refuelling problems and suitable for locating alternative-fuel stations as well as locating any new refuelling stations for network-based, range-limited surface transport mode.

However, the potential locations for refuelling stations in the approach proposed by Kuby and Lim (2005) are still limited to the nodes on the network, which means if the driving range of the vehicle is less than the shortest distance between a pair of contiguous nodes, the vehicle cannot reach the refuelling station even if all nodes are located with refuelling stations. Kuby and Lim (2007) then develops three methods including the ‘mid-path segments’ method, the minimax Added Node Dispersion Problem (ANDP) and the maximin ANDP, which are able to locate refuelling stations between the nodes by adding candidate sites on the arcs. The results show that the ANDP methods provide better solutions.

(8)

that the maximum allowed deviation has a measurable effect on the optimal location decision of refuelling stations. This model could be used as a decision-support tool in the initial deployment of hydrogen refuelling station network.

Another extension of the FRLM, the multi-period flow refuelling location model (M-FRLM), is proposed by Chung and Kwon (2015). They further developed the model in the context of step-by-step introducing alternative-fuel infrastructure into the transport sector and developed three methods to solve the model based on the network expansion technique of FRLM proposed by MirHassani and Ebrazi (2013). The three methods include multi-period optimization (M-opt) method, forward-myopic (F-myopic) method, and backward-myopic (B-myopic) method. The expended network with added source nodes and sink nodes addresses the main difficulty in the original FRLM of finding the valid combinations of facilities, and it is able to reduce the solution time for large-sized networks (MirHassani and Ebrazi, 2013).

The original FRLM and its extensions mentioned above provide different insights for refuelling station location problem, but there is still a nonnegligible factor, the capacity of the refuelling stations, that has not been taking into consideration in these models. Uncapacitated models imply that a single facility can refuel all flows it captured. This assumption is possible for the initial stage of develop an alternative-fuel station network since the demand is small-scaled. However, it becomes less feasible as the demand for alternative-fuel grows over time due to the gradually maturing refuelling infrastructure. Upchurch, Kuby and Lim (2009) therefore proposed the capacitated flow refuelling location model (CFRLM) to limit the number of flows that can be refuelled at each station.

From the perspective of an HRS investor, it is important to consider every activity that involves costs and find a solution to minimize the costs in order to maximize the total profits. However, a gap can be noticed that most models reviewed above consider only the demand side as they aim to locate the facility to capture demand flows as many as possible, while the potential influence of the supplier’s location on the location decisions of the HRSs is neglected. When the delivery cost is affected by the distance between the hydrogen fuel supplier and the HRS, the HRS needs to be located as close to the supplier as possible in order to minimize the delivery cost, meanwhile, its ability to capture maximum demand needs to be kept in order to maximize the revenue. The trade-off can be obvious when locating multiple HRSs with a single supplier. Furthermore, when there is a trade-off between delivery cost and storage cost in the hydrogen fuel sully chain, a classical Economic Order Quantity (EOQ) model can be applied to minimize the costs.

(9)

(10)

Chapter 3 .

_{Model formulation}

3.1 Assumptions

The model developed in this research adopts and combines the idea of multi-period planning proposed by Chung and Kwon (2015) and the classical EOQ model. As the multi-period flow refuelling location model (M-FRLM) provides the base and the classical EOQ model is incorporated to reflect the inventory strategy in each HRS located in the distribution network, assumptions of both models need to be satisfied in the model developed here. Therefore, the following assumptions are necessarily declared. Assumption (1) to (6) are derived from the M-FRLM and assumption (7) to (9) are made to satisfy the EOQ model assumptions. Assumption (10) to (12) is also made for this model.

(1) The unique shortest path for each origin-destination (O-D) pair can be determined in the network and drivers of HFCVs are always use the shortest path rationally.

(2) For each planning period, all HFCVs have the same driving range and same fuel consumption rate, which means every HFCV can run for same distance with same amount of hydrogen fuel. The driving range and fuel consumption rate for the HFCVs can be different for different planning periods.

(3) The hydrogen fuel consumption depends only on the driving distance, which means driving speed is not considered in this model.

(4) The demand of the traffic volume on the shortest path between a pair of O-D nodes is considered as fulfilled when there is an HRS or a combination of HRSs on the path that enables the HFCV to make a round trip (travel from the origin node to the destination node, then travel back to the origin node).

(5) All HFCVs are at least half fuelled at the origin and destination nodes, which ensures that the round trip will be feasible and repeatable for the situation where there is no station located at the origin or destination node. This assumption indicates that there must be at least one HRS within the distance of half the driving range from the origin and the destination nodes.

(6) All HRSs constructed in previous planning period are well and fully functioning (fulfil demands) during all later planning periods. All new stations are immediately operational, which means no construction time is considered in this model.

(11)

in the second planning period.

(8) The HRS always places the order when its inventory level reaches zero, and the lead time between the HRS places an order and receives the order from the hydrogen supplier is assumed to be zero. It is assumed that an HRS receives the order of hydrogen fuel supply immediately after ordering it. This means no stockout in HRS is allowed.

(9) The only hydrogen supplier has sufficient capacity to fulfil all orders from all HRSs in the distribution network at any time.

(10) All nodes in the initial distribution network can be the origin node and/or the destination node. This means if there are 𝑛 nodes in the initial network, 𝑛(𝑛 − 1) pairs of O-D nodes can be determined. This assumption differs from which was made in most pervious flow refuelling location models, where 𝑛(𝑛 − 1)/2 pairs of O-D nodes are defined in a 𝑛-node network. When take into consideration of the real situation of traffic volume on a path, for example the path between node A and node B, the traffic volume from node A to node B in a specific period of time can be different from which from node B to node A. Therefore, this assumption is made to differentiate paths with different directions, for example path AB is different from path BA.

(11) There is no loss of hydrogen fuel during transportation, storage and refuelling.

(12) When a combination of HRSs can fulfil the demand on a path, it is assumed that the total demand of hydrogen fuel on this path is assigned evenly to each HRS located on this path. This means, when there are 𝑛 HRSs located on a path to satisfy all traffic volume on this path, each HRS provides ₎( of the hydrogen fuel demand on this path.

3.2 Network pre-processing

(12)

Figure 1. 25-node network example

3.2.1 Determine the unique shortest path for each O-D pair

The first procedure of network pre-processing is to determine the unique shortest path for each O-D pair in the initial network 𝐺+𝒩-._{, 𝒜1}._{2. When an initial network with a certain}

number of nodes and corresponding link lengths are given, a finite set of O-D pairs can be determined. In the 25-node network example, 25(25-1) = 600 pairs of O-D nodes can be determined. For each pair of O-D nodes, a unique shortest path should be defined. In this research, Dijkstra algorithm is used to define the shortest path for every O-D pair. However, in the 25-node network, there can be multiple paths with the same shortest length for a pair of O-D nodes. For example, the O-O-D pair 1-8, there are three paths, 1-2-4-8, 1-5-4-8 and 1-5-7-8, with the shortest length of 13. In this case, the unique shortest path will be defined using the alphabetical order. Therefore, the path 1-2-4-8 is defined as the unique shortest path for the O-D pair A-H in this research. In the initial network 𝐺+𝒩-._{, 𝒜1}._{2, 𝒩}_-._{is the set of nodes on}

path 𝑞; 𝒜1._{is the set of arcs on path}_{𝑞 formed by two nodes on this path; 𝒬 is defined as}

the set of all candidate paths to be covered in the network; 𝑜𝑟𝑑.(𝑖) defines the ordering index

of the modes on path 𝑞. For example, node B on the path 𝑞: 1-2-4-8 has the ordering index 𝑜𝑟𝑑.(2) = 2. Therefore, the set 𝒜1. can be constructed as follows:

(𝑖, 𝑗) ∈ 𝒜1._if _𝑜𝑟𝑑

.(𝑖) < 𝑜𝑟𝑑.(𝑗) ∀𝑖, 𝑗 ∈ 𝒩-., 𝑞 ∈ 𝒬

3.2.2 Construct the augmented network

(13)

refuelling logic, if the distance between two HRSs is larger than the driving range of HFCVs, the vehicles can never reach the second HRS even if it is fully fulfilled at the first HRS. That means, if the driving range of the HFCV in any planning period 𝑡 is smaller than any arc in the initial network, candidate nodes must be added on the arc to make sure the HFCV can reach the destination. On the other hand, it is also practical to consider nodes on arcs to be the potential location for HRSs, since refueling is not always the purpose for vehicles travelling.

A modified minimax model is used in this research to add candidate nodes on arcs in the network. Minimax model is developed by Kuby, Lim and Upchurch (2005) to solve the added-node dispersion problem (ANDP). It is designed to locate a certain number of added-nodes on arcs to minimize the length of the longest arc in the network. However, the number of nodes to be added is unknown in this research. Therefore, a modified minimax model is developed to add the minimum number of nodes on arcs to make sure no arc in the network is longer than the driving range of HFCV. It changes the objective to objective function (1), which minimize the number of nodes added into the initial network 𝐺+𝒩-._{, 𝒜1}._{2 and adds constrain (2) to make}

sure the maximum arc length in this network will not exceed the driving range of HFCV in planning period 𝑡. The modified minimax model is constructed below:

min 𝑝 (1) subject to: max_𝑑_.(𝑖, 𝑗) ≤ 𝑅I_{∀(𝑖, 𝑗) ∈ 𝒜1}._{, 𝑞 ∈ 𝒬, 𝑡 ∈ 𝒯 (2)} max_𝑑_.(𝑖, 𝑗) ≥ 𝑆MN ∀𝑎, 𝑟 = 0, 1, 2, ⋯ , 𝑝 (3) S 𝑆_MN T NUV = 𝑑_M ∀𝑎 (4) S S 𝑋_MN T NU( M = 𝑝 (5) 𝑋MN ≥ 𝑋M(NZ() ∀𝑎, 𝑟 = 1, 2, ⋯ , 𝑝 − 1 (6) 𝑆MN ≤ 𝑑M𝑋MN ∀𝑎, 𝑟 = 1, 2, ⋯ , 𝑝 (7) 𝑋MN ∈ {0, 1} (8) max_𝑑.(𝑖, 𝑗), 𝑆MN ≥ 0 (9)

where 𝑝 is the number of nodes added into the initial network 𝐺+𝒩-._{, 𝒜1}._{2 ;}

max_𝑑.(𝑖, 𝑗) is the length of the longest arc (𝑖, 𝑗) in the initial network 𝐺+𝒩-., 𝒜1.2; 𝑅I is

the driving range of HFCVs in planning period 𝑡; 𝑎 is the index of arcs; 𝑟 is the index of sub arcs (𝑟 = 0 𝑡𝑜 𝑝) and associated added nodes (𝑟 = 1 𝑡𝑜 𝑝); 𝑆MN is the length of 𝑟th sub

arc of arc 𝑎; 𝑑M is the length of arc 𝑎. Decision variable 𝑋MN = 1 if an 𝑟th node is added

(14)

Constrain (3) ensures the longest arc is longer than or equals to any sub arc in the network. Constrain (4) ensures the total length of all sub arc on each arc 𝑎 equals to the length of arc 𝑎. Constrain (5) ensures that exactly 𝑝 nodes are added into the initial network 𝐺+𝒩-._{, 𝒜1}._2.

Constrain (6) ensures all nodes added are in order of 𝑟. Constrain (7) ensures that all sub arc on arc 𝑎 will not be longer than the length of arc 𝑎. Constrain (8) defines the binary decision variable. Constrain (9) defines the sub arc length and the maximum length of arcs as nonnegative.

The augmented network with added candidate nodes is denoted as 𝐺(𝒩._{, 𝒜}._{). The}

added nodes in the network are represented using lower letters. The total number of paths in the augmented network is the same as in the initial network, since only the nodes in the initial network can be considered as the origin or destination nodes for HFCVs’ travels. If any candidate node 𝑧 is added on a shortest path 𝑞, the sets of nodes on related paths 𝒩._needs

to be updated as follows:

𝒩. _{= 𝒩}_-._{∪ {𝑧} if 𝑑}

.(𝑖, 𝑗) > 𝑅I, 𝑧 is added on arc (𝑖, 𝑗) ∈ 𝒜1., 𝑞 ∈ 𝒬

The ordering index of nodes and the set of arcs on related paths also needs to be updated as follows:

𝒜._{= 𝒜1}._{∪ {(𝑖, 𝑧), (𝑧, 𝑗)}} _if _𝑜𝑟𝑑

.(𝑖) < 𝑜𝑟𝑑.(𝑧), 𝑜𝑟𝑑.(𝑧) < 𝑜𝑟𝑑.(𝑗), 𝑞 ∈ 𝒬

The procedure of the extended minimax method is presented in the following figure.

(15)

3.2.3 Construct the expanded network

The last procedure of network pre-processing is to construct the expanded network 𝐺+𝒩d._{, 𝒜e}._{2 to screen out all feasible combinations of nodes on each path which ensure that}

round trips on the path are possible for HFCVs. The expanded network technique is developed by MirHassani and Ebrazi (2013). This innovative approach constructs a set contains only the feasible siting solutions for each path, substantially reduces the domain for searching the optimal locating solution. The steps of this approach are illustrated in figure:

Step 1. On path 𝑞, add a source node 𝑠 before the origin node O and a sink node 𝑘 after the destination node D, which makes the set of nodes on path 𝑞 in the expanded network 𝒩d. _{= 𝒩}._{∪ {𝑠, 𝑘}. Connect the source node 𝑠 to the origin node O and the destination}

node D to the sink node 𝑘 by adding the arc (𝑠, 𝑂) and arc (𝐷, 𝑘) into the set of arcs on path 𝑞 in the expanded network, which makes 𝒜e._{= {(𝑠, 𝑂), (𝐷, 𝑘)}.}

Step 2. Connect the source node 𝑠 to any node 𝑖 on path 𝑞 that is located within the half range of the HFCV from the origin node O in planning period 𝑡. Add the arc (𝑠, 𝑖) into the set of arcs on path 𝑞 in the expanded network. That is, 𝒜e. _{= 𝒜e}._{∪ {(𝑠, 𝑖)}, if}

𝑑_.(𝑂, 𝑖) ≤j_lk, ∀𝑖 ∈ 𝒩._.

Step 3. Similar to Step 2, connect any node 𝑗 on path 𝑞 that is located within the half range of the HFCV to the destination node D in planning period 𝑡 to the sink node 𝑘. Add the arc (𝑗, 𝑘) into the set of arcs on path 𝑞 in the expanded network. That is, 𝒜e._{= 𝒜e}._∪

{(𝑗, 𝑘)}, if 𝑑_.(𝑗, 𝐷) ≤j_lk, ∀𝑗 ∈ 𝒩._.

Step 4. Connect any two nodes 𝑖 and 𝑗 on path 𝑞, if the node 𝑖 is located before node 𝑗 and the distance between these two nodes is within the range of the HFCV in planning period 𝑡 . That is, 𝒜e. _{= 𝒜e}._{∪ {(𝑖, 𝑗)}, if 𝑜𝑟𝑑}

.(𝑖) < 𝑜𝑟𝑑.(𝑗) ∩ 𝑑.(𝑖, 𝑗) ≤ 𝑅I, ∀𝑖, 𝑗 ∈

𝒩._.

Step 5. Repeat Step 1 to Step 4 for every path 𝑞, 𝑞 ∈ 𝒬. Then construct a new set of arcs 𝒜̅._{of path}_{𝑞 with the arc from the source node 𝑠 to the sink node 𝑘 in the expended}

network. That is, 𝒜̅._{= 𝒜e}._{∪ {(𝑠, 𝑘)}, 𝑞 ∈ 𝒬.}

3.3 Mathematical formulation of the main model

After the network pre-processing, the main model is presented as follows:

(16)

𝑆𝐷𝐹_tI_{= 𝑟}I_S𝑑.𝑓.I+1 − 𝑥•‚.I2 ∑ 𝑦_„I „∈…† .∈𝒬‡ ∀𝑡 ∈ 𝒯, 𝑖, 𝑗 ∈ 𝒩.₍₁₁₎ 𝐸𝑂𝑄_tI_{= ˆ}2𝐶jI𝑆𝐷𝐹tI 𝐶_|I = ˆ 2𝑟I_𝐶 jI 𝐶_|I ∙ Š S 𝑑.𝑓.I+1 − 𝑥•‚ .I₂ ∑ 𝑦_„I „∈…† .∈𝒬‡ ∀𝑡 ∈ 𝒯, 𝑖 ∈ 𝒩.₍₁₂₎

According to the EOQ model, (1) can be reduced to:

𝑚𝑎𝑥 𝑃 = S rS 𝑦_tIr+𝐶_xI_{− 𝐶} yI − 𝑑VtI 𝐶}I − 𝐶VI2𝑆𝐷𝐹tI− ‹2𝐶jI𝐶|I𝑆𝐷𝐹tI• t∈𝒩 • I∈𝒯 (13) (13) equals to: 𝑚𝑎𝑥 𝑃 = S ⎝ ⎜ ⎛ S 𝑦tI ⎝ ⎜ ⎛ (𝐶_xI_{− 𝐶} yI− 𝑑VtI 𝐶}I − 𝐶VI)𝑟I∙ S 𝑑.𝑓.I+1 − 𝑥•‚.I2 ∑ 𝑦_„I „∈…† .∈𝒬‡ − ‹2𝐶jI𝐶|I𝑟I∙ ŠS 𝑑.𝑓.I+1 − 𝑥•‚.I2 ∑ 𝑦_„I „∈…† .∈𝒬‡ ⎠ ⎟ ⎞ t∈𝒩 ⎠ ⎟ ⎞ I∈𝒯 (14) subject to: S 𝑥_t„.I {„|(t,„)∈𝒜̅†_} − S 𝑥_„t.I {„|(„,t)∈𝒜̅†_} = “ 1, 𝑖 = 𝑠.I −1, 𝑖 = 𝑘.I 0, 𝑖 ≠ 𝑠.I_{, 𝑘}.I ∀𝑡 ∈ 𝒯, 𝑞 ∈ 𝒬, 𝑖 ∈ 𝒩d.₍₁₅₎ S 𝑥_„t.I {„|(„,t)∈𝒜̅†_} ≤ 𝑦_tI ∀𝑡 ∈ 𝒯, 𝑖 ∈ 𝒩.₍₁₆₎ 𝑥_t„.I ≥ 0 ∀𝑡 ∈ 𝒯, 𝑞 ∈ 𝒬, ∀(𝑗, 𝑖) ∈ 𝒜̅.₍₁₇₎ S 𝑦_tI t∈𝒩 = 𝑛_I ∀𝑡 ∈ 𝒯 (18) 𝑦_t• _{≥ 𝑦} tI 𝜏 = 𝑡 + 1, … , 𝑇, 𝑡 ∈ 𝒯 (19) 𝑦_tI ∈ {1,0} ∀ 𝑖 ∈ 𝒩, 𝑡 ∈ 𝒯 (20) Sets

𝒬 the set of all candidate paths to be covered (the shortest paths) 𝒬t the subset of 𝒬 that contains all the paths passing node 𝑖

𝒩-. _{the set of nodes of path}_{𝑞 in the initial network 𝐺+𝒩}_-._{, 𝒜1}.₂

𝒜1. _{the set of arcs of path}_{𝑞 in the initial network 𝐺+𝒩}_-._{, 𝒜1}.₂

𝒩 the set of all path-related nodes in the augmented network 𝐺(𝒩._{, 𝒜}._{) ( 𝒩 =}

⋃ 𝒩.

.∈𝒬 )

(17)

𝒜. _{the set of arcs of path}_{𝑞 in the augmented network 𝐺(𝒩}._{, 𝒜}.₎

𝒩d. _{the set of nodes of path}_{𝑞 in the expanded network 𝐺+𝒩}_d._{, 𝒜e}.₂

𝒜e. _{the set of arcs of path}_{𝑞 in the expended network 𝐺+𝒩}_d._{, 𝒜e}.₂

𝒜̅._{the set of arcs of path} _{𝑞 in the expended network 𝐺+𝒩}_d._{, 𝒜e}._{2 with additional arcs}

including the arc from the source node 𝑠 to the sink node 𝑘 𝒯 the set of time period 𝑡

Parameters

𝐶_xI _{the selling price of per unit hydrogen fuel in planning period}_𝑡

𝐶_yI _{the buying price of per unit hydrogen fuel in planning period}_𝑡

𝐶_jI _{the reorder cost of per order of hydrogen fuel in planning period}_𝑡

𝐶_|I _{the holding cost of per unit hydrogen fuel in planning period}_𝑡

𝐶_}I_{the delivery cost of per unit distance per unit hydrogen fuel from hydrogen supplier to an}

HRS in planning period 𝑡

𝐶_VI_{the operation cost of selling per unit hydrogen fuel by an HRS in planning period}_𝑡

𝑑_VtI _{the distance between the hydrogen supplier and the HRS located on node}_𝑖

𝑑_. the length of path 𝑞 (the shortest path length between a pair of O-D nodes) 𝑓_.I _{the traffic volume on path}_{𝑞 in planning period 𝑡 (in number of HFCVs)}

𝑜𝑟𝑑_.I_{(𝑖) the ordering index of node 𝑖 of path 𝑞 in planning period 𝑡}

𝑟I _{the hydrogen fuel consumption rate of HFCV in planning period}_𝑡

𝑅I _{the driving range of HFCVs in planning period}_𝑡

𝑛_I the number of HRSs in the distribution network in planning period 𝑡

Decision variables

𝑦_tI_{= ›}1, 𝑖𝑓 𝑎𝑛 𝐻𝑅𝑆 𝑖𝑠 𝑙𝑜𝑐𝑎𝑡𝑒𝑑 𝑜𝑛 𝑛𝑜𝑑𝑒 𝑖 𝑖𝑛 𝑝𝑙𝑎𝑛𝑛𝑖𝑛𝑔 𝑝𝑒𝑟𝑖𝑜𝑑 𝑡

0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

𝑥_t„.I = the flow on an arc (𝑖, 𝑗) which belongs to 𝒜̅._{for all}_{𝑞 ∈ 𝒬 in planning period 𝑡}

If 𝑥_•‚.I, an arc flow variable on the pseudo arc (𝑠, 𝑘), is 0 on path 𝑞, then path 𝑞 is covered. If 𝑥_•‚.I = 1, then path 𝑞 is not covered.

Other notations

𝑠.I _{the source node of path}_{𝑞 in the expended network 𝐺}£_{in planning period}_𝑡

(18)

𝑆𝐷𝐹_tI_{the demand of hydrogen fuel for the HRS located on node}_{𝑖 in planning period 𝑡}

3.3.1 Objective

The main model developed in this research takes the perspective of HRS investors and aims to maximize the total profits gained over the planning horizon by optimally locating a giving number of HRSs in each planning period. In objective function (10), only one type of income and five types of costs are considered in this model. For each planning period, the income of each HRS 𝐷_tI_𝐶

xI is calculated by the selling price of per unit hydrogen fuel times

the demand of traffic volume that this HRS fulfilled. Five types of costs include the acquisition cost, the reorder cost, the holding cost, the delivery cost and the operation cost. For each HRS in each planning period, the acquisition cost 𝐷_tI_𝐶

yI is the buying price times the fulfilled

demand. The reorder cost only concerns total order time in a certain planning period and can be calculated by the ordering cost of per order 𝐶_jI_{multiplied by order times}_𝐷

tI/𝐸𝑂𝑄tI. The

holding cost relates to the inventory level and can be obtained by multiplying the cost of holding per unit hydrogen fuel in this planning period time 𝐶_|I_{with the average inventory}

level 𝐸𝑂𝑄_tI_{⁄ . The delivery cost is the cost of deliver per unit of hydrogen fuel per unit}₂

distance 𝐶_}I_{times the fulfilled demand}_𝐷

tI and the distance from hydrogen supplier to the

HRS 𝑑_VtI _{. The operation cost}_𝐷

tI𝐶VI concerns all other costs occurred during the business

operation. All prices are considered to be fixed in each planning period but may change in different planning period.

Station Demand Fulfilled (SDF) function (11) has the base of the M-FRLM and adapts the idea of Vehicle Kilometer Travelled (VKT) model which is proposed by Upchurch, Kuby and Lim (2009). Instead of maximizing the total amount of demand fulfilled in the network as VKT model does, SDF focuses on the demand of traffic volume that fulfilled by each HRS in each planning period. Comparing with VKT model, another difference is that longer paths with less HRSs will have the priority in SDF, since it is assumed in assumption that the demand on a path is assigned evenly to HRSs located on the path.

Function (12) calculates the EOQ of each HRS in each planning period. According to the classic EOQ model, the costs related to ordering batch size, that are the reorder cost and the holding cost, can be calculated in combination. That reduces objective function (10) to function (13). The complete objective function is presented in function (14).

3.3.2 Constrains

(19)

path are assigned number 0 to represent that they act as transhipment nodes. The passing flow constrain (16) ensures that a flow passes a node only when there is an HRS located on this node. Constrain (17) is the non-negative flow constrain.

(20)

Chapter 4. Model solution technique

It is very difficult to solve the main model developed in this research due to its NP-Hardness. For each planning period, if 𝑛 HRSs are to be located in a 𝑚-node network, there will be A)_¦_⁄_{𝑛! =} _{𝑚(𝑚 − 1)(𝑚 − 2) ⋯ (𝑚 − 𝑛 − 1) 𝑛(𝑛 − 1)(𝑛 − 2) ⋯ 1}_⁄ _{locating plans.}

If 𝑝 nodes are needed to be added in the augmented network, the scale of locating plans will increases to A)_¦ZT_{⁄ . For example, when 6 HRSs are to be located in the 25-node network}_𝑛!

for a certain planning period, if the driving range of HFCV is 4, 18 candidate nodes are added to the network, which results in more than 6 million locating plans for this planning period. Moreover, the number of locating plans grows almost exponentially when taking into consideration more planning periods. Even though the network expanding technique has reduced the scale of the locating plans dramatically, a heuristic algorithm is needed to solve the model more effectively and efficiently.

In this research, three methods are developed for solving the model based on forward-myopic algorithm which Chung and Kwon (2015) proposed to solve the multi-period planning problem. When using the forward-myopic method, the multi-period planning problem is divided into several single-period planning problems. The optimal locating solution is generated for each planning period separately and all planning periods are solved in sequence. This means, the location selected in previous planning periods will be fixed and become a part of the solution for later planning periods. Based on the logic of forward-myopic method, traversal algorithm, node weight greedy algorithm and proximity greedy algorithm are developed to obtain the optimal solution for a certain planning period.

4.1 Forward-myopic + Traversal

When using traversal algorithm, the candidate locating plans are generated in sequence according to the node index. If 𝑛 nodes are to be located in a planning period, the first candidate plan will determine node 1 as the first location, then the algorithm will make sure the next location selected in this plan has an index larger than the previous one. Therefore, the first set of candidate plans generated for a 𝑚 -node network will be (1, 2, … , 𝑛 − 1, 𝑛), (1, 2, … , 𝑛 − 1, 𝑛 + 1) , (1, 2, … , 𝑛 − 1, 𝑛 + 2) , … , (1, 2, … , 𝑛 − 1, 𝑚) . The next set of candidate plans will be (1, 2, … , 𝑛, 𝑛 + 1) , (1, 2, … , 𝑛, 𝑛 + 2) , (1, 2, … , 𝑛, 𝑛 + 3) , … , (1, 2, … , 𝑛, 𝑚).

(21)

temporary solution, until all plans are checked and then the temporary solution will become the optimal solution for this planning period.

This is a time-consuming method which can obtain the exact optimal solution for a certain planning period. However, it cannot guarantee to obtain the global optimal solution when taking into consideration all planning periods. Also, when there is no plan can satisfy a path in the first planning period, this algorithm always selected the nodes with the least ordering index. This situation can be expected when the driving range of HFCVs is too small and the number of HRSs to be located is too few.

4.2 Forward-myopic + Node weight greedy

A more efficient solving technique designed in this research is the node weight greedy algorithm. When using the node weight algorithm, a weight related to VKT and the number of nodes needed for a path will be assigned to every node in the augmented network. All paths will be checked and determined if they can be fulfilled with the number of HRSs to be located in this planning period. For every path, there is a minimum number of HRSs needed to fulfil it. If the number of nodes needed for a path is less than the number of HRSs to be located in this planning period, extra nodes will be added according to its weight to generate a candidate plan. The procedures for assigning weights to nodes and generate candidate plans are as follows: Step 1. Select all feasible combinations for a path with the minimum number of nodes,

allocate the total vehicle kilometre travelled (the sum of traffic volume of both direction times the path length) of this path equally to each feasible combination selected. For each feasible combination, assign its VKT evenly to nodes in this combination. For each node in the augmented network, the node weight is the sum of VKT it is assigned from feasible combinations of paths which contains this node.

Step 2. Sort all nodes in the augmented network in descending order according to the node weight.

Step 3. For each path, check if its minimum number of nodes needed, for example ℎ, is equal to or smaller than the number of HRSs to be located in this planning period 𝑛_I− 𝑛_I¨(. a) If ℎ > 𝑛_I− 𝑛_I¨(, check next path.

b) If ℎ ≤ 𝑛_I− 𝑛_I¨(, select all feasible combinations of this path which contain less than or equal to ℎ nodes, merge each feasible combination with nodes located with HRSs in last planning period.

(22)

b) If 𝑔 < 𝑛_I− 𝑛_I¨(, check if the node with the maximum weight, for example node 𝑖, is in the combination. If node 𝑖 is contained in the combination, check next node. Otherwise, add node 𝑖 to the combination. Keep adding nodes until the combination contains 𝑛_I− 𝑛_I¨( nodes and save the combination with added nodes as a candidate plan.

Step 5. Organize all candidate plans and Remove duplicates.

When a set of candidate plans are generated, calculate the profits for each plan and save the plan with greater profits as the temporary solution until all plans are checked and the optimal solution is generated. Instead of checking nodes in sequence, the node weight greedy algorithm designed in this research starts from sets of nodes on the same path, which increases the probability of fulfilling a path. This algorithm puts nodes on longer paths with greater traffic volume in priority. When adding one more node cannot fulfil more paths or catch more traffic volume, this algorithm will select the node with maximum weight among the rest nodes, expecting to fulfil path with greater traffic volume in later planning periods.

4.3 Forward-myopic + Proximity greedy

Another logic to increase the probability of fulfilling more paths is to locate HRSs on nodes link to a feasible combination. When a feasible combination of a path contains nodes less than the number of HRSs to be located in a planning period, this proximity greedy heuristic algorithm adds the closest nodes to the feasible combination. When there are no paths can be satisfied with the number of HRSs to be located in the first planning period, proximity greedy algorithm will consider the supplier location as the initial candidate node, considering lowering the delivery cost in later planning periods.

The procedures to generate candidate plans using proximity greedy algorithm are as follows:

Step 1. For each path, check if its minimum number of nodes needed, for example ℎ, is equal to or smaller than the number of HRSs to be located in this planning period 𝑛_I− 𝑛_I¨(. a) If ℎ > 𝑛_I− 𝑛_I¨(, check next path.

b) If ℎ ≤ 𝑛_I− 𝑛_I¨(, select all feasible combinations of this path which contain less than or equal to ℎ nodes, merge each feasible combination with nodes located with HRSs in last planning period.

Step 2. Check the number of nodes (𝑔) in each combination. a) If 𝑔 = 𝑛_I− 𝑛_I¨(, save the combination as a candidate plan.

(23)

add (𝑛_I− 𝑛_I¨(− 𝑔) nodes to the combination. Save the combination with added nodes as a candidate plan.

Step 3. Organize all candidate plans and Remove duplicates.

(24)

Chapter 5. Computational experiments

In this chapter, several computational experiments are conducted to evaluate the performance of two greedy algorithms ‘F-myopic + Node weight’ and ‘F-myopic + Proximity’ developed in this research by comparing with the exact optimal solution generated by ‘F-myopic + Traversal’ algorithm. Furthermore, the influence of different parameter values and different change patterns of parameters on the locating solutions are tested.

All three algorithms are coded in MATLAB R2018a and all computational experiments are conducted using computers in Groningen University. The solution time may vary due to the changing of user population on university network.

5.1 performance test

In this section, the efficiency and effectiveness of the node weight greedy algorithm and the proximity algorithm are compared with traversal algorithm in terms of different network sizes, number of HRSs located in each planning period and total number of planning periods. The efficiency of greedy algorithms is evaluated using solving time and the number of candidate plans generated. The effectiveness is analyzed in terms of profits, HRS location, traffic volume coverage rate and path coverage rate.

5.1.1 Network size

Firstly, the algorithms are tested in 9-node network, 13-node network and 21-node network. Each network tested is a part of the widely used 25-node network. For each test, there are 3 planning periods and 3 HRSs to be located in each planning period. The supplier is located on node 5 and the driving range is 4. All other parameters are the same for all tests and all planning periods. The results are shown in Table 1 and Table 2.

Table 1.

Comparison of efficiency in different network size settings. 𝑛 = (3, 6, 9)©_,_𝑅I_{= 4, supplier location: node 5}

Network size Solving time No. of candidate plans

Original Augmented No. of paths F-myopic Traversal F-myopic Node weight F-myopic Proximity F-myopic Traversal F-myopic Node weight F-myopic Proximity 9-node 16 72 0.2997 0.0487 0.1493 966 43 152 13-node 28 156 3.2770 0.1016 0.4206 7116 75 292 21-node 38 420 19.9987 0.3189 1.4476 19941 166 589

(25)

of time compared with ‘F-myopic + Traversal’. The difference of solving time between different algorithms largely due to their ways of generating candidate plans. It can be seen in the table that the number of candidate plans generated by ‘F-myopic + Traversal’ is extremely larger than which generated by other two algorithms.

Table 2.

Comparison of effectiveness in different network size settings. 𝑛 = (3, 6, 9)©_,_𝑅I_{= 4, supplier location: node 5}

Network size

F-myopic + Traversal vs. F-myopic + Node weight No. of different HRSs

Profits Traffic volume coverage Path coverage in No. Path coverage in length T1 T2 T3 9-node 0 0 1 137678.4 2.34% 5.56% 3.67% 13-node 0 0 1 187810.57 1.41% -1.28% -0.79% 21-node 0 0 3 215721.06 3.20% 1.90% 0.41%

In all three network size settings, ‘F-myopic + Proximity’ successfully obtained the optimal solutions as ‘F-myopic + Traversal’ did. However, ‘F-myopic + Node weight’ generated inferior solutions although it showed the advantage in covering more paths and longer paths in the 13-node network setting.

5.1.2 Number of HRSs located in each planning period

Another key influencer of the solving time is the number of HRSs located in each planning period. Three algorithms are tested in the 25-node network with the supplier located on node 11. All other parameters are the same for all tests and all planning periods. The results are shown in Table 3 and Table 4.

Table 3.

Comparison of efficiency when locating different number of HRSs in each planning period. Network size: 25-node network, 𝑅I_{= 4, supplier location: node 11}

No. of HRSs located in each planning period

Solving time No. of candidate plans

F-myopic _Traversal _{Node weight}F-myopic _ProximityF-myopic F-myopic _Traversal _{Node weight}F-myopic _ProximityF-myopic

𝑛 = (1, 2, 3)© _0.1515 _0.1266 _0.1397 ₁₂₆ ₁₅ ₃₄

𝑛 = (2, 4, 6)© _2.0247 _0.1900 _0.2616 ₂₄₆₄ ₇₃ ₁₁₂

𝑛 = (3, 6, 9)© _38.7603 _0.4249 _3.1135 ₂₉₉₉₁ ₁₇₂ ₈₇₇

𝑛 = (4, 8, 12)© _407.4938 _0.6264 _4.1126 ₂₅₈₀₂₁ ₂₃₃ ₁₁₂₂

𝑛 = (5, 10, 15)© _2929.0388 _1.0867 _6.0091 _1701876 ₃₂₉ ₁₂₇₃

(26)

Table 4.

Comparison of effectiveness when locating different number of HRSs in each planning period. Network size: 25-node network, 𝑅I_{= 4, supplier location: node 11}

No. of HRS in each planning period No. of different HRSs Profits Traffic volume coverage Path coverage in No. Path coverage in length T1 T2 T3 F-myopic Traversal vs. F-myopic Node weight 1 0 0 0 0.00 0.00% 0.00% 0.00% 2 0 0 0 0.00 0.00% 0.00% 0.00% 3 0 1 1 1323955.46 4.78% 2.67% 0.63% 4 0 4 1 452444.34 -0.77% -1.33% -0.40% 5 2 4 3 4111088.96 -8.08% -5.33% -2.55% F-myopic Traversal vs. F-myopic Proximity 1 0 0 0 0.00 0.00% 0.00% 0.00% 2 0 0 0 0.00 0.00% 0.00% 0.00% 3 0 0 0 0.00 0.00% 0.00% 0.00% 4 0 4 4 109812.67 6.72% 1.33% 0.34% 5 0 0 1 404697.82 -2.26% 0.67% 0.29%

In Table 4, both ‘F-myopic + Node weight’ and ‘F-myopic + proximity’ can obtain the optimal locating plans when 1 and 2 HRSs are to be located in each planning period. When 3 HRSs are to be located in each planning period, ‘F-myopic + Node weight’ cannot perform well as ‘F-myopic + proximity’ does. When locating 4 and 5 HRSs in each planning period, ‘F-myopic + Node weight’ obtains greater coverage rate in traffic volume, number of paths satisfied, and length of paths covered, while a considerable amount of profits is lost. In comparison to ‘F-myopic + Node weight’, ‘F-myopic + proximity’ performs better in settings where more HRSs are needed to be located in each planning period.

5.1.3 Number of planning periods

Several computational experiments are also conducted to test the performance of three algorithms in settings with different number of planning periods. The 25-node network is used for the experiments. The supplier is located on node 11 and all other parameters are kept he same for all planning periods. The results are shown in Table 5.

Table 5.

Comparison of efficiency in different number of planning periods setting.

Network size: 25-node network, 𝑅I_{= 4, supplier location: node 11, number of HRSs located in each planning period: 3}

Total number of planning periods

Solving time No. of candidate plans

(27)

It can be noticed in Table 5 that, for all three algorithms, adding planning periods in a certain network setting will increase the solving time. However, the influence of number of planning periods on solving time is not as significant as that in pervious experiments. In all network settings with different number of planning periods, ‘F-myopic + Proximity’ can obtain the optimal solution found by ‘F-myopic + Traversal’. As can be seen in Table 6, ‘F-myopic + Node weight’ does not perform perfectly for all planning periods, but for some planning periods, for example planning periods 4 and 5 in Table 6, it can obtain a solution which provides more profits and greater coverage rate than ‘F-myopic + Traversal’ does.

Table 6.

Comparison of effectiveness in 𝑛 = (3, 6, 9, 12,15,18)©

Network size: 25-node network, 𝑅I_{= 4, supplier location: node 11}

Planning period

F-myopic + Traversal vs. F-myopic + Node weight No. of different HRS Profits Traffic volume coverage Path coverage in No. Path coverage in length T1 0 0.00 0.00% 0.00% 0.00% T2 3 505943.04 -0.56% -0.67% -0.49% T3 1 818012.42 4.78% 2.67% 0.63% T4 2 -156677.63 -5.03% 1.00% 0.46% T5 2 -161730.47 16.26% 7.33% 2.35% T6 1 313015.05 -0.66% -2.00% -0.88% Sum 1318562.42

From the experiments results presented above, it can be concluded that ‘F-myopic + Node weight’ is the most efficient algorithm in terms of solving time, however, it cannot always obtain the optimal solution. There also a possibility that ‘F-myopic + Node weight’ obtains better solution with more profits and greater coverage rate than ‘F-myopic + Traversal’ does, but this is a rather random situation. Overall, ‘F-myopic + Proximity’ can perform better than ‘F-myopic + Node weight’ in most network settings while not losing the efficiency. There is a high probability for ‘F-myopic + Proximity’ to obtain the optimal solution. ‘F-myopic + Traversal’ can always find the exact optimal locating plan for current planning period, however the solving time can be extremely long when the problem scale is big.

5.2 Scenario test

(28)

5.2.1 Supplier location

Supplier location is a new parameter that is introduced in this research. To investigate the influence of supplier location on the locating decisions, 9 nodes are selected as locations of supplier in different scenarios. In the 25-node network setting with driving range equals 4 and 𝑛 = (3, 6, 9)©_{, the optimal solutions obtained by ‘F-myopic + Traversal’ for scenarios where}

supplier is located on node 6 and node 11 are presented in Figure 4. It can be observed that node close to the supplier location tend to be chosen in early planning period. This situation can be explained by the existing of delivery cost. The shorter the distance between the supplier and HRSs, the less needs to be spent on delivery.

Figure 4. Optimal solution when supplier is located on node 6 and node 11 Table 7.

Scenario test: supplier location. Supplier

F-myopic + Traversal vs. F-myopic + Node weight

Location Weight rank No. of links Time in seconds No. of different HRS Profit Traffic volume coverage Path coverage in No. Path coverage in length No. of candidate plans T1 T2 T3 14 1 5 40.9538 0 0 0 0.00 0.00% 0.00% 0.00% 29806 20 6 3 40.0125 0 0 0 0.00 0.00% 0.00% 0.00% 29839 17 11 3 41.2786 0 1 3 1123822.85 4.50% -1.33% -0.09% 29838 11 16 4 39.5460 0 3 1 1323955.46 4.78% 2.67% 0.63% 29819 5 21 4 37.3068 0 0 2 215721.06 2.24% 1.33% 0.25% 29814 22 2 2 40.5572 0 0 0 0.00 0.00% 0.00% 0.00% 29806 13 13 5 38.6508 0 0 0 0.00 0.00% 0.00% 0.00% 29830 6 24 2 37.1374 0 0 0 0.00 0.00% 0.00% 0.00% 29819 25 25 1 35.8924 2 3 4 351764.87 0.93% 0.00% 0.04% 29821

(29)

weight. Table 7 shows the comparison between solutions obtained by ‘F-myopic + Traversal’ and ‘F-myopic + Node weight’. ‘F-myopic + Node weight’ performs undesirably when the supplier is located on node 25 which has the least node weight in the 25-node network. The number of links, that is the number of nodes directly connected to the supplier node, does not have a major influence on the effectiveness of algorithms tested.

In the main model developed in this research, there is only one supplier with fixed location in all planning period. This is not always the case in real situation. If the hydrogen supplier and hydrogen refuelling infrastructure are owned and operated by different organization, for the purpose of cost saving, HRSs investor will select the most economical strategy to purchase hydrogen. Multiple hydrogen supplier could be chosen, and the demand occurred in HRSs could be allocated to different supplier. In the case where the supplier and HRSs achieve cooperation, other ordering strategy could be made by HRS investor in order to maximize the supply chain profits. On the other hand, if the hydrogen supplier and the refuelling infrastructure are vertically integrated, both the supplier and HRSs could be located in the network in order to achieve long-term profits for the organisation.

5.2.2 Change of driving range

Another set of scenarios are tested with different driving range values and different change patterns of driving range. With the network setting where supplier located on node 11, 𝑛 = (2, 4, 6, 8,10)©_{, all three algorithms are able to generate the optimal locating plans in all}

scenarios. Figure 5 shows the profits obtained by the optimal locating plan in each planning period. It is obvious that more profits can be created when the driving range of HFCVs is larger in a certain period.

Figure 6 shows the profits made the optimal solution when the driving range of HFCVs increases over time. In this set of tests, the driving range starts with the same value in planning period 1 and increases to the same value in planning period 5. In the situation that driving range grows with a decreasing rate, that is, the increment of driving range is larger in earlier planning period, more profits can be made in this network setting. That is expected since the sooner a larger driving range can be obtained, the more paths can be satisfied with existing HRSs. In another situation where driving range is decreasing over time, more profits are generated with an increasing change rate. In all three scenarios, the driving range decreases from a same value and ends in the last planning period with a same small value. The results shown in Figure 7 is expected since the decrement of driving range is smaller in the setting with increasing change rate.

(30)

Figure 5. Scenario test: different driving ranges

Figure 6. Scenario test: driving range increases over time

Figure 7. Scenario test: driving range decreases over time

5.2.3 Change of traffic volume and fuel consumption rate

Next, scenarios with different traffic volume and different change patterns of traffic volume are tested. All three algorithms are able to generate the optimal locating plans in all scenarios. Results are shown in Figure 8, 9 and 10. It can be observed form Figure 8 that no matter of large or small traffic volume, profits will grow in later planning period. The difference in the first planning period is not significant, but the profits grow quicker in the scenario with larger traffic volume. Conclusions similar with scenarios of driving range can be drawn according to results shown in Figure 9 and Figure 10. In scenarios where traffic volume increases with a decreasing rate and decrease with an increasing rate, greater profits can be expected.

Similar results can be obtained in scenarios testing for fuel consumption rate. Results are shown in Figure 11, Figure 12 and Figure 13.

0 10000000 20000000 30000000 40000000 Change rate increasing

Change rate constant Change rate decreasing

Driving range increases over time

T1 T2 T3 T4 T5

Driving range decreases over time

T1 T2 T3 T4 T5 0 10000000 20000000 30000000 T1 T2 T3 T4 T5

Profits in settings of different driving range

(31)

Figure 8. Scenario test: different traffic volumes

Figure 9. Scenario test: traffic volume increases over time

Figure 10. Scenario test: traffic volume decreases over time

Figure 11. Scenario test: different fuel consumption rates

Figure 12. Scenario test: fuel consumption rate increases over time 0

10000000 20000000

T1 T2 T3 T4 T5

Profits in settings of different traffic volume

50% 100% 150%

Traffic volume increases over time

T1 T2 T3 T4 T5

0 2000000 4000000 6000000 8000000 10000000 Change rate increasing

Traffic volume decreases over time

T1 T2 T3 T4 T5

Fuel consumption rate increases over time

T1 T2 T3 T4 T5 0

10000000 20000000

T1 T2 T3 T4 T5

profits in settings of different fuel consumption rate

(32)

Figure 13. Scenario test: fuel consumption rate decreases over time

5.2.4 Change of delivery price

In the scenario tests of delivery price, interesting results can be obtained. In Figure 15, when delivery cost increases over time, the profits will decrease dramatically in later planning periods. It is expected since the initial delivery cost is low in earlier planning periods, nodes far from supplier location might be selected which could result in huge delivery costs in later planning period. However, Figure 14 shows that if the delivery price is already high in earlier planning periods and keeps stable over time, greater profits can be achieved.

Figure 14. Scenario test: different delivery price

Figure 15. Scenario test: delivery price increases over time Figure 16. Scenario test: delivery price decreases over time 0 20000000 40000000 60000000 80000000 100000000 T1 T2 T3 T4 T5

Delivery price decreases over time

Change rate increasing Change rate constant Change rate decreasing -60000000 -40000000 -20000000 0 20000000 T1 T2 T3 T4 T5

Delivery price increases over time

Change rate increasing Change rate constant Change rate decreasing

0 10000000 20000000 30000000 40000000 50000000 60000000 70000000 T1 T2 T3 T4 T5

profits in settings of different delivery price

0.01 1 2

0 4000000 8000000 12000000 Change rate increasing

Fuel consumption rate decreases over time

(33)

Chapter 6. Conclusion and future research

This research investigated the multi-period locating problem for hydrogen refuelling stations in the transport sector. A mathematical model is developed to maximize the total profits gained over the planning horizon by optimally locating a giving number of HRSs in each planning period. The main model introduces the hydrogen fuel supplier into the network and the classical EOQ model is incorporated to generate the optimal ordering strategy for each HRS.

Two customised greedy heuristics algorithms, node weight greedy algorithm and proximity greedy algorithm, which based on Forward-myopic method are developed in order to provide effective and efficient solution for the model. Result shows that both the proximity greedy algorithm performs well in most situations. The node weight greedy algorithm can also provide satisfying solutions within short time.

(34)

Reference

Balinski. M. (1965) ‘Integer programming, methods uses and computation’, Management Science, 12, pp.253-265.

BP p.l.c. (2018) BP Energy Outlook, 2018 Edition. London, United Kingdom. Available at:

https://www.bp.com/content/dam/bp/en/corporate/pdf/energy-economics/energy-outlook/bp-energy-outlook-2018.pdf [Accessed: May 2nd_{, 2018].}

Campíñez-Romero, S., Colmenar-Santos, A., Pérez-Molina, C. and Mur-Pérez, F. (2018) ‘A hydrogen refuelling stations infrastructure deployment for cities supported on fuel cell taxi roll-out’, Energy, 148, pp.1018-1031.

Chopra, S. and Meindl, P. (2007) Supply Chain Management: Strategy, Planning and

Operations. Prentice Hall, New Jersey.

Chung, S.H. and Kwon, C. (2015) ‘Multi-period planning for electric car charging station locations: A case of Korean Expressways’, European Journal of Operational Research, 242(2), pp.677-687.

Church, R.L. and ReVelle, C. (1974) ‘The maximal covering location problem’, Papers in

Regional Science, 32, pp.101-118.

Cockroft, C.J. and Owen, A.D. (2007) ‘The economics of hydrogen fuel cell buses’, Economic

Record, 83, pp.359-370.

Daskin, M.S., Snyder, L.V. and Berger, R.T. (2005) ‘Facility location in supply chain design’. In: Langevin, A. & Riopel, D. (eds) Logistics systems: Design and operation. Springer, New York, pp.39–66.

Demirbas, A. (2008) ‘Present and future transportation fuels’, Energy Sources, Part A:

Recovery, Utilization, and Environmental Effects, 30(16), pp.1473-1483.

Dincer, I. (2007) ‘Environmental and sustainability aspects of hydrogen and fuel cell systems’,

International Journal of Energy Research, 31, pp.29-55.

Farahani, R.Z. and Hekmatfar, M (Eds.) (2009) Facility location: Concepts, models,

algorithms and case studies, Physica Verlag, Heidelberg, Germany.

Hakimi, S.L. (1964) ‘Optimum locations of switching centres and the absolute centres and medians of a graph’, Operations Research, 12, pp.450-459.

Hodgson, M.J. (1990) ‘A flow capturing location-allocation model’, Geographical Analysis, 22, pp.270-279.

Hodgson, M.J. and Berman, O. (1997) ‘A billboard location model’, Geographical and

(35)

Hodgson, M.J., Rosing, K.E. and Zhang, J. (1996). ‘Locating vehicle inspection stations to protect a transportation network’, Geographical Analysis, 28, pp.299–314.

Höök, M. and Tang, X. (2013) ‘Depletion of fossil fuels and anthropogenic climate change — A review’, Energy Policy, 52, pp.797-809.

Karlsson, C. (2016) Research Methods for Operations Management, 2nd_{edition, Routledge:}

New York.

Kim, J.G. and Kuby, M.J. (2012) ‘The deviation-flow refueling location model for optimizing a network of refueling stations’, International Journal of Hydrogen Energy, 37(6), pp.5406-5420.

Kuby, M.J. and Lim, S. (2005) ‘The flow-refuelling location problem for alternative-fuel vehicles’, Socio-Economic Planning Sciences, 39, pp.125–45.

Kuby, M.J. and Lim, S. (2007) ‘Location of alternative-fuel stations using the flow-refueling location model and dispersion of candidate sites on arcs’, Networks and Spatial

Economics, 7(2), pp.129-152.

Kuby, M.J., Lines, L., Schultz, R., Xie, Z., Kim, J.G. and Lim, S. (2009) ‘Optimization of hydrogen stations in Florida using the Flow-Refueling Location Model’. International

journal of Hydrogen Energy, 34(15), pp.6045–6064.

Kuby, M., Lim, S., & Upchurch, C. (2005) ‘Dispersion of nodes added to a network’,

Geographical Anal, 37(4), pp.383-409.

Lim, S. and Kuby, M.J. (2010) ‘Heuristic algorithms for siting alternative-fuel stations using the flow-refueling location model’, European Journal of Operational Research, 204(1), pp.51–61.

Lin, Z., Ogden, J., Fan, Y. and Chen, C.W. (2008) ‘The fuel-travel-back approach to hydrogen station siting’, International Journal of Hydrogen Energy, 33(12), pp.3096–3101.

Marbán, G. and Valdés-Solís, T. (2007) ‘Towards the hydrogen economy?’, International

Journal of Hydrogen Energy, 32(12), pp.1625-1637.

Melaina, M.W. (2003) ‘Initiating hydrogen infrastructures: Preliminary analysis of a sufficient number of initial hydrogen stations in the U.S.’, International Journal of Hydrogen

Energy, 28(7), pp.743–755.

Melo, M.T., Nickel, S. and Saldanha-da-Gama, F. (2009) ‘Facility location and supply chain management – A review’, European Journal of Operational Research, 196(2), pp.401-412.

Multi-period Planning for Hydrogen Distribution Network in the Transport Sector