On the medication distribution system for home health care through convenience stores,lockers,and home delivery

Abstract

Medication distribution service can be delivered based on a combination of home delivery and customer pickup. That is, medications are delivered either to customers’ homes directly or to the pickup facilities (e.g. lockers) close to customers’ homes. In Taiwan, there are more than 11,000 convenience stores that provide a 24-h service for customers to pick up the ordered items from e-commerce, which is unique to the world. In the medication distribution system, convenience stores can provide a unique opportunity for customers to more conveniently collect medications at stores, and also can reduce the operating cost for a logistics company providing the medication delivery service. Therefore, this work proposes a medication distribution system through convenience stores, lockers, and home delivery. Under this system, this work investigates how to simultaneously determine employment of convenience store chains, the convenience store locations to be visited, locations of lockers, vehicle routes for convenience stores and lockers, and vehicle routes for customers’ homes, so that the total operating cost is minimized. This work further proposes a genetic algorithm to solve the medication distribution problem. Through simulation, the experimental results show that the proposed algorithm is able to solve the problem efficiently.

Keywords

genetic algorithm home delivery home health care logistics location routing problem locker medication distribution

Introduction

Home health care has been implemented in various medical care services,^1,2 and has been widely supported by the medication distribution based on home delivery services, which delivers medications from a pharmacy directly to customers’ homes.^3,4 Generally, this convenient medication distribution service is much suitable for patients with chronic diseases taking maintenance medications⁵ and patients with other regular-use prescriptions. However, such a medication distribution service based on home delivery requires visiting all patients’ homes by vehicles, which cause enormous transportation cost. In business, logistics companies have been striving to reduce the transportation cost to increase their competitiveness in a competitive market. To reduce the transportation cost, some logistics companies have developed efficient distribution systems that deliver items through pickup facilities (e.g. lockers) instead of distributing to customers’ homes directly,^6–8 and thus vehicles do not need to visit all patients.

A recent work by Veenstra et al.⁶ has investigated the health care logistics of medication distribution in which lockers (which are automatic self-service cabinets and are available for 24 h) enable customers to receive medications if customers’ homes are within the coverage distance of those lockers; otherwise, medications are delivered to customers’ homes directly. Lockers are visited by vehicles to replenish medications with an ideal situation that it must be replenished before a given time period of the day and before customers return home. Thus, the vehicle route that visits lockers is separated from the vehicle route that needs to visit customers’ homes when they are at homes. The work by Veenstra et al.⁶ considered a number of potential locations where lockers can be installed, and then investigated a joint facility location and vehicle routing problem which simultaneously determines the locations of installing lockers, the vehicle routes visiting all installing lockers, and the vehicle routes visiting all the customers that are not covered by any lockers.

To the best of our understanding, the 24-h convenience store (CS) pickup service in Taiwan is unique to the world. With this service, customers can pick up the ordered items at CSs whenever they are available after they receive pickup messages. In addition to providing flexible time for picking up items, this service is advantageous because an enormous number of CSs can support the arising number of customers. In 2019, the major CS chains in Taiwan (with an area of 35,808 km²) include 7-Eleven™ (having 5579 branches), FamilyMart™ (having 3406 branches), Hi-Life™ (having 1380 branches), and OK Mart™ (having 902 branches). On average, there is a CS every 0.31 km² in Taiwan. The CS density in Taiwan is only lower than that in South Korea, and is ranked second in the world. Therefore, it has been popular and extensively used for customers in Taiwan to pick up items through CSs, especially for the items ordered from e-commerce. A unique opportunity for using CSs in the medication distribution system is that customers can more conveniently collect their medications at 24-h CSs instead of going to hospitals or pharmacies in person. Multiple medication packages for different customers with close home locations can be jointly distributed to a CS just once, so that the transportation cost can be reduced. Furthermore, the total cost can be decreased when existing CSs are employed instead of constructing new pickup facilities.

In light of the above, this work proposes a novel medication distribution system through CSs, lockers, and home delivery (Figure 1). The proposed system includes a depot (i.e. a logistics center of a pharmacy), customers’ homes, CSs, and lockers. For convenience of explanation, both CSs and lockers are called pickup facilities through the rest of this article. Although a lot of medication packages can be collected together within a pickup facility, this medication distribution system is secure. For using the locker pickup service, the customer requires providing a pickup code to open the locker. Similarly, for using the CS pickup service, the customer requires providing an identification document to be verified by the CS clerk. The pickup facilities are associated with a coverage distance. If a customer is within the coverage of a certain pickup facility, then this system delivers the customer’s medication to this facility; otherwise, medication is delivered to the customer’s home directly. Note that, although both lockers and CSs provide 24-h pickup services, they differ as follows: (1) CSs have existed, but lockers are considered to be installed; and (2) CSs are divided into multiple chain brands, and can join this service only when their chain has a contact with the pharmacy. Hence, employment of each CS chain (with a set of CSs) requires payment of a contract fee; whereas employment of each locker requires an installation cost, which mainly depends on the installation site.

Figure 1.

Illustration of the proposed medication delivery system through CSs, lockers, and home delivery.

Under the above assumptions, this work investigates the medication distribution problem which integrates the location routing problem (LRP), consideration of employment of CS chain brands, and the covering concept for serving customers through pickup facilities. For convenience of notation, this problem is called the medication distribution problem with multiple delivery methods (MD² for short). The MD² problem simultaneously determines how to employ CS chains with contracts, the CS locations to be visited, the locations of installing lockers, the vehicle routing for pickup facilities, and the vehicle routing for uncovered customers, so that the total operating cost (including the total contract fee, cost of installing lockers, and the total routing cost) is minimized. This work further proposes a genetic algorithm (GA) to solve the MD² problem, and then implements and evaluates the GA on various problem instances.

The main contributions of this work are as follows:

Different from the previous work by Veenstra et al.,⁶ this work proposes additionally considers CSs, which have existed and are available for being employed without further investment in construction, in the medication distribution system based on lockers and home delivery.

This work provides a GA solution to solve the MD² problem. In addition, the performance of the proposed GA is evaluated on different-scale problem instances.

Literature review

The LRP generally joins the facility location problem (FLP) and the vehicle routing problem (VRP). Because both FLPs and VRPs are NP-hard, the LRP is also NP-hard.^9–12 The essential characteristics of a standard LRP consist of the deterministic data, one planning period, a finite set of potential locations for deploying facilities, a single-objective problem, satisfaction of customers’ demand, no load transfer, visiting only once for each customer by only one vehicle, and no inventory.⁹ Furthermore, the standard LRP can be extended to solve various problems that consider additional characteristics,¹³ for example, inclusion of uncertain information,¹⁴ subcontracting options,^10,15 capacitated constraints,^16,17 pickup and delivery,¹⁸ inventory decision,¹⁹ fuzzy demand,²⁰ price-sensitive demand,²¹ and restricted time window constraints.^18,20,22

