Abstract

1. Introduction

2. Mathematical formulation

The problem is defined on a complete undirected graph G = (V₀, E) with node set V₀ and edge set E. Therein, V₀ = V ∪ {0}, where 0 denotes the depot, and V denotes the set of population nodes and/or potential distribution centers (DCs). In the envisaged application context, a population node is typically a village; therefore, for shortness, we shall use here the term “village” instead of “population node”.¹ A village i ∈ V that is not a potential DC can be viewed alternatively as a potential DC with capacity zero. Conversely, a potential DC that is not a village can be viewed alternatively as a village with zero inhabitants. Thus, in the graph representation, it is not necessary to make a fundamental distinction between villages and potential DCs. Nevertheless, we shall preferably use the index i for a node when referring to its role as a village, and the index j when referring to its role as a DC.

As usual in the vehicle routing literature, we extend the given transportation network to a complete graph by setting E = {(i, j)|i ∈ V₀, j ∈ V₀, i < j}. Note that edges are undirected. The distance between two nodes i ∈ V₀ and j ∈ V₀ is denoted by d_ij. The driving cost of a vehicle from node i to node j is assumed as proportional to d_ij, i.e., it can be expressed as τ · d_ij.

We restrict ourselves to a single commodity to be supplied to clients (i.e., inhabitants of the villages). This consideration also covers the situation where different commodities are to be supplied, but where their composition is homogeneous and there are no large items that could cause divisibility problems.

For the transportation of the commodity, a fleet K of vehicles is assumed to be available. The vehicles k ∈ K can have different capacities, but driving costs for a certain route are assumed to be identical, which represents a situation where velocities do not differ and the main cost component is the salary of the driver.² The (load) capacity of vehicle k is denoted by Q_k.

By the delivery period, we understand the time interval during which each (opened) DC has to be visited once. Depending on the application context, the length of a delivery period can be 1 day or also 1 week or even more. Let us start with the consideration of the planning process for a single delivery period; the case where periodic deliveries take place will be outlined in Remark 3 below.

The cost for opening a DC in node j for the current delivery period is denoted by c_j. In each node, at most one DC can be opened. Possible cost components can be costs for leasing a building, setup costs, fitment costs, surveillance costs or other. We always let opening costs refer to the delivery period under consideration. E.g., if at the beginning of the delivery period, no DC has yet been installed in node j, then the full costs for establishing a DC in this node arise if it is opened. If a DC in node j was already open in the last delivery period, then c_j only consists of the costs for keeping it operative in the current period. In this way, also the cost effect of changing the locations of opened DCs can easily be captured.

Furthermore, to each (potential) DC j, a capacity γ_j is assigned, which can depend on the location j. The quantity γ_j is to be interpreted as an upper bound for the throughput, i.e., for the amount of the commodity that can be supplied to clients at node j within one delivery period.

The demand occurring in village i ∈ V for the considered delivery period is denoted by W_i. The quantities W_i (i ∈ V) are assumed to be random variables with a known joint distribution. The random variables need not to be independent. It is convenient to decompose W_i as W_i = ξ_i · w_i, where w_i = E(W_i) is the expected value of W_i, and the random variables ξ_i characterize the effect of uncertainty. The numbers w_i can be obtained as baseline estimates, derived from the population sizes of the villages.

Covering tour models (based on facility location models using “coverage” measures) assume that the distance between a client and the nearest tour stop determines whether or not this client can be counted as “covered” by the chosen supply plan: if the distance is smaller or equal to some threshold value d₀, the client is covered, otherwise she/he is uncovered (cf. Current and Schilling [2], Gendreau et al. [6]). In the present paper, we introduce a flexible generalization of this basic concept by assuming that (i) each client considers the nearest DC for possibly supply, (ii) the share of clients actually traveling from their village i to a (nearest) DC in node j is given by a non-increasing function ψ of the distance d_ij. The special case where ψ(d) = 1 if d ≤ d₀ and ψ(d) = 0 otherwise leads us back to the assumptions of the traditional covering tour model. Possible alternatives are exponential decay functions of the form ψ(d) = r^d with 0 < r < 1, or step functions of the form

with 0 < d₀ < d₁ <  ⋯  < d_R and 1 > b₁ > b₂ >  ⋯  > b_R > 0.

The structure of the decision process is supposed to be the following:

• • Before the delivery period starts, a decision has to be made on which DCs to open, and on the tours by which the available vehicles visit the opened DCs. Split deliveries are not allowed, i.e., each opened DC has to be visited by exactly one vehicle. We call the decision on DCs and tours the first stage decision.
• • During the delivery period, actual demands become known. More precisely, the drivers can now observe the request in each opened DC j, where the request in DC j is defined as that part of the demand that is claimed by the people who have come to DC j. If the request in a DC j exceeds the capacity γ_j of the DC, the maximum supply in DC j is γ_j; otherwise, the maximum supply in this DC is identical to the request. If the load of a vehicle k (limited by the capacity Q_k) is sufficient to cover the maximum supplies in all DCs assigned to this vehicle, the driver delivers in each DC just this maximum supply. (The commodity is assumed as available in the depot to an unlimited amount.) If, on the other hand, the load of the vehicle is not sufficient, the driver delivers the entire load Q_k to the visited DCs, choosing an arbitrary distribution over the DCs. Since for our objective functions (which will be explained below), it does not matter the demand of which clients is satisfied, the special choice of the distribution is irrelevant. As a particular consequence, it makes no difference for our model whether the requests in each DC become known at once (e.g., communicated to the depot by phone), which gives the driver more freedom in choosing the distribution, or gradually (when going from DC to DC, the driver observes how many people are waiting there for delivery). Typically, in case of shortage, the driver will try to distribute the available load as evenly as possible for equity reasons, but this aspect is not represented in our objective functions. We call the decision on the amounts of actual deliveries to each DC on each tour the second stage decision.

Two objective functions are considered:

• • Objective function f₁ measures the total cost, i.e., the cost for opening DCs plus the total driving cost.
• • Objective function f₂ measures the expected uncovered demand. Uncovered demand is composed of (i) the demand of those clients who do not go to the nearest DC because the distance is too large for them, (ii) the demand that is not satisfied in a DC because of capacity limits of the DC, and (iii) the demand that is not satisfied in a DC because of vehicle capacity limits. The expectation is taken with respect to the outlined stochastic model.

We formulate the problem as a bi-objective, two-stage stochastic optimization problem with recourse. The following decision variables are used:

• • x_ijk denotes the number of times vehicle k drives through edge (i,j).
• • y_ij = 1 if village i is assigned to a DC in node j, and y_ij = 0 otherwise.
• • z_jk = 1 if vehicle k visits node j, and z_jk = 0 otherwise.
• • u_jk denotes the supply delivered by vehicle k to a DC in node j.

The variables x_ijk, y_ij and z_jk determine the first-stage decision, whereas the variables u_jk constitute the second-stage decision. We shall use the abbreviation δ(S) for the set of edges in E with exactly one end-node in S ⊆ V₀, i.e., δ(S) = {(i, j) ∈ E|i ∈ S, j ∈ V₀ ⧹ S or j ∈ S, i ∈ V₀ ⧹ S}. Instead of δ({j}), we shall shortly write δ(j). Moreover, the abbreviation x_k(E′) = ∑_(i,j)∈E′x_ijk for a set E′ ⊆ E will be used. The 3-dimensional array (x_ijk) and the matrices (y_ij), (z_jk) and (u_jk) will be denoted by x, y, z and u, respectively.

The bi-objective stochastic covering tour problem is now given as follows:

First stage:

Second stage:

Eq. (2) characterizes the first-stage problem as a bi-objective minimization problem with objective functions f₁ and f₂. Eq. (3) defines f₁ as the sum of driving costs and opening costs; observe that for j ∈ V, a DC in node j is opened if ∑_k∈Kz_jk = 1. Eq. (4) introduces f₂ as the expected uncovered demand, where the expectation is taken with respect to the distribution of the vector ξ of random variables ξ_i. The actual uncovered demand R(y, z, ξ) results as the solution of the second-stage problem (15)–(19). In the terminology of stochastic programming, we can call f₂ the recourse function; it expresses the expected cost of the decision to be made after the unknown parameter values become known. In our case, this decision consists of the choice of the amounts of the commodity to be unloaded in each DC, based on information on the actual requests. Since the first objective function is already determined by the solution (x,y,z) of the first-stage problem, the recourse function only captures the value of the second objective function, i.e., of the expected uncovered demand.