A variety of algorithms and heuristic approaches to address various LRPs have been developed. The work by Zhang et al.¹⁴ proposed a hybrid intelligent algorithm integrating uncertain simulation and GA to solve the sustainable multi-depot emergency facilities LRP with uncertain multi-objective. A hybrid heuristic method integrating simulated annealing (SA) and variable neighborhood search (VNS) was represented by Stenger et al.¹⁰ to solve the LRP with capacitated vehicles of small package shippers, in which the proposed idea of subcontract depot operations is similar to the work by Stenger et al.,¹⁵ which employed an adaptive VNS algorithm that combines an approach for routing and customer selection in the shaking step. The adaptive VNS algorithm was also developed to determine the swap station locations and the vehicle routes with capacitated electric vehicles in Hof et al.¹⁶

To solve the capacitated LRP with a tight capacity constraint on both depots and vehicles was proposed by Yu et al.¹⁷ using a hybrid GA. They developed the algorithm that allows both feasible solutions and infeasible solutions kept in two subpopulations, and then combines them into the population subsequently. The work by Capelle et al.¹⁸ presented a column generation method to solve the LRP involved the vehicle route for both picking up the items and delivering items to the destination. The decision of inventory management at the facilities was integrated with the LRP in the work of Hiassat et al.,¹⁹ which adopted the GA to determine the number and the locations of required warehouses, the inventory level at each retailer, and the routes. To solve the LRP with fuzzy demands and time windows, the work by Fazayeli et al.²⁰ adopted a two-part GA for solving a mixed-integer mathematical fuzzy model. The work by Ahmadi-Javid et al.²¹ studied the profit-maximization LRP with price-sensitive demands, and further tackled the problem by a branch-and-price algorithm. The work by Schiffer and Walther²² solved an LRP within time window by an exact approach that simultaneously determined the routing of electric vehicles and the charging infrastructure.

In general, the LRP is considered with only one level of routing within the network from depots to customers. Furthermore, the LRP can be extended to the so-called two-echelon location routing problem (2E-LRP), which has two levels of routing: the routes from depots to satellites and the routes from satellites to customers.²³ The work by Breunig et al.²⁴ tackled the 2E-LRP arising in the city logistics by a hybrid metaheuristic approach combining local and large neighborhood search. Moreover, the work by Pichka et al.²⁵ addressed a special case of the 2E-LRP, that is, the two-echelon open location routing problem (2E-OLRP) operated by the third-party logistics providers, and hence the vehicle routes did not need to return to the depot and to satellites after finishing the service. They proposed a mixed-integer linear programming method and a hybrid SA heuristic to solve this problem. The work by Zhou et al.⁷ solved the multi-depot two-echelon VRP which involves two levels of VRPs including the delivery option for the last mile distribution. They proposed a hybrid multi-population GA to improve the search efficiency and speed up the evolution process obtained by the multi-population strategy. They included two subpopulations derived from feasible solutions and infeasible solutions similar to the work by Yu et al.¹⁷

LRPs have been particularly applied in the real-world case with practical instances to increase efficient distribution networks and to increase competitiveness in the recent market. Some real-world case studies were presented by Veenstra et al.⁶ and Zhu and Ursavas,²⁶ which solved the LRP for the medication distribution in the Netherlands. The proposed model by Zhu and Ursavas²⁶ determined the facility and routing decisions simultaneously by considering three node types: depots, satellites, and customers. They formulated and solved the problem by Lagrangian relaxation, branch and cut algorithm, and upper bound heuristic approach. In their work, each customer had to be visited by vehicles, whereas the work by Veenstra et al.⁶ did not need. That is, some customers were assigned to receive their medicines at lockers instead of receiving at home directly. The function of lockers by Veenstra et al.⁶ was a temporary medication storage that is then distributed to customers. The work by Veenstra et al.⁶ focused on the medication distribution from the depot (i.e. pharmacy) to patients in order to minimize the routing cost. They considered two types of routes: locker routes (from the depot to lockers) and patient routes (from the depot to patients’ homes). Their work further proposed a branch-and-bound algorithm and a hybrid heuristic algorithm for solving this LRP.

Our problem is also related to the covering problem. In most of the covering problems, customers receive services from facilities depending on the distance between customers and facilities. Basically, the decisions on the location of a set of facilities and a set of demand points are considered.^27,28 The solutions to the covering location problem were proposed by Marín et al.²⁷ and Zarandi et al.²⁹ with exact methods for small instances and heuristics for large instances. The Lagrangian heuristic was employed to tackle the covering location problem under multi-period stochastic by Marín et al.²⁷ while the SA algorithm was proposed by Zarandi et al.²⁹ to solve the maximal covering location problem with numerous demand nodes and potential facilities. In addition, the covering tour problem, in which each facility might be visited more than once by a vehicle, was solved in the studies.^8,30–33

The work by Karaoğlan et al.³⁰ formulated an integer non-linear programming problem under probabilistic coverage, and solved it by a branch-and-cut algorithm and a local search heuristic based on VNS. The work by Pham et al.⁸ proposed a branch-and-cut algorithm comparing to GA. The formulation of the covering tour problem arising in humanitarian logistics was solved by a heuristic approach in the works of Nedjati et al.³¹ and Naji-Azimi et al.³² and a greedy randomized adaptive search procedure.³³ In Naji-Azimi et al.,³² the covering tour method was employed to locate the satellite distribution centers for the humanitarian aid distribution in a disaster area. The decision on facilities was determined to cover all customers. The covering concept by Veenstra et al.⁶ was also similar to Naji-Azimi et al.³² (i.e. the items were distributed to customers through facilities in the case that customers were covered by those facilities). However, the work by Veenstra et al.⁶ proposed that the opened lockers do not need to cover all customers because the remaining customers can receive their parcels at their homes.

Methodology

Problem description

This work proposes the MD² problem, in which medications are distributed through multiple delivery methods: CSs, lockers, and home delivery. The MD² problem aims to determine the employed CS chains with contracts, the CS locations to be visited, the locations of installing lockers, the vehicle routing for pickup facilities, and the vehicle routing for uncovered customers, so that the total operating cost (including the total contract fee, cost of installing lockers, and the total routing cost) is minimized.

Consider a problem instance with one depot, n_B CS chains/brands, n_CS CSs, n_L potential locations for installing lockers, and n_C customers. The instance can be represented as a complete graph G = (V, A), in which the node set V is partitioned as V = {0} ∪ N_CS ∪ N_L ∪ N_C, where {0} represents the depot; set N_CS includes n_CS CSs denoted by S₁, S₂,. . ., $S_{n_{C S}}$ ; set N_L includes n_L potential locker locations denoted by L₁, L₂,. . ., $L_{n_{L}}$ ; set N_C includes n_C customers denoted by C₁, C₂,. . ., $C_{n_{C}}$ ; each arc (i, j) ∈ A is associated with a distance d_ij, travel cost c_ij, and travel time t_ij. Each facility j ∈ N_CS ∪ N_L has a coverage distance r_j for serving customers. Each node j ∈ V has a service time g_j required to operate the service.