Constraints (5) ensure that to each village, exactly one DC is assigned. Constraints (6) express that a village can only be assigned to an opened DC. Constraints (7), wherein M denotes a sufficiently large number, make sure that each village is assigned to the nearest opened DC. Constraints (8) state that each DC can have no more than one visiting vehicle, and inherently there can only be at most one open DC in a village. Constraint (9) ensures that all vehicles visit the depot. Constraints (10) are the degree constraints for the tours. Constraints (11) are the usual subtour elimination constraints. Constraints (12) state that a link between two villages can be traversed by a vehicle at most once, whereas a link between a village and the depot can be traversed at most twice (the possible value 2 accounts for the feasibility of a forth-and-back tour between the depot and a single village). Constraints (13) and (14), finally, define the variables y_ij and z_jk as binary variables.

The first-stage problem combines the mathematical programming model for the CVRP with that of a (bi-objective) location problem: the location aspect is given by the choice of the DCs. For a fixed choice of the nodes in which DCs are to be opened, a CVRP remains to be solved. In our formulation, Eqs. (10) and (11) are similar to the usual connectivity and subtour elimination constraints in a CVRP, respectively, but they are extended in such a way that only opened DCs are allowed (and required) to be visited by vehicles.

In the second-stage problem, Eq. (15) defines the objective function (to be minimized) as the uncovered demand, where the uncovered demand results as the difference between total demand and total supply. Constraints (16) state that the supply u_jk delivered by vehicle k in DC j must not exceed the total request in DC j. Note that the total request in a DC is obtained as the sum of the requests from all villages that are assigned to this DC, and that the request from a village i addressed to DC j is given as the demand ξ_iw_i in village i, multiplied by the factor ψ(d_ij). Constraints (17) ensure that the supply delivered by a vehicle to a DC does not exceed the capacity of this DC, and that vehicle k only delivers to DC j if it stops there. Constraints (18) make sure that the total supply provided by a vehicle does not exceed its load capacity, and constraints (19) define the supplies as nonnegative reals.

3. Solution technique

4. Computational experiments

4.2. Examples

In order to provide some insight into the problem, we now take two of those instances as examples and study selected solutions. The first instance is called Mbayene and has 11 villages. Fig. 1 shows two solutions for this instance. These solutions are successive elements on the Pareto front. Each of the two solutions has two routes; each route visits exactly one village. Thick curvy lines represent these routes, following the real road network, and starting from the square representing the depot. (We do not represent the overall network in order to keep the pictures readable.) Full and dashed lines represent villagers coming from another village to a DC; full lines are for the case where 100% of the villagers come, dashed lines for the case where 50% only of the villagers come. If we look at the first solution, we observe that the village called Niandoul Ouolof is open as a DC, since it is visited by a vehicle. All villagers from Diemoul Ouolof go to this DC, as well as half of the villagers from Sine Amar and from Ndiourky. Now if we look at the second solution, Niandoul Ouolof is no longer open as DC; Diemoul Ouolof is open instead, and all villagers from Niandoul Ouolof as well as half of the villagers of Sine Amar and Ndiourky go to Diemoul Ouolof. In terms of cost, the second solution is worse: the routing cost is slightly higher, and the number of open DCs is the same. In terms of coverage, the second solution is better, although this is less obvious to see: it is due to the facts that (i) the capacity at Niandoul Ouolof is less than what would be needed to satisfy all inhabitants using it as a DC, and (ii) the capacity at Diemoul Ouolof is higher than that at Niandoul Ouolof because it is a larger village. In the end, the second solution has a slightly better coverage but also costs a little bit more, because the larger village is farther away from the depot.

Open in a separate window

Fig. 1

Example of two adjacent solutions on the Pareto front for instance Mbayene. (a) First solution to instance Mbayene and (b) Second solution to instance Mbayene.