For example, in Figure 2, there are three CS chains (i.e. n_B = 3) and eight CSs (i.e. n_CS = 8) denoted by S₁–S₈, in which the first CS chain includes CSs S₁–S₃; the second CS chain includes CSs S₄–S₆; and the third CS chain includes CSs S₇ and S₈. There are eight locker locations (i.e. n_L = 8) denoted by L₁–L₈, and 23 customers (i.e. n_C = 23) denoted by C₁–C₂₃.

Figure 2.

Illustration of an example of the medication distribution problem.

The decision on the MD² problem involves with employing multiple CS chains, and hence we consider that set N_B represent n_B CS chains. Note that only the CSs from the employed CS chains can be considered to visit or not. This problem considers two types of vehicle routes: routes consisting of pickup facilities (i.e. CSs and lockers) and routes for customers’ homes. Hence, two corresponding vehicle sets are denoted by M and K, in which each vehicle m ∈ M delivers items from the depot to pickup facility locations; whereas each vehicle k ∈ K delivers items from a depot to customers’ homes.

Each customer i ∈ N_C is assigned once to pick the medication up at either a CS, a locker, or the customer’s home. That is, if customer i is within the coverage distance r_j of a pickup facility (i.e. a CS or a locker), then the customer will be assigned to receive the medication item from that facility; otherwise, the customer’s home will be identified as a node to be visited by some vehicle k ∈ K. However, the travel cost for vehicle k is penalized by a factor π that is greater than one, because the delivery through pickup facilities is preferred to the direct delivery to customers’ homes. In addition, the routing distance of each vehicle m ∈ M (respectively k ∈ K) must not exceed the maximum routing distance $r_{\max}^{1}$ (respectively $r_{\max}^{2}$ ). Similarly, the travel duration of each vehicle m ∈ M (respectively k ∈ K) (composed of the total travel time and the total service time occurred at each visited node) must not exceed the maximum routing duration T₁ (respectively T₂).

The MD² problem aims to determine both the employed facilities and the vehicle routes of both types to serve all customers. Five types of decision variables in this problem include as follows: (1) which CS chains will be employed, (2) which CSs (each of which belongs to some CS chain) will be utilized, (3) which locker locations will be selected to install lockers, (4) the vehicle routes visiting the employed facility locations, and (5) the vehicle routes visiting the uncovered customers’ homes.

With the above notations, the objective of the MD² problem is to minimize the total cost formulated as follows

\begin{array}{l} Minimize π \cdot \sum_{i, j \in N_{C} \cup {0}; i \neq j} \sum_{k \in K} c_{i j}^{c} \cdot x_{i j}^{k} + \sum_{i, j \in N_{C S} \cup N_{L} \cup {0}; i \neq j} \sum_{m \in M} c_{i j}^{f} \cdot y_{i j}^{m} \\ + \sum_{j \in N_{L}} F_{j}^{open} \cdot v_{j} + \sum_{j \in N_{B}} F_{j}^{fee} \cdot u_{j} \end{array}

(1)

The above objective function is to minimize four types of costs. The first term is the total transportation cost of the routes visiting the uncovered customers, where π is the penalty factor that is greater than one; $c_{i j}^{c}$ is the travel cost between customers i and j; $x_{i j}^{k}$ is the binary variable deciding whether vehicle k ∈ K travels from customer i to customer j. The second term represents the transportation cost of the routes visiting the employed facilities, where $c_{i j}^{f}$ is the travel cost between facilities i and j; $y_{i j}^{m}$ is the binary variable deciding whether vehicle m ∈ M travels from customer i to customer j. The third term represents the total cost of installing lockers, where $F_{j}^{open}$ is the fixed cost of installing a locker at locker location j; and v_j is the binary variable deciding whether a locker is installed at location j. The fourth term represents the total contract fee for employing CS chains, where $F_{j}^{fee}$ is the contract fee for employing CS chain j; and u_j is the binary variable deciding whether CS chain j is employed.

Overview of the proposed GA

This work proposes a GA approach to solve the MD² problem. The GA is a metaheuristic approach inspired by the process of natural selection. The GA is a powerful search technique that performs the global search characteristic, and hence is suitable to solve large-scale optimization problems with the capability of reaching near-optimum values. The proposed GA is given in Algorithm 1, which is explained as follows. Consider a population of multiple chromosomes, each of which represents a candidate solution for the problem. The GA starts from randomly creating the initial population (Line 1), and its cost is evaluated (Line 2). Then, until the maximal number of iterations is achieved, a main loop in Lines 3–16 repeatedly selects the parent chromosomes from the current population (Line 4), executes GA operators on these parent chromosomes, that is, crossover (Line 5) and mutation (Lines 6–13), and remains the better chromosomes in P (Line 15). Finally, the chromosome with the best cost in P is outputted as the solution of the GA (Line 17).

The main components are detailed in the following sections.

Algorithm 1. The proposed GA
1:	Randomly generate the initial population P
2:	Evaluate the cost of each chromosome in P
3:	while the maximal number of iterations is achieved do
4:	Apply the binary tournament selection mechanism to select a set of parent chromosomes P_parent from P
5:	Randomly generate a number from {1, 2,. . ., 5}. If the number is no greater than 3 (respectively greater than 3), the one-point crossover operator (respectively uniform crossover operator) is conducted on P_parent to generate a set of offspring chromosomes P_offspring under the crossover rate P_c
6:	Randomly select some chromosomes in P_offspring under the mutation rate P_m
7:	for each selected chromosome s do
8:	Randomly generate a number γ from {1, 2,. . ., 5}
9:	if γ is 1, 2, or 3 then
	Conduct the type-I mutation operation on the γ th part of chromosome s
10:	else (i.e. γ is 4 or 5)
11:	Conduct the type-II mutation operation on the γth part of chromosome s
12:	end if
13:	next for
14:	Evaluate the cost of each chromosome in P
15:	Replace the chromosomes with better cost in P by the chromosomes with better cost in P_offspring
16:	end while
17:	Output the best chromosome in P

Solution representation

A chromosome (i.e. a candidate solution for the MD² problem) in the proposed GA is represented as (b₁, b₂,. . ., $b_{n_{B}}$ | h₁, h₂,. . ., $h_{n_{C S}}$ | l₁, l₂,. . ., $l_{n_{L}}$ | f₁, f₂,. . ., $f_{n_{C S} + n_{L}}$ | s₁, s₂,. . ., $s_{n_{C}}$ ), consisting of the following five types of decision variables:

For i ∈ {1, 2,. . ., n_B}, binary variable b_i is one if CS chain i is contracted.

For j ∈ {1, 2,. . ., n_CS}, binary variable h_j is one if CS j is employed, but it depends on the results derived from the results of b_i (i.e. only the CSs of contracted CS chains can be employed and visited by vehicles).

For k ∈ {1, 2,. . ., n_L}, binary variable l_k is one if a locker is installed at location k.

〈f₁, f₂,. . ., $f_{n_{C S} + n_{L}}$ 〉 is a permutation of {S₁, S₂,. . ., $S_{n_{C S}}$ , L₁, L₂,. . ., $L_{n_{L}}$ } used for deciding the vehicle routes visiting facilities.

〈s₁, s₂,. . ., $s_{n_{C}}$ 〉 is a permutation of {1, 2,. . ., n_C} used for deciding the vehicle routes visiting the uncovered customers.

For example, a solution for this instance in Figure 2 shows a problem instance with n_B = 3, n_CS = 8, n_L = 8, and n_C = 23 is encoded in Figure 3. From Part 1, the first CS chain (corresponding to three CSs S₁–S₃ in Figure 2) and the second CS chain (corresponding to three CSs S₄–S₆) are contracted, but the third CS chain (corresponding to two CSs S₇ and S₈) is not. From Part 2, CSs S₂, S₃, S₄, S₅, and S₆ are employed, but S₈ is not employed and visited, because its chain is not contracted (from b₃ = 0 in Part 1). From Part 3, lockers are installed at locations L₁, L₂, L₄, L₆, and L₈. Part 4 is a permutation of {S₁, S₂,. . ., S₈, L₁, L₂,. . ., L₈}. Based on Parts 1–3, the closed facilities are marked by “×” below the permutation of Part 4 in Figure 3. Hence, the permutation of the remaining facilities is L₁–S₃–L₂–S₅–L₄–S₂–L₆–S₄–L₈–S₆. Considering some constraints for vehicles, Figure 2 shows three facility routes L₁–S₃–L₂–S₅, L₄–S₂–L₆, and S₄–L₈–S₆. Part 5 is a permutation of {1, 2,. . ., 23}. Based on the three facility routes, the customers covered by the visited facilities are marked by “v” below the permutation of Part 5 in Figure 3. Hence, the permutation of the unvisited customers is C₁–C₂–C₃–C₁₁–C₁₂–C₁₃–C₁₄–C₁₅. Considering some constraints for vehicles, Figure 2 shows two customer routes C₁–C₂–C₃ and C₁₁–C₁₂–C₁₃–C₁₄–C₁₅.

Figure 3.

Solution representation for the problem instance in Figure 2.

Cost evaluation

This work sets the cost function in GA as the objective function (1) to minimize the total cost. The fittest chromosome in GA is corresponding to a solution with the minimal cost.

Given a chromosome (b₁, b₂,. . ., $b_{n_{B}}$ | h₁, h₂,. . ., $h_{n_{C S}}$ | l₁, l₂,. . ., $l_{n_{L}}$ | f₁, f₂,. . ., $f_{n_{C S} + n_{L}}$ | s₁, s₂,. . ., $s_{n_{C}}$ ), the cost of the chromosome is evaluated in Algorithm 2, which is explained as follows. Line 1 initializes all the n_C customers as uncovered customers. Then, Lines 2–6 refer to the chromosome encoding to determine which customers are covered by some CSs or lockers. Line 7 initializes the total cost ϕ to be zero, and also initializes two variables that will be used later. Then, Lines 8–18 calculate the cost for multiple routes for uncovered patients. Similarly, Lines 19–30 calculate the cost for multiple routes for facilities. Then, Line 31 calculates the cost of opening lockers, and Line 32 calculates the cost of employing CS chains. Finally, Line 33 outputs the total cost.

Algorithm 2. Cost evaluation (b₁,. . ., $b_{n_{B}}$ \| h₁,. . ., $h_{n_{C S}}$ \| l₁,. . ., $l_{n_{L}}$ \| f₁,. . ., $f_{n_{C S} + n_{L}}$ \| s₁,. . ., $s_{n_{C}}$ )
1:	All the n_C customers are initialized as uncovered customers
2:	for each t = 1 to n_C do
3:	if customer C_t is within the coverage r_j of an opened CS S_j of a contracted CS chain i (i.e. b_i = h_j = 1 and d(C_t, S_j) < r_j) for all i ∈ {1, 2,. . ., n_B} and j ∈ {1, 2,. . ., n_CS}, or within the coverage r_k of an opened locker L_k (i.e. l_k = 1 and d(C_i, L_k) < r_k) for all k ∈ {1, 2,. . ., n_L} do
4:	Customer C_t is marked as a covered customer
5:	end if
6:	next for
7:	Set the total cost ϕ = 0, the distance so far d_s = 0, and the previous index p = 0
8:	for each t = 1 to n_C do
9:	if customer $C_{s_{t}}$ is uncovered then
10:	if d_s + d( $C_{s_{p}}$ , $C_{s_{t}}$ ) < $r_{\max}^{2}$ then
11:	Let ϕ = ϕ + π · $c_{p t}^{c}$ and d_s = d_s + $d_{p t}^{c}$
12:	else
13:	Let ϕ = ϕ + π · ( $c_{p 0}^{c}$ + $c_{0 t}^{c}$ ) and d_s = $d_{0 t}^{c}$
14:	end if
15:	p = t
16:	next for
17:	next for
18:	If at least one route for uncovered patients exists, then ϕ = ϕ + π· $c_{p 0}^{c}$
19:	Set the distance so far d_s = 0 and the previous index p = 0
20:	for each t = 1 to n_CS + n_L do
21:	if the facility corresponding to f_t is an opened locker or a CS of some employed chain then
22:	if d_s + d(f_p, f_t) < $r_{\max}^{1}$ then
23:	Let ϕ = ϕ + $c_{p t}^{f}$ and d_s = d_s + $d_{p t}^{f}$
24:	else
25:	Let ϕ = ϕ + $c_{p 0}^{f}$ + $c_{0 t}^{f}$ and d_s = $d_{0 t}^{f}$
26:	end if
27:	p = t
28:	end if
29:	next for
30:	If at least one route for facilities exists, then ϕ = ϕ + $c_{p 0}^{f}$
31:	For each k = 1 to n_L, if l_k = 1, then ϕ = ϕ + $F_{k}^{open}$
32:	For each i = 1 to n_B, if b_i = 1, then ϕ = ϕ + $F_{i}^{fee}$
33:	Output the total cost ϕ

Initializing the population

This section presents initializing the population as follows. Each gene (i.e. binary variable) in Parts 1, 2, and 3 of each chromosome in the population is randomly assigned 0 and 1. Part 4 in the chromosome is assigned a random permutation of {S₁, S₂,. . ., $S_{n_{C S}}$ , L₁, L₂,. . ., $L_{n_{L}}$ }, and Part 5 in the chromosome is assigned a random permutation of {1, 2,. . ., n_C}.