The second example refers to the communauté of Meouane. Fig. 2 shows a solution close to the end of the Pareto front (corresponding to full coverage). For this reason, most villages in this solution are fully covered. Two of them are only partially covered though, and one of them is not covered at all. Moreover, it could be that some of the fully covered villages have in fact unsatisfied demand, due to capacity restrictions at the DCs.

Open in a separate window

Fig. 2

A solution close to the end of the Pareto front for instance Meouane.

4.3. Experimental results

The algorithms described in Section 3 are implemented in C++ and compiled with g++ 4.3.2. CPLEX 12.1 is used as MIP solver and multi-threading is disabled. The program is run on a 2.67 GHz Intel Xeon CPU using the Linux operating system.

4.3.1. Exact algorithm

For each instance, we run the algorithm described in Section 3 with a CPU time limit of 3 days. We now report, for each instance, the number of nodes, the total number of solutions found, and the CPU effort in seconds. If the solution method cannot find the whole Pareto front within the 3 day budget, we report “ > 3d” for CPU time. Table 1 summarizes these results.

Table 1

Problem size, Pareto front size and CPU effort for Senegal instances.

Instance	# nodes	Front size	CPU time
Cherif Lo	10	13	11.77
Diender Guedj	22	34	> 3d
Fandane	24	51	> 3d
Fissel	21	107	> 3d
Koul	16	29	569.11
Malicounda Wolof	12	27	36.08
Mbayene	12	31	81.75
Mboro	31	112	> 3d
Meouane	17	31	231.13
Merina Dakhar	19	46	2451.19
Mont Roland	15	17	128.31
Ndiagagniao	24	92	> 3d
Ndiakhene	9	6	1.42
Ndiass	18	55	6924.66
Ndieyene Shirak	11	27	49.21
Neugeniene	17	36	883.28
Ngandiouf	19	75	13261.8
Nguekhokh	19	36	573.29
Notto Gouye Diama	11	16	11.51
Notto	28	43	> 3d
Pekesse	12	24	54.72
Pire Goureye	18	55	3663.46
Pout	29	55	> 3d
Sandira	15	28	198.01
Tassette	21	33	4973.5
Thiadiaye	13	21	58.33
Thienaba	10	14	14.06
Thilmanka	14	39	268.81
Tiba Ndiaye	14	22	104.49
Touba Toul	20	43	12,706.5

Open in a separate window

There seems to be a tractability limit closely related to the instance size: instances with 20 nodes or less can all be solved, one instance with 21 nodes can also be solved (Tassette), but one instance with 21 nodes cannot be solved (Fissel) and all instances with more than 21 nodes cannot be solved either.

4.3.2. Heuristic algorithm

We also solved all instances of Table 1 by the heuristic algorithm H1 of Section 3.3. For this purpose, we applied a schedule with the following specification:

• • The first period lasts half of the total invested CPU time B and has a gap tolerance set to the default value of 0.0001 (this is the exact method);
• • the second period is split into 10 subperiods, each increasing the gap tolerance of the previous period by 10%.

In other words, we choose the parameters of H1 as 𝒯 = 11, c_t = 0.5B + 0.05(t − 1), and g_t = 0.0001 + 0.1(t − 1) (t=1,…,11). Using this specification and a total CPU budget B of 8 h, H1 was able to solve all instances. In all cases, H1 stopped before reaching the CPU budget limitation of 8 h: the highest gap tolerance ever reached in any of the runs was 70%.

Let us now look at a first example, the Diender Guedj instance. This instance could not be solved by the exact method, although the last found solution has a very low uncovered demand; so it is likely that only few solutions on the true Pareto front have been missed by the exact method after the invested runtime of 3 days. In Fig. 3, we plot in the same diagram the points on the Pareto front found by the exact method (red “+” signs) and the image points of the solutions proposed by the heuristic H1 (green “×” signs). A vertical blue line shows the end of the first period for the heuristic H1, i.e., the period where H1 works as the exact algorithm. The solutions corresponding to the points ×to the left of this line were found using the parameter setting for the exact method, whereas those corresponding to the points to the right of this line were found using an increased gap tolerance. (Since the CPU budgets were different for the exact method and for H1, the blue line has no interpretation for the execution of the exact algorithm.)