Parent selection

The chromosomes in the population are chosen for GA operations based on the binary tournament selection mechanism. The fitter chromosomes providing the lower cost are favored. The concept of binary tournament selection is that, two chromosomes are randomly selected from the population, and then the winner being a chromosome with the better cost value is chosen to be inserted into the mating pool (i.e. parent chromosomes). This mechanism can retain some good chromosomes while giving the chance for other weaker individuals to take part in mating.

Crossover operator

Crossover is an operator in the GA that recombines the gene codes of two parents, which partially contribute characteristics to new chromosomes (or called offspring). In this process, pairs of parent chromosomes in the mating pool are mated randomly with a crossover rate P_c to create new offspring chromosomes. The chromosome representation in this work includes two encoding types (i.e. the binary encoding in Parts 1, 2, and 3 for FLPs, and the permutation encoding in Parts 4 and 5 for VRPs), and hence the proposed GA conducts two types of crossover operators for the two chromosome encoding types, respectively.

For the binary encoding in Parts 1, 2, and 3, we apply the one-point crossover operator. With loss of generality, consider the one-point crossover operator at Part 1 on two parent chromosomes x₁ and x₂. First, randomly cuts each of the two concerned parent chromosomes x₁ and x₂ at Part 1 into two binary substrings, in which the two substrings of x₁ (respectively x₂) is x₁₁ and x₁₂ (respectively x₂₁ and x₂₂), and then interchanges the two substrings of the two parent chromosomes to generate two offspring chromosomes: x₁₁ x₂₂ and x₂₁ x₁₂.

For the permutation encoding in Parts 4 and 5, we apply the uniform crossover operator. First, a random binary string called mask is generalized. In the mask, if the ith bit value is 1, then the corresponding ith genes in the two parents are exchanged. The other genes masked by 0 in either parent are copied by the remaining genes in the order of the other parent to maintain the chromosome feasibility. That is, each offspring chromosome inherits some genes from one parent (when mask bits are 1) and also inherits remaining genes from the other parent (when mask bits are 0).

Mutation operator

Mutation is executed by changing a gene(s) in a chromosome to prevent the algorithm from becoming trapped in local optima and to explore new solutions in the solution space so as to increase diversity of solutions. In this mutation operator, a number of chromosomes in the population are randomly selected to be mutated with a mutation rate P_m. In each selected chromosome, we randomly determine a part to be mutated. If the selected part is one of Parts 1, 2, or 3, then the Type-I mutation operator is conducted; otherwise (i.e. Parts 4 or 5), the Type-II mutation operator is conducted. The Type-I mutation operator is to randomly choose a gene from the concerned chromosome, and then to modify this gene to 0 (respectively 1) if the original gene value is 1 (respectively 0). The Type-II mutation operator performs a swap operator, that is, two genes in the concerned chromosome are randomly chosen to be swapped with each other, while other genes keep to remain their positions.

Experimental results

This section evaluates the performance of the GA in an experimental environment. The performance on small-scale and large-scale instances is investigated, and the experimental analyses including a sensitivity analysis are conducted.

Experimental environment and parameter setting

The proposed GA is implemented in the C++ programming language, and runs on a PC with an Intel^® Core™ i7-7700, 3.60 GHz (four core) CPU, and 8 GB RAM. The parameters used in the proposed GA are set as follows: the population size is set as 40; the crossover rate P_c is set as 0.5; the mutation rate P_m is set as 0.01; and the maximal number of iterations is set as 2000.

Analysis on small-scale instances

This section describes the performance of the proposed GA on small-scale problem instance. The locations of CSs, potential lockers, and customers are randomly generated within a two-dimensional (2D) geographical region [0, 100] × [0, 100]. The depot is randomly located in the range [25, 75] × [25, 75]. The distance between all nodes and covering radius of facilities are computed as the Euclidean measure. We assume that the travel time and the travel cost are set equal to the travel distance. The other parameter setting of the problem instances is given in Table 1, in which U is the uniform distribution function.

Table 1.

Parameter setting for small-scale instances.

Parameter	Value
Number of CS chains (n_B)	4
Number of CSs (n_CS)	70
Number of CSs of each CS chain	{25, 20, 15, 10}
Contract fee of employing each CS chain $(F_{j}^{fee})$	{40, 35, 30, 25}
Cost of installing a locker $(F_{j}^{open})$	U(50, 100)
Coverage distance of a CS (r_j)	20
Coverage distance of a locker (r_k)	U(10, 20)
Maximum routing distance for a facility route $(r_{\max}^{1})$	70
Maximum routing distance for a customer route $(r_{\max}^{2})$	70
Service time at each node (g_j)	10
Penalty factor (π)	10
Maximum routing duration for a facility route (T₁)	100
Maximum routing duration for a customer route (T₂)	100

CS: convenience store.

The convergence analysis under four various parameter settings of P_c and P_m (in which n_L = 10 and n_C = 50) is shown in Figure 4. From Figure 4, the case with P_c = 0.5 and P_m = 0.01 performs best. That is, in this case, the improvement during 2000 iterations occurs up to the point where no further improvement is possible, and up to the point where the solution is in steady state, and hence this setting is applied in the rest of this work.

Figure 4.

Convergence analysis of the proposed GA.

Subsequently, we execute 100 runs of the four cases, and the results are shown in Figure 5. Although the 100 runs of each case obtain different cost values, these cost values show stability within a reasonable range.

Figure 5.

The cost values of 100 runs of the proposed GA under four parameter settings of P_c and P_m.

The experimental results of the proposed GA on 54 instances under various settings of n_C and n_L (in which n_CS is fixed) are compared in Table 2, in which the “best cost” is the best value among 10 runs of the proposed GA under the same setting, and the “time” is the CPU time taken to find the best solution. From Table 2, the proposed GA can search for the optimal solution within short CPU time. Overall, the CPU time increases slightly when the problem size grows up.

Table 2.

Comparison of the experimental results of the proposed GA on 54 small-scale instances under various settings of n_C and n_L.

Instance	n_C	n_CS	n_L	Best cost	Gap (%)	Time (s)
R1	20	70	5	260.61	11.35	4.41
R2	20	70	10	234.05	0.00	4.81
R3	20	70	15	305.73	30.62	5.25
R4	20	70	20	292.99	25.18	6.02
R5	20	70	25	298.86	27.69	6.59
R6	20	70	30	300.32	28.31	6.60
R7	30	70	5	456.84	3.27	5.59
R8	30	70	10	456.09	3.10	6.04
R9	30	70	15	519.81	17.51	6.80
R10	30	70	20	553.70	25.17	7.23
R11	30	70	25	648.07	46.50	7.80
R12	30	70	30	442.36	0.00	8.55
R13	40	70	5	490.31	13.26	6.69
R14	40	70	10	449.39	3.81	7.16
R15	40	70	15	656.47	51.65	8.25
R16	40	70	20	432.89	0.00	8.61
R17	40	70	25	448.90	3.70	9.80
R18	40	70	30	537.31	24.12	10.38
R19	50	70	5	692.29	44.02	7.72
R20	50	70	10	494.36	2.84	8.43
R21	50	70	15	533.14	10.91	9.53
R22	50	70	20	480.71	0.00	10.31
R23	50	70	25	539.35	12.20	11.18
R24	50	70	30	582.13	21.10	12.50
R25	60	70	5	457.81	0.00	8.89
R26	60	70	10	517.64	13.07	10.01
R27	60	70	15	500.24	9.27	10.85
R28	60	70	20	584.45	27.66	11.82
R29	60	70	25	485.95	6.15	11.96
R30	60	70	30	562.03	22.76	13.37
R31	70	70	5	505.54	11.90	10.20
R32	70	70	10	451.77	0.00	11.62
R33	70	70	15	523.30	15.83	11.45
R34	70	70	20	488.49	8.13	13.43
R35	70	70	25	468.14	3.62	14.49
R36	70	70	30	530.98	17.54	16.14
R37	80	70	5	528.85	4.58	11.52
R38	80	70	10	505.68	0.00	12.72
R39	80	70	15	531.16	5.04	13.12
R40	80	70	20	588.15	16.31	15.36
R41	80	70	25	543.15	7.41	16.21
R42	80	70	30	545.30	7.83	17.41
R43	90	70	5	575.89	6.76	12.64
R44	90	70	10	545.09	1.05	14.23
R45	90	70	15	652.38	20.94	15.14
R46	90	70	20	639.77	18.60	16.15
R47	90	70	25	560.15	3.84	18.46
R48	90	70	30	539.45	0.00	20.43
R49	100	70	5	555.03	7.93	14.37
R50	100	70	10	554.41	7.81	15.16
R51	100	70	15	514.25	0.00	16.74
R52	100	70	20	638.57	24.18	19.03
R53	100	70	25	590.23	14.78	19.94
R54	100	70	30	555.40	8.00	22.40

In Table 2, every six instances form a group with the same n_C value but with different n_L values. For each group of the experimental results, a gap of performance for an instance i is calculated as

Gap = \frac{(Cost of instance i_the best cost)}{(The best cost)} \times 100 %

where “the best cost” is the best cost among the six instances in the same group. From Table 2, the best cost result in each group may not occur in the case with the same n_L value. However, the best cost values derived from using n_L = 10 are often the best in their respective group as three times (i.e. the results of instances R2, R32, and R38). According to the experiment results, there is no conclusion for the effect of changing numbers of lockers because differences of random instances may lead to different related costs such as uncertain locations of customers and pickup facilities.

Analysis on large-scale instances

This section analyzes the experimental results on large-scale instances. The locations of CSs, potential lockers, and customers are randomly generated within a 2D geographical region [0, 10,200] × [0, 10,200]. The depot is randomly located at the center of this region. The other parameter setting of the large-scale problem instances is given in Table 3. The other settings are similar to those in the last section.

Table 3.

Parameter setting for large-scale instances.

Parameter	Value
Number of CS chains (n_B)	4
Number of CSs (n_CS)	308
Number of CSs of each CS chain	{138, 88, 64, 18}
Contract fee of employing each CS chain $(F_{j}^{fee})$	{40,000, 35,000, 30,000, 25,000}
Cost of installing a locker $(F_{j}^{open})$	U(10,000, 20,000)
Coverage distance of a CS (r_j)	500
Coverage distance of a locker (r_k)	U(300, 500)
Maximum routing distance for a facility route $(r_{\max}^{1})$	14,400
Maximum routing distance for a customer route $(r_{\max}^{2})$	14,400
Service time at each node (g_j)	180
Penalty factor (π)	10
Maximum routing duration for a facility route (T₁)	18,000
Maximum routing duration for a customer route (T₂)	18,000

CS: convenience store.

The results of the 54 instances under various settings of n_C and n_L (in which n_CS is fixed) solved by the proposed algorithm are compared in Table 4. This table shows that the proposed GA is able to search for the optimal solution during 10 runs within reasonable CPU time. Moreover, the CPU time increases slightly when the size of problem grows up.

Table 4.

Comparison of the experimental results of the proposed GA on 54 large-scale instances under various settings of n_C and n_L.

Instance	n_C	n_CS	n_L	Best cost	Gap (%)	Time (s)
L1	20	308	5	277,998	0.51	55.01
L2	20	308	10	305,260	10.37	57.37
L3	20	308	15	302,204	9.26	65.08
L4	20	308	20	366,154	32.39	62.65
L5	20	308	25	358,796	29.72	63.20
L6	20	308	30	276,582	0.00	71.88
L7	30	308	5	725,038	3.77	76.65
L8	30	308	10	706,998	1.19	77.95
L9	30	308	15	698,667	0.00	77.77
L10	30	308	20	726,931	4.05	84.75
L11	30	308	25	739,573	5.85	88.88
L12	30	308	30	785,172	12.38	99.20
L13	40	308	5	695,864	0.00	91.83
L14	40	308	10	742,809	6.75	89.15
L15	40	308	15	801,550	15.19	92.96
L16	40	308	20	807,035	15.98	94.65
L17	40	308	25	776,359	11.57	92.81
L18	40	308	30	893,978	28.47	97.94
L19	50	308	5	819,366	0.00	110.71
L20	50	308	10	880,949	7.52	112.61
L21	50	308	15	887,628	8.33	121.52
L22	50	308	20	929,580	13.45	127.97
L23	50	308	25	867,590	5.89	123.14
L24	50	308	30	997,599	21.75	136.83
L25	60	308	5	1,016,659	0.00	125.43
L26	60	308	10	1,061,941	4.45	130.85
L27	60	308	15	1,100,571	8.25	136.93
L28	60	308	20	1,101,812	8.38	142.44
L29	60	308	25	1,109,730	9.15	149.50
L30	60	308	30	1,253,534	23.30	158.26
L31	70	308	5	970,020	0.00	133.79
L32	70	308	10	1,154,079	18.97	143.01
L33	70	308	15	1,174,357	21.07	149.04
L34	70	308	20	1,294,203	33.42	155.33
L35	70	308	25	1,380,436	42.31	161.88
L36	70	308	30	1,342,194	38.37	166.91
L37	80	308	5	1,259,803	0.00	153.37
L38	80	308	10	1,369,819	8.73	160.97
L39	80	308	15	1,371,764	8.89	165.45
L40	80	308	20	1,382,363	9.73	183.98
L41	80	308	25	1,410,407	11.95	182.48
L42	80	308	30	1,442,512	14.50	191.75
L43	90	308	5	1,397,062	2.32	175.08
L44	90	308	10	1,377,890	0.92	178.37
L45	90	308	15	1,403,991	2.83	185.16
L46	90	308	20	1,365,327	0.00	186.27
L47	90	308	25	1,425,397	4.40	189.51
L48	90	308	30	1,481,756	8.53	204.71
L49	100	308	5	1,466,167	0.21	191.86
L50	100	308	10	1,509,689	3.18	191.33
L51	100	308	15	1,463,091	0.00	205.69
L52	100	308	20	1,484,684	1.48	206.02
L53	100	308	25	1,497,318	2.34	211.97
L54	100	308	30	1,488,898	1.76	235.96