Open in a separate window

Fig. 3

Comparison of exact and H1 for instance Diender Guedj. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We can see that in this example, if H1 is given 8 h, it finds all solutions that the exact method has found within 3 days of CPU effort, plus an additional solution with zero uncovered demand.

As a second example, let us consider the Ndiagagniao instance (see Fig. 4). Although less spectacular, the comparison shows again that H1 produces high-quality solutions: up to the largest f₁ value for which the exact algorithm has found a solution, all heuristic solutions are Pareto-optimal or have at least image points very close to Pareto-optimal points. To the right hand side of the largest f₁ value for which a solution was still found by the exact approach, H1 delivers additional solutions that are non-dominated by any of the solutions obtained by the exact method. This example also shows, however, that as a consequence of the nonzero gap tolerance given to the MIP solver, it can occur that H1 proposes a solution that is even dominated by another solution proposed by H1: In the considered instance, H1 proposes a solution with cost 62,411 and uncovered demand 647.33, which is dominated by another proposed solution (found in a later iteration during the execution of H1) with cost 62,083 and uncovered demand 578.01. It is easy to see that this can happen although the epsilon constraint is always respected. Nevertheless, we do not consider this as an important issue in a heuristic context, as it is always possible to filter out dominated solutions.

Open in a separate window

Fig. 4

Comparison of exact and H1 for instance Ndiagagniao.

Figs. 3 and 4 also show that in addition to supported solutions, the Pareto front can contain a considerable number of unsupported solutions as well (cf. Section 3.2). In Fig. 3, e.g., solutions no. 4, 7, 9, 10 and 11 from left are obviously unsupported: None of these solutions can be optimal under any weighted average of the two objectives. In the right part of the Pareto front, more than half of all solutions are unsupported. Omitting unsupported solutions in Fig. 4 would leave a larger gap in the cost range between 14,000 and 18,000 and may thus hide meaningful solutions for a decision maker with a low budget.

5. Conclusions

We have introduced a bi-objective covering tour model with stochastic demand where the two objectives are given by (i) cost (opening cost for distribution centers plus routing cost for a fleet of vehicles) and (ii) expected uncovered demand, respectively. In our model, we take account of the fact that only a certain percentage of clients is willing to go from their homes to the nearest distribution center, depending on the distance. An application in humanitarian logistics has been outlined, but other application areas are also possible.

For the computational solution of the resulting bi-objective two-stage stochastic program with recourse, we use a branch-and-cut technique within an epsilon-constraint algorithm, applied to a sample-average version of the problem obtained from a fixed random sample of demand vectors. For this counterpart version, our algorithm provides the exact set of Pareto-optimal solutions. Although already the single-objective, deterministic boundary case of the considered problem is an NP-hard Covering Tour Problem and the complexity is further increased by the bi-objective and the stochastic generalization, we were able to compute solutions for real-world test instances from rural communities in Senegal up to instance sizes of about 20 villages per community.

Future research should address at least three questions: first, we work with a fixed sample of scenarios and have to keep the number of scenarios comparably low for reasons of computational efficiency. To obtain a better approximation of the sampled counterpart model to the real stochastic structure of the demand, it would be desirable to apply variance reduction techniques such as Importance Sampling, and to switch to a variable-sample approach which allows the consideration of a much larger set of random scenarios without compromising runtime too much. The Adaptive Pareto Sampling (APS) technique [7], a variable-sample method for multi-objective stochastic optimization, could be useful in this context.

Secondly, for tackling larger instances, the application of multi-objective heuristic techniques seems indispensable. Well-known multi-objective metaheuristics such as the Nondominated Sorting Genetic Algorithm II (NSGA-II) by Deb et al. [3] could be used in this context. This may again be based on either a fixed-sample or a variable-sample approach.

Third, in our model, rather simple assumptions on the choice of distribution centers by clients are made. These assumptions could be refined, e.g., by adopting concepts from consumer behavior in location models (cf. [4]).

Acknowledgments

Footnotes

References