In Table 4, every six instances form a group with the same n_C value but with different n_L values. For each group, a gap of performance for an instance is calculated as mentioned in the explanation of Table 2. The best cost result in each group may not occur in the case with the same n_L value. However, the best cost values derived from using n_L = 5 are often the best in their respective group as five times (i.e. the results of instances L13, L19, L25, L31, and L37).

According to the experiment results, it seems that there is no conclusion for the effect of changing numbers of lockers because differences of random instances may lead to different related costs such as uncertain locations of customers and pickup facilities. Furthermore, the number of existing CSs comparing to the number of defined lockers is enormous to serve customers, thus changing the number of lockers may not significantly affect cost. In practical, lockers should be installed at the places where have customers but no existing CS such as in rural area in order to support the medication delivery efficiently.

Sensitivity analysis

This section presents the impacts of changing the values of four parameters (i.e. the coverage distance of a pickup facility, the number of CSs, the contract fee of a CS chain, and the cost of installing a locker) on the cost. Each parameter is modified on four instances (i.e. 50 customers-10 lockers, 50 customers-20 lockers, 70 customers-10 lockers, and 70 customers-20 lockers). The other parameter setting of these instances follows the parameter setting of large-scale instances in the last section, except the number of CSs follows the parameter setting of small-scale instances to reduce the computational time. The parameter values are varied by multiplying the base case value (mentioned in previous sections) by 0.8 up to 2 with step of 0.2 for the coverage distance and the number of CSs, and by 0.5 up to 3.5 with step of 0.5 for the contract fee and the cost of installing a locker. The results of the sensitivity analysis are presented in Figure 6.

Figure 6.

Sensitivity analysis on the cost of four input parameters: (a) the coverage distance of a pickup facility, (b) the number of CSs, (c) the contract fee of a CS chain, and (d) the cost of installing a locker.

Figure 6(a) shows that when the coverage distance of a pickup facility increases, the cost value decreases consistently for all scenarios. It is clear that the coverage distance considerably affects the cost. Each facility can serve more customers when its coverage distance increases, so that both the number of uncovered customers and the number of employed facilities decrease. However, adjustment of this factor should be considered with customer satisfaction. Figure 6(b) shows the effect of changing the number of CSs. Even if the plot in Figure 6(b) seems to be fluctuate, the cost in all scenarios tends to drop slightly when the number of CSs increases. This indicates that the number of CSs slightly affects change of the cost.

Figure 6(c) and (d) shows the impact of changing the contract fee of CS chains and the cost of installing a locker, respectively. In this work, the contract fee is defined as a flat-rate pricing (i.e. whether the number of CSs changes or not, the price does not change); whereas the cost of installing a locker is calculated depending on the number of lockers. The cost slightly increases if the coefficient on these two factors increases for instances C3, C4, D2, D3, and D4; whereas the cost is slightly stable for instances C1, C2, and D1. The reason why the cost has a little change is that these two types of cost are small as compared with other components of the cost function.

Discussion

Although the delivery for medications has been increasingly implemented, based on the recent literature reviews, we found that only a few studies—for example, Fathollahi Fard et al.³, Veenstra et al.⁶, and Zhu and Ursavas²⁶—have solved the logistics problem of medication delivery to determine the locations of potential facilities and the vehicle routing simultaneously. In this work, we extend the framework derived from Veenstra et al.,⁶ which investigates the joint facility location and VRP for medication delivery through lockers and customers’ homes. The main difference between the previous works by Veenstra et al.⁶ and our work is that the CS pickup service is added to the medication distribution system, that is, some customers can be served to collect their medications at CSs. CSs have existed, and have a lot of branches available to serve customers. Even though the main function of CSs is similar to lockers in the work of Veenstra et al.,⁶ a medication distribution service provider does not need to invest in constructing new pickup facilities. The covering concept is also considered in this work (i.e. this work focuses on both coverage distances of CSs and lockers to serve customers). Each customer is assigned to receive medications at either a CS or a locker if the customer is located within the service coverage of the facility; otherwise, the customer is assigned to be visited directly. To the best of our knowledge, this is a first work to propose the medication distribution system through CSs, lockers, and home delivery.

The experimental results show that the proposed GA has been designed to be suited to solve our MD² problem, and it can solve this problem within reasonable computational time. The proposed GA can solve the smallest-scale instance R1 (in which n_C = 20, n_CS = 70, and n_L = 5) within 4.41 s; while it can also solve the largest-scale instance L54 (in which n_C = 100, n_CS = 308, and n_L = 30) within 235.96 s or less than 4 min. However, the experimental results are not compared with other methods. Hence, a line of our further work is to explore an approach that can provide better performance to solve this problem and achieve better solution quality, and conduct a comprehensive experimental comparison with other methods.

The feasibility to serve the medication distribution through the three pickup approaches in real-world implementation can be supported with the following reasons: (1) A lot of CS branches have existed in Taiwan, and they enable customers to access the service easily. It leads to a huge number of customers who have used the CS pickup service, especially for e-commerce. (2) In Taiwan, as compared with CSs, lockers are not the main stream to receive parcels. However, lockers enable customers who live far away from CSs to receive parcels through installed lockers. (3) In Taiwan, a company that provides the home delivery service for medications has existed. Hence, if one would like to implement other pickup facility delivery methods, the home delivery method must be included. Home delivery can support customers who live in places without CSs; and it can support the case that it is not worth investing in installing lockers when compared with the home deliver method. (4) The only skill that customers require to use the proposed service is to order medications through the Internet, that is, customers only need to upload their prescriptions and make the order through an online channel. Generally, most inhabitants in Taiwan can access the Internet, and are able to use the Internet for health care purposes. For instance, the outpatient department in some hospitals requires patients to make an online reservation in advance, and hence, they may be familiar with the online health care service.

In the future, this work can be applied to practical case studies with real data. We plan to consider more factors affecting the real-world vehicle routing situations, for example, vehicle type, alternative path selection, and traffic conditions. In addition, the customers’ satisfaction may be considered in the further framework to choose the delivery method which they are satisfied with.

Conclusion

This work has proposed the medication distribution problem with multiple delivery methods: CSs, lockers, and home delivery. Customers can receive medications through either CSs, lockers, or customers’ homes. The main contribution of this work is to additionally employ the 24-h CS pickup service serving customers to receive their medications at CSs. Although both CSs and lockers perform similar features as pickup facilities, they still have some different characteristics on the facility existence and the investment cost. This work further proposes a GA approach to solve the problem. The proposed GA can be used to evaluate the effects of the cost when the coverage distance of a pickup facility, the number of CSs, the contract fee of a CS chain, and the cost of installing a locker are changed. In experimental results, various scenarios are created and observed when the number of potential locker locations and the number of customers are adjusted whereas the number of CSs is fixed. Experimental results clearly indicate the effectiveness of the proposed GA.

Footnotes

Acknowledgements

The authors thank the anonymous referees for comments that improved the content as well as the presentation of this article.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported, in part, by the Ministry of Science and Technology, Taiwan, under grants MOST 106-2221-E-009-101-MY3 and MOST 108-2628-E-009-008-MY3.

ORCID iD

Chun-Cheng Lin

References

Winge

Johannesson

Perjons

, et al. The coordination hub: toward patient-centered and collaborative care processes. Health Informatics J 2015; 21(4): 284–305.

Hamann

Bezboruah

KC.

Outcomes of health information technology utilization in nursing homes: do implementation processes matter?

Health Inform J 2020; 26(3): 2249–2264.

Fathollahi Fard

Govindan

Hajiaghaei-Keshteli

, et al. A green home health care supply chain: new modified simulated annealing algorithms. J Clean Prod 2019; 240: 118200.

Watson

Ahmed

Mccall

, et al. Home delivery of medication—the role of a patient information leaflet on maximising service uptake. Int J STD AIDS 2014; 25(14): 1035–1037.

Anosike

Adibe

Isah

, et al. Willingness to pay for pharmacist-provided home telemonitoring among patients with chronic diseases in Enugu metropolis. Health Inform J 2020; 26(2): 829–840.

Veenstra

Roodbergen

Coelho

, et al. A simultaneous facility location and vehicle routing problem arising in health care logistics in the Netherlands. Eur J Oper Res 2018; 268(2): 703–715.

Zhou

Baldacci

Vigo

, et al. A multi-depot two-echelon vehicle routing problem with delivery options arising in the last mile distribution. Eur J Oper Res 2018; 265(2): 765–778.

Pham

Hà

Nguyen

XH.

Solving the multi-vehicle multi-covering tour problem. Comput Oper Res 2017; 88: 258–278.

Schneider

Drexl

A survey of the standard location-routing problem. Ann Oper Res 2017; 259(1–2): 389–414.

10.

Stenger

Schneider

Schwind

, et al. Location routing for small package shippers with subcontracting options. Int J Prod Econ 2012; 140(2): 702–712.

11.

Ghaffari Nasab

Ghazanfar Ahari

Ghazanfari

. A hybrid simulated annealing based heuristic for solving the location-routing problem with fuzzy demands. Sci Iran 2013; 20(3): 919–930.

12.

Zhang

Yang

, et al. Application of hierarchical facility location-routing problem with optimization of an underground logistic system: a case study in China. Math Probl Eng 2018; 2018: 1–10.

13.

Drexl

Schneider

A survey of variants and extensions of the location-routing problem. Eur J Oper Res 2015; 241(2): 283–308.

14.

Zhang

, et al. Sustainable multi-depot emergency facilities location-routing problem with uncertain information. Appl Math Comput 2018; 333: 506–520.

15.

Stenger

Vigo

Enz

, et al. An adaptive variable neighborhood search algorithm for a vehicle routing problem arising in small package shipping. Transport Sci 2013; 47(1): 64–80.

16.

Hof

Schneider

Goeke

Solving the battery swap station location-routing problem with capacitated electric vehicles using an AVNS algorithm for vehicle-routing problems with intermediate stops. Transport Res B-meth 2017; 97: 102–112.

17.

Zhou

Liu

X-F.

A novel hybrid genetic algorithm for the location routing problem with tight capacity constraints. Appl Soft Comput 2019; 85: 105760.

18.

Capelle

Cortés

Gendreau

, et al. A column generation approach for location-routing problems with pickup and delivery. Eur J Oper Res 2019; 272(1): 121–131.

19.

Hiassat

Diabat

Rahwan

A genetic algorithm approach for location-inventory-routing problem with perishable products. J Manuf Syst 2017; 42: 93–103.

20.

Fazayeli

Eydi

Kamalabadi

IN.

Location-routing problem in multimodal transportation network with time windows and fuzzy demands: presenting a two-part genetic algorithm. Comput Ind Eng 2018; 119: 233–246.

21.

Ahmadi Javid

Amiri

Meskar

. A profit-maximization location-routing-pricing problem: a branch-and-price algorithm. Eur J Oper Res 2018; 271(3): 866–881.

22.

Schiffer

Walther

The electric location routing problem with time windows and partial recharging. Eur J Oper Res 2017; 260(3): 995–1013.

23.

Cuda

Guastaroba

Speranza

MG.

A survey on two-echelon routing problems. Comput Oper Res 2015; 55: 185–199.

24.

Breunig

Schmid

Hartl

, et al. A large neighbourhood based heuristic for two-echelon routing problems. Comput Oper Res 2016; 76: 208–225.

25.

Pichka

Bajgiran

Petering

, et al. The two echelon open location routing problem: mathematical model and hybrid heuristic. Comput Ind Eng 2018; 121: 97–112.

26.

Zhu

Ursavas

Design and analysis of a satellite network with direct delivery in the pharmaceutical industry. Transp Res E Logist Transp Rev 2018; 116: 190–207.

27.

Marín

Martínez-Merino

Rodríguez-Chía

, et al. Multi-period stochastic covering location problems: modeling framework and solution approach. Eur J Oper Res 2018; 268(2): 432–449.

28.

Farahani

Asgari

Heidari

, et al. Covering problems in facility location: a review. Comput Ind Eng 2012; 62(1): 368–407.

29.

Zarandi

MHF

Davari

Sisakht

SAH.

The large-scale dynamic maximal covering location problem. Math Comput Model 2013; 57(3-4): 710–719.

30.

Karaoğlan

Erdoğan

Koç

Ç.

The multi-vehicle probabilistic covering tour problem. Eur J Oper Res 2018; 271(1): 278–287.

31.

Nedjati

Izbirak

Arkat

Bi-objective covering tour location routing problem with replenishment at intermediate depots: formulation and meta-heuristics. Comput Ind Eng 2017; 110: 191–206.

32.

Naji Azimi

Renaud

Ruiz

, et al. A covering tour approach to the location of satellite distribution centers to supply humanitarian aid. Eur J Oper Res 2012; 222(3): 596–605.

33.

Flores-Garza

Salazar-Aguilar

Ngueveu

, et al. The multi-vehicle cumulative covering tour problem. Ann Oper Res 2017; 258: 761–780.