Combining Heuristics with Simulation and Fuzzy Logic to Solve a Flexible-Size Location Routing Problem under Uncertainty

Abstract

The location routing problem integrates both a facility location and a vehicle routing problem. Each of these problems are NP-hard in nature, which justifies the use of heuristic-based algorithms when dealing with large-scale instances that need to be solved in reasonable computing times. This paper discusses a realistic variant of the problem that considers facilities of different sizes and two types of uncertainty conditions. In particular, we assume that some customers’ demands are stochastic, while others follow a fuzzy pattern. An iterated local search metaheuristic is integrated with simulation and fuzzy logic to solve the aforementioned problem, and a series of computational experiments are run to illustrate the potential of the proposed algorithm.

1. Introduction

When designing and managing supply chains, one of the most relevant problems is the simultaneous location of distribution facilities and the routing of vehicles to deliver products to a set of geographically dispersed customers. The former is considered a strategic decision, while the latter is operational. This problem is known in the scientific literature as the location routing problem (LRP). The LRP addresses these two types of decisions in an integrated manner. From the formal view of the operational research community, the LRP is known to be NP-hard, since it can be reduced to either the facility location problem (FLP), the vehicle routing problem (VRP) or the multidepot VRP, which are all known to be NP-hard. This computational complexity means that optimal solutions are really difficult to obtain in a reasonable computational time. Thus, heuristic approaches are required to solve medium- and large-sized instances. Due to its complexity, some of the first studies tackled the problem by splitting it into the corresponding subproblems [1,2]. Nevertheless, this approach might lead to suboptimal solutions.

Due to the increase in computational power and the development of fast heuristic approaches, the LRP has been studied in an integrated way, which clearly has improved the obtained results [3]. One of the most studied versions of the LRP is the capacitated LRP, in which both depot and vehicle capacity constraints must be satisfied (the acronym LRP will henceforth refer to this version). However, all previous works consider the depot capacity as a fixed value for each location. This could not be a suitable approach when dealing with realistic problems, since it is usual that decision-makers can select the size of a facility from a discrete set of known available sizes, or even freely. For real-world problems, this set is usually associated with investment activities, such as building facilities [4], purchasing equipment [5] or qualifying workforce [6]. From an academic point of view, despite the increasing number of published works on the LRP, the consideration of flexible sizes for facilities has been rarely addressed in the literature. Nevertheless, real-life examples from both LRP [4,7,8] and non-LRP [5,6] contexts show the relevance of considering a variety of facility sizes to select from.

Traditional LRP approaches consider that parameters are deterministic or crisp, i.e., they assume that inputs are known in advance. This assumption is far from reality in many applications, such as waste collection, humanitarian logistics and urban freight distribution, where uncertainty is a key factor to consider. Despite this, the literature on the LRP addressing uncertain parameters is still scarce. In order to overcome this problem, articles employing stochastic approaches can be found in the literature. Customers’ demand is one of the most addressed stochastic parameters [9,10,11,12,13]. Other parameters might also be considered as stochastic, such as transportation costs and travel speeds [14] or logistic costs and travel distance [15]. In general, many articles addressing stochasticity in routing problems hybridize simulation models with heuristic or metaheuristic algorithms to tackle efficiently both uncertainty and NP-hardness. In many real-life situations, however, it might not be possible to accurately model all uncertainty sources as stochastic variables following a probability distribution. This might be the case, for instance, when the volume of observations is low or the available data does not have enough quality [16]. Hence, uncertainty in the LRP has also been tackled through the use of fuzzy sets. Parameters such as customers’ demands [17,18,19,20], travel times [21,22] or time windows [23] have been modeled as fuzzy in several studies. Notice that, whenever possible, modeling uncertainty as stochastic variables might allow a deeper statistical analysis of the results.

To the best of our knowledge, there are no works in the literature simultaneously addressing stochastic and fuzzy approaches to model demand uncertainty in a flexible-size LRP. This is a realistic scenario, since many companies might have historical data on trustworthy customers and not enough data on new or unreliable ones. Hence, the main contributions of this paper are two-fold: on the one hand, a new variant of the location routing problem is studied, where facility sizing decisions and hybrid fuzzy-stochastic demands are simultaneously considered. On the other hand, this paper proposes a competitive solution approach based on the hybridization of a metaheuristic algorithm with both simulation and fuzzy logic, i.e., a so-called fuzzy simheuristic, to solve the aforementioned problem. Indeed, simheuristics have been traditionally proposed to deal with stochastic issues in hard combinatorial optimization problems [24]. However, their hybridization with fuzzy logic has been rarely studied.

The remainder of this paper is organized as follows: Section 2 describes previous works on the topic. Section 3 presents a description of our addressed problem. Section 4 explains the fuzzy simheuristic approach used to solve the problem. Section 5 describes a series of computational experiments. Section 6 analyzes our obtained results. Finally, Section 7 draws some conclusions and future research perspectives.

2. Literature Review

This section presents a summary of the published manuscripts on the main topics addressed by this work. Thus, Section 2.1 outlines works related to the location routing problem in both its deterministic version and the variant including uncertain parameters. Additionally, Section 2.2 summarizes the main contributions on the field of simheuristics and fuzzy logic as methodologies to handle uncertain parameters in routing problems.

2.1. The Location Routing Problem

Perhaps the first work related to the location routing problem is the one by Maranzana [25], who analyze the influence of transportation costs on location decisions. Moreover, Salhi and Rand [1] quantify for the first time the benefits of considering routing decisions when locating facilities. They also state that solving each subproblem (location and routing) independently does not provide optimal solutions. Multiple variants of the LRP have been proposed over time. These variants depend on the characteristics of depots (capacitated or not), vehicles (capacitated or not, homogeneous or heterogeneous fleet), costs (symmetric or asymmetric) or the consideration of uncertain parameters. All capacitated variants addressed by these authors assume that depot sizes are fixed and cannot be changed.

Considering the limited computational power available at that time, the initial works on the LRP firstly solved the underlying location problem and used the obtained solution as a starting point to handle the corresponding routing problem. However, as the computational power has notably increased in recent years, the newest approaches deal with the LRP in an integrated manner [26,27]. Among the recently published works on the deterministic LRP, Escobar et al. [28] propose a granular tabu search within a VNS framework to speed up computational times without decreasing the solutions quality. A biased-randomization-based metaheuristic of two phases is developed by Quintero-Araujo et al. [27] to solve the capacitated version of the problem. Ferdi and Layeb [29] propose a GRASP with a novel technique used to create clusters around the open depots. Traditional applications of the LRP include horizontal cooperation [30], electric vehicle routing problems [31,32,33], city logistics [34], humanitarian logistics [35] or supply chain network design [36]. Moreover, most recent applications are related to environmental issues [37], cold supply chains [38] or waste management [39].

When dealing with uncertainty, most works have focused on the use of stochastic modeling. One of the utilized approaches has been the hybridization of simulation techniques with metaheuristics. For instance, Quintero-Araujo et al. [9] propose a simheuristic to solve an LRP with stochastic demands, by hybridizing Monte Carlo simulation with an iterated local search metaheuristic. A similar approach is employed by Tordecilla et al. [13], who address an LRP where the sizes of facilities to locate are also a variable to consider. Rabbani et al. [10] also propose a simheuristic approach that combines a nondominated sorting genetic algorithm-II (NSGA-II) and Monte Carlo simulation. They tackle a multiobjective multiperiod LRP in the context of the hazardous waste management industry. Both generated waste and number of people at risk are stochastic. Inventory decisions are also taken into account. Sun et al. [11] address a real-world case from an express delivery company in Shanghai. These authors tackle an LRP in which demand is stochastic and can be split for self-pickup. Then a simulation-based optimization model is proposed and two heuristics results are compared.

Other parameters are also considered to be uncertain. For instance, Herazo-Padilla et al. [14] hybridize an ant colony optimization metaheuristic with discrete-event simulation to solve an LRP in which both transportation costs and vehicle travel speeds are considered stochastic. Authors demonstrate that their proposed approach is not only efficient but is able to find statistical interactions among the different parameters. Zhang et al. [15] present an approach that hybridizes a genetic algorithm with simulation to solve a sustainable multiobjective LRP in the context of emergency logistics. The authors consider the travel distance, the demand and the cost of opening a depot as uncertain variables. Additionally, the emergence of new technologies introduces new challenges. This is the case of Zhang et al. [12], who address the problem of locating battery swap stations and routing electric vehicles with stochastic demands. This problem is solved using a hybrid approach that combines a variable neighborhood search with a binary particle swarm optimization algorithm. The problem’s complexity increases when considering the low autonomy of this type of vehicles, since route failures can frequently be present when demands are not known in an accurate manner.

Uncertainty in the LRP has been studied using either stochastic or fuzzy parameters.

Table 1. Recent works related to the location routing problem with uncertain parameters.

Reference	Type of Uncertainty	Uncertain Parameter	Mathematical Approach	Solving Approach	Objective Criterion
Quintero-Araujo et al. [9]	Stochastic	Demand	Mixed-integer linear programming	Simheuristic Iterated local search Monte Carlo simulation	Minimize cost
Rabbani et al. [10]	Stochastic	Demand Number of people at risk	Mixed-integer non- linear programming	Simheuristic NSGA-II Monte Carlo Simulation	Minimize cost Minimize environmental risk
Sun et al. [11]	Stochastic	Demand	Mixed-integer linear programming	Biogeography-based optimization Adaptive large neighborhood search	Minimize cost
Zhang et al. [12]	Stochastic	Demand	—	Variable neighborhood search Particle swarm optimization	Minimize cost
Tordecilla et al. [13]	Stochastic	Demand	—	Simheuristic Iterated local search Monte Carlo simulation	Minimize cost
Herazo-Padilla et al. [14]	Stochastic	Transportation cost Travel speed	Mixed-integer linear programming	Ant colony optimization Discrete-event simulation	Minimize cost
Zhang et al. [15]	Stochastic	Demand Travel distance Depot opening cost	Mixed-integer non- linear programming	Genetic algorithm Uncertain simulation	Minimize travel time Minimize emergency relief cost Minimize CO${}_2$ emissions
Zhang et al. [17]	Fuzzy	Demand	A fuzzy chance constrained model	Particle swarm optimization Variable neighborhood search Stochastic simulation	Minimize cost Minimize additional travel distance due to route failures
Mehrjerdi and Nadizadeh [18]	Fuzzy	Demand	A fuzzy chance constrained model	A greedy clustering method Ant colony system Stochastic simulation	Minimize cost
Fazayeli et al. [19]	Fuzzy	Demand	Mixed-integer non- linear programming	Exact approach Genetic algorithm	Minimize cost
Nadizadeh and Kafash [20]	Fuzzy	Demand	A fuzzy chance constrained model	A greedy clustering method Ant colony system Stochastic simulation	Minimize cost
Zarandi et al. [21]	Fuzzy	Travel time	A fuzzy chance constrained model	Simulated annealing Fuzzy simulation	Minimize cost
Zarandi et al. [22]	Fuzzy	Demand Travel time	A fuzzy chance constrained model	Simulated annealing Fuzzy simulation	Minimize cost Minimize additional travel distance due to route failures
Ghezavati and Morakabatchian [23]	Fuzzy	Time windows	Mixed-integer linear programming	Exact approach	Minimize cost Minimize risks
Ghaffari-Nasab et al. [40]	Fuzzy	Demand	A fuzzy chance constrained model	Simulated annealing Stochastic simulation	Minimize cost Minimize additional travel distance due to route failures
Nadizadeh and Nasab [41]	Fuzzy	Demand	A fuzzy chance constrained model	A hybrid heuristic algorithm Ant colony system Stochastic simulation	Minimize cost Minimize additional travel distance due to route failures
Wei et al. [42]	Fuzzy	Transportation cost Number of people that may be at risk	A fuzzy chance constrained model	Genetic algorithm Fuzzy simulation	Minimize cost Minimize risks

shows an overview of works addressing this topic, which includes: (i) whether the uncertainty is addressed stochastically or in a fuzzy fashion; (ii) the considered uncertain parameter; (iii) the mathematical modeling approach; (iv) the approach used to solve the problem; and (v) the objective function. Analyzed works show a clear preference for considering an uncertain demand, as well as for using fuzzy chance constrained models. Given both the considered uncertainty and the combinatorial nature of the LRP, most works employ a hybrid approach combining simulation with a metaheuristic algorithm. Finally, cost minimization is the prevalent objective, although a few works also consider the minimization of risk or the minimization of the additional travel distance due to route failures. Regarding works on fuzzy parameters, Zhang et al. [17] propose a hybrid particle swarm optimization (PSO) algorithm to solve a capacitated LRP with fuzzy triangular demands (CLRP-FD). The hybrid PSO algorithm is composed of three phases including a local search method and stochastic simulation. In addition, the authors propose a chance-constrained programming model for the CLRP-FD. Zarandi et al. [21] consider a multidepot LRP with fixed depot capacity and fuzzy travel times. Mehrjerdi and Nadizadeh [18] present a fuzzy chance constrained programming model where demands are modeled as fuzzy numbers. A four-phase method called “greedy clustering” is proposed, in which both an ant colony system metaheuristic and stochastic simulation are included. Fazayeli et al. [19] propose an LRP with time windows and fuzzy demands as the delivery part of a multimodal transport network. The mixed integer mathematical fuzzy model is coded and solved using GAMS and compared to the results provided by a genetic algorithm. Nadizadeh and Kafash [20] analyze a LRP with simultaneous pick-up and delivery in the context of reverse logistics. Both types of demands (pick-up and deliveries) are fuzzy variables. A fuzzy chance constrained programming model is proposed to represent the problem, and a greedy clustering method is used to solve it.

The analyzed works show that uncertainty in the LRP has been addressed either by using stochastic or fuzzy demands but never considering both types of uncertainty at the same time—e.g., that some customers’ demands are modeled as stochastic variables while others are modeled as fuzzy values. In addition, to the best of our knowledge, there are no previous studies on the LRP with facility sizing decisions and hybrid fuzzy-stochastic demands. Only Tordecilla et al. [13] have studied a similar LRP variant, although considering all customers’ demands as stochastic. Thus, our work aims to fulfill the existing gap by considering a flexible-size LRP and two different types of uncertain parameters: stochastic and fuzzy demands.

2.2. Simheuristics and Fuzzy Logic for Vehicle Routing Problems under Uncertainty

When dealing with combinatorial optimization problems subject to uncertain parameters, one of the most recommended approaches is the combination of simulation (to handle stochasticity) with heuristic-based methods (to deal with the optimization part of the problem) [43,44]. In that sense, a simheuristic approach is a relatively new and efficient technique to tackle combinatorial optimization problems under uncertainty [24,45]. In general, a simheuristic algorithm works as follows: (i) given a stochastic problem, the random variables are transformed into their deterministic counterpart by using expected values; (ii) an approximated framework (heuristic or metaheuristic) is used to generate high-quality solutions for the transformed deterministic instance that can also be “promising” solutions for the stochastic version of the problem; (iii) these promising solutions are sent to a simulation engine in order to estimate its quality in a stochastic environment. The simulation engine, in addition, provides feedback to better guide the search used by the approximated procedure; and (iv) an improved estimation of the quality of the solutions is obtained for a subset of “elite” solutions using a longer simulation process. Different simheuristic algorithms have been presented in the literature to solve routing problems.

Stochastic demands in vehicle routing problems are addressed by Quintero-Araujo et al. [46] and Gruler et al. [47]. Moreover, stochastic demands are also studied in arc routing problems [48]. Stochastic versions of the inventory routing problem can be found in Gruler et al. [49]. Real applications like the waste collection problem with stochastic demands are analyzed in Gruler et al. [50]. Intermodal routing problems have also considered other stochastic parameters, such as capacity [51] or travel times [52,53]. Additionally, the need of using fuzzy logic in vehicle routing problems arises when there are some vague or uncertain parameters. The literature presents various works in which fuzzy logic is added, for instance, to model uncertain demands [54,55,56,57], travel times [58,59], capacity [57,59], and service times [60]. Additional aspects are also considered by these works, such as time windows [57,61,62], environmental aspects [59], multiple objectives [59], intermodal transportation [57,59] and an open VRP [63]. Additional applications of metaheuristics combined either with Monte Carlo simulation or fuzzy logic can be found in several fields, such as scheduling [64,65], controller optimization [66,67] parameter estimation [68], finance [69], facility location [70], etc.

3. Problem Description

The location routing problem is a well-known problem in which three main decisions must be made: (i) locating one or more facilities; (ii) allocating customers to open facilities without exceeding their capacity; and (iii) designing a number of routes whose aggregated customers’ demand does not exceed a vehicle capacity. Each route must start and finish at the same facility. Furthermore, we consider a location routing problem with facility sizing decisions, where the size of each open facility is also a variable to decide on. Furthermore, we also consider both stochastic and fuzzy demands. Hence, a percentage of the vehicles’ capacity is reserved as a safety stock (SS), in case the demand is higher than expected. Therefore, the main decision variables in this problem are related to the number of facilities to open, the facilities’ size and location, which customers must be allocated to each open facility, how many vehicles must be used and how to design the associated routes. This problem is NP-hard since it contains, as special cases, the capacitated vehicle routing problem, the multidepot VRP and the facility location problem, all of them known to be computationally hard. Figure 1 provides an example of a complete LRP solution. Facilities are represented by diamonds and customers by circles. Black (solid) diamonds are the open facilities, while noncolored diamonds correspond to nonopen facilities. For each open depot, a set of routes starting and finishing at the corresponding depot location is designed to serve all customers’ demands. Each route is assigned a single vehicle.

Formally speaking, the LRP can be defined on a complete, weighted, and undirected graph $G(V,E,C)$, in which $V=J\cup I$ is the set of nodes (comprising the subset J of potential facility locations and subset I of customers), E is the set of edges, and C is the cost matrix of traversing each edge. Delivery routes are performed by a set K of unlimited homogeneous vehicles with limited capacity. This problem also assumes that all vehicles are shared by all facilities (i.e., no depot has a specific fleet) and each edge $e\in E$ satisfies the triangle inequality. The customers’ demands are uncertain and are modeled using stochastic values for a subset of customers $I_1$, and fuzzy values for a subset of customers $I_2$, such that $I_1\cup I_2=I$. The variant of the LRP considered in this paper is the one in which a decision must be made about the size of the facilities to open. Hence, a set L of alternative sizes for each facility and associated fixed and variable opening costs are provided as inputs. Depots might have equal or different capacities. Each customer node must be served by exactly one vehicle that starts and finishes its route in the facility to which it has been allocated (i.e., split deliveries are not allowed). The following notation is used to describe our problem:

Parameters

$s_l$ = Available size of type $l\in L$

$D_i$ = Uncertain demand of customer $i\in I$

$f_j$ = Fixed opening cost of depot $j\in J$

$o_{jl}$ = Variable opening cost of depot $j\in J$ with size of type $l\in L$

$c_e$ = Cost of traversing arc $e\in E$

q = Capacity of each vehicle

$\%SS$ = Safety stock percentage

Decision variables

$y_{jl}$ = Binary variable that indicates whether the depot $j\in J$ is open with size $l\in L$ or not.

$x_{ij}$ = Binary variable that indicates whether customer $i\in I$ is assigned to the depot $j\in J$ or not

$w_{ek}$ = Binary variable that indicates whether arc $e\in E$ is used in the route performed by vehicle $k\in K$ or not

The objective is to minimize the total cost ($TC$), which includes opening facilities costs ($OC$), routing costs ($RC$), and failure costs ($FC$), i.e., $TC=OC+RC+FC$. These parts are defined in Equations (1)–(3).

$$OC=\sum _{j\in J}\sum _{l\in L}(f_j+o_{jl})y_{jl}$$

(1)

$$RC=\sum _{e\in E}\sum _{k\in K}{c}_e{w}_{ek}$$

(2)

$$FC=min\{c_{reac},c_{prev}\}$$

(3)

$FC$ represents the cost incurred whenever the actual demand of a route is greater than the vehicle capacity, where $c_{reac}$ and $c_{prev}$ depend on the corrective action considered, namely:

• A reactive strategy with a cost $c_{reac}$, in which a vehicle must perform a round-trip to its assigned facility for a replenishment if the actual current-customer demand is higher than the vehicle’s current load.
• A preventive strategy with a cost $c_{prev}$, in which a vehicle must perform a detour to the facility before visiting the next customer. The decision about performing this detour depends on the type of demand of the next customer. If the demand is stochastic, the detour is carried out whenever the expected demand of the next customer is higher than the current capacity of the vehicle. Alternatively, if the demand is fuzzy, this decision depends on the comparison between the fuzzy values of both the demand of the next customer and the current capacity.

Let $\varnothing \ne S\subset V$ be a subset of nodes, ${\delta }^{+}\left(S\right)$ the set of edges leaving S, ${\delta }^{-}\left(S\right)$ the set of edges entering S, and $A\left(S\right)$ the set of edges with both ends in S. Hence, the location routing problem with facility sizing decisions and uncertain demands can be modeled as the following integer program:

$$MinimizeTC$$

(4)

                              subject to:

$$\sum _{k\in K}\sum _{e\in {\delta }^{-}\left(i\right)}{w}_{ek}=1\forall i\in I$$

(5)

$$\sum _{i\in I}\sum _{e\in {\delta }^{-}\left(i\right)}{D}_i{w}_{ek}\le (1-\%SS)q\forall k\in K$$

(6)

$$\sum _{e\in {\delta }^{+}\left(n\right)}{w}_{ek}=\sum _{e\in {\delta }^{-}\left(n\right)}{w}_{ek},\forall k\in K,\forall n\in V$$

(7)

$$\sum _{e\in {\delta }^{+}\left(J\right)}{w}_{ek}\le 1\forall k\in K$$

(8)

$$\sum _{e\in A\left(S\right)}{w}_{ek}\le \left|S\right|-1\forall S\subseteq I,\forall k\in K$$

(9)

$$\sum _{e\in {\delta }^{+}\left(j\right)}{w}_{ek}+\sum _{e\in {\delta }^{-}\left(i\right)}{w}_{ek}\le 1+x_{ij}\forall i\in I,\forall j\in J,\forall k\in K$$

(10)

$$\sum _{j\in J}{x}_{ij}=1\forall i\in I$$

(11)

$$\sum _{i\in I}{D}_i{x}_{ij}\le \sum _{l\in L}{s}_l{y}_{jl}\forall j\in J$$

(12)

$$\sum _{l\in L}{y}_{jl}\le 1\forall j\in J$$

(13)

$$\forall y_{jl},x_{ij},w_{ek}\in \{0,1\}$$

(14)

The objective function (4) minimizes the total cost. Constraint (5) ensures that each customer is served by a single route and a single vehicle. Constraint (6) guarantees that the total demand served by a vehicle in a route does not exceed its capacity. This limit is reduced by a safety stock, which is a percentage of the vehicle capacity reserved to respond more effectively to the uncertain demand. Constraint (7) guarantees the continuity of each route. Constraint (8) ensures the return of each vehicle to its starting depot. Constraint (9) guarantees the subtour elimination. Constraint (10) ensures that a customer is served by a route departing from an open depot only if this customer is allocated to this depot. Constraint (11) guarantees that a customer is assigned to only one depot. Constraint (12) ensures that the total demand served from a depot does not exceed its assigned size. Constraint (13) guarantees that only one size is assigned to an open depot. Finally, Constraint (14) determines that all decision variables are binary.

4. Solution Approach

Since the problem described in Section 3 is known for being NP-hard, the formulated mathematical model is not employed to find an optimal solution but just to provide a better understanding of the problem details Hence, we propose a fuzzy simheuristic approach [24] for minimizing the expected total cost. Traditionally, simheuristics have been used to solve optimization problems with stochastic components, such as arc routing problems with stochastic demands [48], stochastic waste collection problems [50] or team orienteering problems with stochastic travel times [71]. We have extended the simheuristic framework by including fuzzy components in order to deal with combinatorial optimization problems with uncertainty components of both stochastic and nonstochastic nature. In particular, our methodology combines an iterated local search (ILS) metaheuristic with Monte Carlo simulation and fuzzy inference systems (FIS) to deal with stochastic and fuzzy variables, respectively. As discussed in Ferone et al. [72], several metaheuristic frameworks offer a well-balanced combination of efficiency and relative simplicity and can be easily extended to a fuzzy simheuristic. In general, our approach is composed of three stages. During the first stage, a set of promising LRP solutions are generated using a constructive heuristic, which employs biased-randomization techniques [73]. In the second stage, the ILS metaheuristic tries to improve each of these promising solutions by iteratively exploring the search space and conducting a short number of simulations. Finally, in the third stage, a refinement procedure using a larger number of simulation runs is applied to these elite solutions, which allows one to obtain a more accurate estimation of the expected total cost.

Algorithm 1 outlines the main components of Stage 1. It generates quickly a ranked list of “promising” LRP solutions. The main input parameters of this heuristic are: the list of customers with both their demand and location in Cartesian coordinates, the list of facilities including their opening costs and the vehicle capacity. The algorithm procedure is as follows: initially, the minimum and maximum ($nbDepot{s}_0$ and $maxNbDepots$, respectively) numbers of facilities required to serve the total demand are computed. Both bounds are calculated by dividing the total demand by the maximum available facility size, and the minimum available facility size, respectively, and they are rounded up to the next integer number. Then we run our algorithm for each number of facilities between $nbDepot{s}_0$ and $maxNbDepots$ (line 3). Later, for each iteration of the line 4 loop, a new set of random locations are generated (line 5). This is stored in $usedOpenDepots$ to avoid repeating. Next, if the available capacity of facilities in $openDepots$ is enough to satisfy customers demand, customers’ allocation and routing procedures are carried out; otherwise, $openDepots$ is rejected. The customers’ allocation procedure is performed by producing a new $map$ (line 9) where each facility has a list of all customers sorted by savings. These savings represent the benefit of allocating each customer to the current depot instead to the best alternative facility. Then a facility in $openDepots$ is selected randomly, and a biased-randomized procedure is used to allocate a customer of the list to the current depot. This procedure ends when all customers have been allocated. In the step in line 10 a VRP is solved for each subset facility-customers in the map. Finally, a feasible LRP solution is yielded and stored in the pool of solutions $poolSol$. The algorithm ends returning a top list of complete LRP solutions, assessed in terms of opening and routing costs.

Algorithm 1 Constructive heuristic ($cust$, $depots$, $vehCap$, $\beta$, $ite{r}_{max}$)

1: $usedOpenDepots$←∅ 2: $\langle nbDepot{s}_0$, $maxNbDepots\rangle$←$computeDepotsBound(depots$) 3: for$nbDepots$←$nbDepot{s}_0$ to $maxNbDepots$ do 4:       for $iter$← 1 to $ite{r}_{max}$ do 5:             $openDepots$←$depotsToOpen\left(nbDepots\right)$ 6:             if $openDepots$∉$usedOpenDepots$ then 7:                 if $capacity\left(openDepots\right)$≥$demand\left(cust\right)$ then 8:                     $usedOpenDepots$←$add(usedOpenDepots,openDepots)$ 9:                     $map$←$allocateCustomers(openDepots,cust)$ 10:                   $lrpSol$←$CWS(map,\beta ,vehCap)$ 11:                   $poolSol$←$add(poolSol,lrpSol)$ 12:               end if 13:            end if 14:        end for 15: end for 16: return$sortingByCost\left(poolSol\right)$

Algorithm 2 outlines Stages 2 and 3. During the second stage, each “promising” map generated by the constructive heuristic is processed by the simulation and the fuzzy components to estimate its safety stock (line 4). This procedure is carried out by performing a low number of runs, where a new value is assigned to each random or fuzzy element based on its probability distribution or fuzzy function, respectively. We use Monte Carlo simulation in order to estimate the stochastic variables, whilst a fuzzy inference system is used to estimate the fuzzy variables. Then, the objective function and the constraints are evaluated under the random/fuzzy generated values to compute the expected cost of each promising map. Next, the ILS metaheuristic tries to improve the set of “promising” maps by iteratively exploring the search space and conducting a second process of fuzzy/simulation runs. We start the process by perturbing the current base solution baseSol (line 8). In this phase we use two different strategies. In the first one, the algorithm randomly selects a set of customers and tries to reassign them in a random way to another facility without violating its capacity. Regarding the second strategy, the algorithm randomly exchanges the allocation of a percentage of customers among facilities. This process is dependent on the value of k, which represents the degree of exchange to be applied. This value is updated in each iteration between $K_{min}$ and $K_{max}$, i.e., it is reset to $K_{min}$ whenever a new solution newSol outperforms the baseSol, and it is increased whenever the algorithm fails to improve the current solution until a maximum value $K_{max}$. The strategy to be used in each iteration of the algorithm is randomly selected.

Algorithm 2 ILS-based Fuzzy Simheuristic ($inputs$, $\alpha$, $\beta$, $\lambda$, $Inc$, $T_0$, $K_{min}$, $K_{max}$, $I_0$, $t_{max}$)

1: $initSol$← genInitSol(inputs,$\alpha$, $\beta$ ) 2: $baseSol$←$initSol$ 3: $bestSol$←$baseSol$ 4: fastSimulation(baseSol)% Fuzzy and Monte Carlo Simulation 5: T ←$T_0$ 6: while (time ≤$t_{max}$) do % ILS stage 7:     k ←$K_{min}$ 8:     perturbationSol ← perturbation(baseSol, k, $\alpha$, $\beta$) 9:     newSol ← localSearch(perturbationSol) 10:    if (detCost(newSol) < detCost(baseSol)) then 11:         fastSimulation(newSol) % Fuzzy and Monte Carlo simulation 12:         if (expCost(newSol) < expCost(baseSol)) then 13:              baseSol ← newSol 14:              if (expCost(newSol) < expCost(bestSol)) then 15:                  bestSol ← newSol 16:                  insert(poolBestSol,bestSol) 17:              end if 18:              k ←$K_{min}$ 19:           end if 20:       else% SA-based acceptance criterion 21:           temperature ← updateTemperature(detCost(newSol), detCost(baseSol), T) 22:           if ($\mathcal{U}$(0,1) ≤ temperature) then 23:                baseSol ← newSol 24:                k ←$K_{min}$ 25:           else 26:                k ← min(k * Inc,$K_{max}$) 27:           end if 28:       end if 29:       T ←$\lambda$T 30: end while 31: for (sol ∈ poolBestSol) do % Refinement stage - Fuzzy and Monte Carlo simulation 32:       longSimulation(sol) 33:       if (expCost(sol) < expCost(bestSol)) then 34:           bestSol ← sol 35:       end if 36: end for 37: return bestSol

Afterwards, the algorithm starts a local search around the perturbed solution in order to improve it (line 9). This stage consists in a two-opt inter-route operator, which interchanges two chains of randomly selected customers between different facilities. A newSol is returned whenever no more improvements are achieved. Later, whenever the deterministic cost of the baseSol is improved (line 10), the newSol is processed by the simulation and the fuzzy components to deal with the uncertainty of the proposed problem, using a low number of runs to compute the expected cost of the solution (line 11). Notice that this procedure does not only provide estimated values to the expected cost associated with the solutions generated by our approach, but it also reports feedback to the metaheuristic search process. If the newSol is also able to improve the expected cost of the baseSol (line 12), the latter is updated. In the same way, if the expected cost of the newSol improves the cost of the best solution (bestSol) found so far (line 14), the latter is updated and added to the pool of elite solutions (line 16). This pool contains the best stochastic/fuzzy solutions found so far. The number of solutions in this pool is a known parameter that depends on the available computational time. Moreover, by limiting the size of this pool we ensure that we only keep track of the top solutions as the algorithm evolves. In order to further diversify the search, the algorithm might occasionally accept nonimproving solutions following an acceptance criterion (lines 20–28). Specifically, we have used a simulated-annealing acceptance criterion, which contains a decaying probability that is regulated by a dynamic temperature parameter (T).

Finally, a refinement procedure using a larger number of simulation runs is executed in the third stage for each elite solution (lines 31–36). Hence, a more accurate summary of output variables can be obtained. As before, both probability distributions and fuzzy functions are employed in this simulation, depending on whether the element has a stochastic or fuzzy nature. Finally, the "best" solution (or pull of best alternative solutions) is returned, considering that the decision maker might be not only interested in the average value associated with a solution but also in its variability level. Particularly, the main output variables in our experiments are: the opening and routing costs, the cost incurred whenever a route fails and the safety stock.

5. Computational Experiments

Multiple sets of instances are found in the literature to test the algorithms designed to solve the LRP [74,75,76]. Nevertheless, these sets do not consider characteristics such as parameters uncertainty and flexible facility sizes, i.e., instances must be adapted to our problem’s features. Therefore, we use Akca’s [74] instances and introduce the following modifications:

• Traditional LRP instances consider that a single fixed size is available to assign to open depots. We extend this unit set to five alternative sizes, so that our algorithm selects one of them for each open depot. If $s_j$ is the size proposed by the original instance for each potential depot $j\in J$, and L is the set of available sizes, our approach’ alternative sizes are $s_{jl}\in \{(1-2r)s_j,(1-r)s_j,s_j,(1+r)s_j,(1+2r)s_j\}$, where $l\in L$, $0.0<r<0.5$, and r is the range of difference between available sizes. When $r=0$, the case is the same as the traditional LRP. We consider that $r=0.25$.
• Traditional LRP instances consider a fixed cost ($f_j$) incurred whenever a depot $j\in J$ is open. We keep this parameter unaltered. Additionally, we introduce a variable cost ($o_{jl}$) depending on $f_j$ and $s_{jl}$, namely: $o_{jl}={\displaystyle \frac{(s_{jl}-s_j)}{2s_j}}{\displaystyle \frac{{\sum }_j{f}_j}{\left|J\right|}}$. This formula preserves $o_{jl}$ in the same order as $f_j$ for each depot $j\in J$. Besides, it yields negative costs whenever $s_{jl}<s_j$, positive costs whenever $s_{jl}>s_j$, and a null cost when $s_{jl}=s_j$. Thus our results can be compared with those found in the LRP literature.
• An uncertain demand $D_i$ for each customer $i\in I$ is considered. The demand of half of the customers is assumed to follow a log-normal probability distribution. If ${\varphi }_i$ is the deterministic demand in the Akca’s set, then $E\left[D_i\right]={\varphi }_i$. In addition, three different values of variance are considered: low, medium and high, i.e., for $\lambda \in \{0.05,0.10,0.20\}$, $Var\left[D_i\right]=\lambda {\varphi }_i$. These variability values are preserved identical to the ones used by Tordecilla et al. [13], in order to perform a suitable results comparison. The demand of the other half of the customers is considered to be fuzzy. In this case, $D_i$ can be estimated as low (DL), medium (DM) or high (DH). The demand in each of these fuzzy sets is represented by a triangular fuzzy number $D_i=(d_{1i},d_{2i},d_{3i})$. If q is the vehicle total load capacity, all fuzzy demand values are expressed as a proportion of q in order to perform an appropriate comparison between the demand and the vehicle available capacity, i.e., $0\le D_i\le 1$. The membership function of these fuzzy sets are displayed in Figure 2.

Fuzzy Approach for the Demand and the Vehicle Available Capacity

When considering customers with stochastic demands, the decision about visiting the next customer in a route is made simply by comparing its expected demand with the vehicle’s current capacity. If this demand is greater, the vehicle will perform a detour to the depot for a replenishment. Nevertheless, when the next customer demand is fuzzy, the decision about serving it is made employing a preference index $p_i$ [77]. It indicates the strength of our inclination to visit the next node in a route. This index depends on both the estimated demand of the next node $D_{i+1}$ and the vehicle capacity $C_i$ that remains available after serving the customer $i\in I$. $C_i$ is expressed as a proportion of q, i.e., $0\le C_i\le 1$. It also can be treated as low (CL), medium (CM) or high (CH), and it is represented by a triangular fuzzy number $C_i=(c_{1i},c_{2i},c_{3i})$. The membership function of the capacity fuzzy sets are displayed in Figure 3.

The preference index is defined between 0 and 1, i.e., $0\le p_i\le 1$. When $p_i=1$, we will definitely visit the next node in a route since the vehicle available capacity can for sure meet its demand. When $p_i=0$, we are sure that $D_{i+1}$ exceeds $C_i$ and the vehicle must return to the depot for a replenishment. We consider that the preference can be very low (PVL), low (PL), medium (PM), high (PH) or very high (PVH). Each of these categories is represented by a fuzzy set, whose membership function is depicted in Figure 4. Additionally, we define a set of reasoning rules (

Table 2. Reasoning rules determining the visit preference strength.

Demand	Available Capacity
	CL	CM	CH
DL	PM	PH	PVH
DM	PL	PM	PH
DH	PVL	PL	PM

) to determine the preference to visit the next node depending on the levels of both the demand and the vehicle available capacity.

Figure 5 displays the procedure used to compute the preference index $p_i$ after serving the customer $i\in I$. This procedure is described as follows:

• Simulate the actual demand of each customer employing a fuzzy simulation approach. Based on the works by Teodorović and Pavković [77], Sun et al. [59] and Sun [78], we follow the steps described below:

  (a)

  Generate a random demand $d_i$ between a lower bound and an upper bound. Since the objective is preserving the variability conditions similar to the stochastic demands, the lower and upper bounds are given by the expressions $\frac{{\varphi }_i-\sqrt{3\lambda {\varphi }_i}}{q}$ and $\frac{{\varphi }_i+\sqrt{3\lambda {\varphi }_i}}{q}$, respectively.

  (b)

  Calculate the membership degree $\mu \left(d_i\right)$ of this demand. Notice that $\mu \left(d_i\right)\in [0,1]$.

  (c)

  Generate a random number $\rho \in [0,1]$.

  (d)

  Compare $\rho$ and $\mu \left(d_i\right)$. If $\rho \le \mu \left(d_i\right)$, then assume the actual demand of the customer i as $d_i$; otherwise, repeat steps (a)–(d) until this condition is fulfilled.

• Calculate the vehicle available capacity subtracting from q the sum of the simulated demand of the first m customers visited in the current route, including the customer i. Whenever the route fails and the vehicle must perform a trip to the depot for a replenishment, the counting of m starts again from 1.
• Estimate the fuzzy demand and the fuzzy available capacity according to the categories previously defined: low, medium or high.
• Determine the membership function of the preference index using the reasoning rules defined in Table 2.
• Calculate a crisp preference index using the center of gravity as defuzzification method. Additional methods can be found in Klir and Yuan [79], and Opricovic and Tzeng [80].

We define a known threshold $p^{*}$, such that $0\le p^{*}\le 1$. The computed preference index $p_i$ must be compared with $p^{*}$ in order to make a decision about the vehicle next destination. If $p_i\ge p^{*}$, the vehicle should visit the next customer directly; otherwise, we estimate that the vehicle available capacity cannot meet the next customer demand. In this case, both preventive ($c_{prev}$) and reactive ($c_{reac}$) costs are calculated (see Section 3). If $c_{prev}<c_{reac}$, the vehicle should perform a detour to the depot for a preventive replenishment; otherwise, it should visit the next customer directly and react to its real demand. The lower the threshold level, the greater the inclination to unload the vehicle as much as possible before making a replenishment trip to the depot. In this case, less preventive detours are performed. Hence, the number of times that a reactive round-trip must be carried out increases. Previous tests using modified Akca’s instances yielded lower costs when $p^{*}=0.45$.

The following parameters are used by our algorithm to run the experiments: (i) 350 iterations for map perturbations; (ii) 150 iterations for the biased-randomized savings heuristic; (iii) 150 iterations for splitting; (iv) a random value between $0.05$ and $0.80$ for ${\beta }_1$, the parameter of the geometric distribution associated with the biased-randomized selection during the allocation map process; (v) a random value between $0.07$ and $0.23$ for ${\beta }_2$, the parameter of the geometric distribution associated with the biased-randomized heuristic for routing; (vi) $n=100$ runs for the initial simulation stage; (vii) $N=5000$ runs for the intensive simulation stage; and (viii) 100 iterations to estimate the safety stock (SS), testing only discrete values between 0% and 10%. Our proposed algorithm was coded as a Java application. All experiments were executed on a standard Windows PC with a Core $i5$ processor and 6 GB RAM. A total of ten different random seeds were used for each instance.

6. Results and Discussion

Table 3. Comparative results between our flexible solutions under different uncertainty levels.

Instance	Best Deterministic Solution			Best Stochastic Solution [13]					Best Hybrid Solution					Best Fuzzy Solution
	OC	RC	TC	OC	RC	FC	TC	SS	OC	RC	FC	TC	SS	OC	RC	FC	TC	SS
Low variability
Cr30x5a-1	200.00	575.14	775.14	200.00	575.14	2.37	777.51	0%	200.00	575.14	3.31	778.45	2%	200.00	575.14	5.86	781.00	2%
Cr30x5a-2	200.00	607.28	807.28	200.00	607.28	0.04	807.32	3%	200.00	607.28	0.12	807.40	3%	200.00	607.28	0.12	807.40	3%
Cr30x5a-3	187.50	507.92	695.42	187.50	509.25	10.99	707.74	3%	187.50	509.25	17.48	714.22	3%	187.50	509.25	25.50	722.25	3%
Cr30x5b-1	225.00	623.22	848.22	225.00	623.22	9.37	857.59	0%	225.00	623.22	14.59	862.81	0%	225.00	623.22	22.85	871.07	1%
Cr30x5b-2	187.50	625.32	812.82	187.50	625.32	0.00	812.82	2%	187.50	625.32	0.00	812.82	2%	187.50	625.32	0.00	812.82	2%
Cr30x5b-3	187.50	684.58	872.08	187.50	684.58	2.25	874.33	1%	187.50	684.58	6.35	878.43	1%	187.50	684.58	9.50	881.58	1%
Cr40x5a-1	162.50	731.84	894.34	162.50	731.84	0.03	894.37	1%	162.50	731.84	0.07	894.41	1%	162.50	731.84	0.59	894.93	1%
Cr40x5a-2	225.00	637.26	862.26	225.00	639.02	0.10	864.12	0%	225.00	639.02	0.81	864.83	1%	225.00	642.02	0.03	867.05	3%
Cr40x5a-3	162.50	752.88	915.38	162.50	752.88	0.97	916.35	0%	162.50	752.88	3.26	918.64	0%	162.50	752.88	6.82	922.21	1%
Cr40x5b-1	162.50	852.04	1014.54	162.50	852.04	6.90	1021.45	1%	162.50	852.04	12.24	1026.78	1%	162.50	852.04	20.79	1035.33	1%
Cr40x5b-2	225.00	690.57	915.57	225.00	690.57	0.08	915.65	1%	225.00	690.57	0.62	916.18	1%	225.00	690.57	1.23	916.79	1%
Cr40x5b-3	175.00	764.33	939.33	175.00	772.87	0.07	947.93	2%	175.00	772.87	0.29	948.16	2%	175.00	772.87	0.35	948.22	2%
Average	191.67	671.03	862.70	191.67	672.00	2.76	866.43	1.17%	191.67	672.00	4.93	868.59	1.42%	191.67	672.25	7.80	871.72	1.75%
Medium variability
Cr30x5a-1	200.00	575.14	775.14	200.00	575.14	7.63	782.77	2%	200.00	575.14	9.67	784.81	2%	200.00	575.14	12.91	788.05	2%
Cr30x5a-2	200.00	607.28	807.28	200.00	607.28	0.46	807.74	3%	200.00	607.28	1.94	809.22	3%	200.00	607.28	1.43	808.71	3%
Cr30x5a-3	187.50	507.92	695.42	187.50	509.25	18.50	715.25	3%	187.50	509.25	24.10	720.85	3%	187.50	509.25	29.73	726.48	3%
Cr30x5b-1	225.00	623.22	848.22	225.00	623.22	14.63	862.85	0%	225.00	623.22	18.32	866.53	3%	225.00	623.22	24.23	872.45	3%
Cr30x5b-2	187.50	625.32	812.82	187.50	625.32	0.00	812.82	2%	187.50	625.32	0.00	812.82	2%	187.50	625.32	0.00	812.82	2%
Cr30x5b-3	187.50	684.58	872.08	187.50	684.58	10.21	882.28	0%	187.50	684.58	12.79	884.87	1%	187.50	684.58	12.88	884.96	1%
Cr40x5a-1	162.50	731.84	894.34	162.50	739.24	0.01	901.75	3%	162.50	739.24	0.01	901.75	3%	162.50	739.24	0.00	901.74	3%
Cr40x5a-2	225.00	637.26	862.26	225.00	643.52	3.07	871.59	1%	225.00	642.02	0.24	867.26	3%	225.00	642.02	0.57	867.59	3%
Cr40x5a-3	162.50	752.88	915.38	162.50	752.88	4.46	919.85	1%	162.50	752.88	8.57	923.95	1%	162.50	752.88	11.83	927.22	1%
Cr40x5b-1	162.50	852.04	1014.54	162.50	858.58	4.54	1025.62	2%	162.50	858.58	8.01	1029.09	2%	237.50	795.18	0.00	1032.68	4%
Cr40x5b-2	225.00	690.57	915.57	225.00	690.57	2.06	917.63	1%	225.00	690.57	3.77	919.33	0%	225.00	690.57	5.80	921.37	1%
Cr40x5b-3	175.00	764.33	939.33	175.00	772.87	1.42	949.29	2%	175.00	772.87	2.53	950.40	2%	175.00	772.87	2.96	950.82	2%
Average	191.67	671.03	862.70	191.67	673.54	5.58	870.79	1.67%	191.67	673.41	7.50	872.57	2.08%	197.92	668.13	8.53	874.57	2.33%
High variability
Cr30x5a-1	200.00	575.14	775.14	200.00	575.14	19.66	794.80	2%	200.00	575.14	19.82	794.96	0%	200.00	575.14	24.25	799.38	1%
Cr30x5a-2	200.00	607.28	807.28	200.00	607.74	0.02	807.76	5%	200.00	611.41	0.02	811.43	7%	200.00	607.74	0.04	807.78	5%
Cr30x5a-3	187.50	507.92	695.42	187.50	509.25	27.86	724.61	2%	187.50	509.25	29.95	726.70	4%	187.50	509.25	33.41	730.16	3%
Cr30x5b-1	225.00	623.22	848.22	225.00	623.22	19.99	868.21	10%	225.00	623.22	20.73	868.95	10%	225.00	623.22	24.86	873.08	10%
Cr30x5b-2	187.50	625.32	812.82	187.50	625.32	0.10	812.92	3%	187.50	625.32	0.20	813.02	5%	187.50	625.32	0.15	812.97	3%
Cr30x5b-3	187.50	684.58	872.08	187.50	684.58	24.93	897.00	1%	187.50	684.58	29.03	901.11	5%	187.50	684.58	34.01	906.09	5%
Cr40x5a-1	162.50	731.84	894.34	162.50	737.20	2.85	902.55	2%	162.50	735.84	7.83	906.17	1%	162.50	735.84	9.38	907.71	1%
Cr40x5a-2	225.00	637.26	862.26	225.00	642.02	1.79	868.82	3%	225.00	642.02	1.48	868.50	3%	225.00	642.02	2.25	869.27	3%
Cr40x5a-3	162.50	752.88	915.38	162.50	763.69	5.78	931.97	2%	162.50	763.69	7.76	933.96	2%	162.50	752.88	18.65	934.04	1%
Cr40x5b-1	162.50	852.04	1014.54	237.50	786.00	4.65	1028.14	3%	237.50	792.36	2.84	1032.70	4%	237.50	786.00	8.47	1031.97	3%
Cr40x5b-2	225.00	690.57	915.57	225.00	690.57	9.35	924.91	2%	225.00	690.57	12.59	928.15	2%	225.00	690.57	14.96	930.53	2%
Cr40x5b-3	175.00	764.33	939.33	175.00	780.62	4.14	959.76	3%	175.00	780.62	4.90	960.52	3%	175.00	780.62	5.86	961.48	3%
Average	191.67	671.03	862.70	197.92	668.78	10.09	876.79	3.17%	197.92	669.50	11.43	878.85	3.83%	197.92	667.77	14.69	880.37	3.33%

shows our obtained results for 12 Akca’s instances. Five main indicators are computed: depots opening costs (OC), which is formed by both fixed and variable costs; routing costs (RC); failure costs (FC), which is incurred whenever the vehicle must perform either a detour or a round-trip to the depot; total costs (TC); and the estimated safety stock (SS) level. Four types of solutions are compared. All of them are flexible, i.e., they consider facility sizing decisions. Firstly, our best deterministic solutions are shown, i.e, there is no uncertainty in the customers’ demand and its realization is exactly as expected. In this case, a safety stock is not necessary and there are no failure costs. Secondly, we show the best stochastic solutions reported by Tordecilla et al. [13], in which the exact customers’ demand is not known. Instead, all of them follow a log-normal distribution with known mean and standard deviation. Thirdly, our best hybrid fuzzy-stochastic solutions are displayed, in which half of the customers’ demand follows a log-normal distribution, and half of the customers’ demand is considered to be fuzzy. Finally, our best fuzzy solutions are shown, in which all customers’ demand is considered to be fuzzy, due to a high level of uncertainty. Additionally, results for three levels of variability ($\lambda$) are shown. Clearly, our best deterministic solutions are the same regardless of the variability level, given the total absence of uncertainty.

Results in Table 3 show a slight average increase in total costs when increasing the variability level for all types of solutions, except the best deterministic solution. This growth is caused mainly by the rise in failure costs, since a greater number of detours and round-trips is expected when the demand variability level is higher. Additionally, total costs also increase when the uncertainty level is higher regardless of the variability level, i.e., the deterministic solution is the cheapest one, and the fuzzy solution is the most costly. If we compare only the average deterministic cost of each set of solutions, formed by the sum of OC and RC, we obtain values with negligible differences. Hence, the contrasts in total costs are caused mainly by failure costs. For example, for the instance Cr30x5a-3 in the low variability scenario, 1.6% of total costs are failure costs in the best stochastic solution. However, in the best fuzzy solution this percentage rises to 3.5%. Most instances show this steady growth when increasing the uncertainty level, which confirms that fuzzy scenarios have a higher uncertainty level when compared with deterministic and stochastic scenarios. Finally, the average safety stock increases when both variability and uncertainty levels rise, since more protection against uncertainty is necessary in both cases.

Results corresponding to our best deterministic solution in Table 3 were yielded assuming that the realized demand is deterministic. Hence, an additional experiment has been performed, in which this solution (called henceforth OBD) is tested in a hybrid fuzzy-stochastic environment, using 0% of safety stock protection against uncertainty. Figure 6 compares this solution’s results with our best-found hybrid fuzzy-stochastic solution (OBF) in terms of failure costs. Results for 12 Akca’s instances are depicted for each demand variability scenario. Extreme points in dashed lines indicate the average cost for each set of data. As expected, average failure costs show an increasing trend when the variability grows, regardless of the type of solution. Conversely, Figure 6 shows that OBF outperforms OBD when tested under uncertainty conditions. This fact demonstrates the quality of our fuzzy simheuristic approach, especially in scenarios where the demand variability is high.

Table 4. Comparative results between our hybrid solutions when considering facility sizing decisions.

Instance	Best Nonflexible Hybrid Solution					Best Flexible Hybrid Solution					Gap TC
	OC	RC	FC	TC	SS	OC	RC	FC	TC	SS
Low variability
Cr30x5a-1	200.00	619.51	3.45	822.96	1%	200.00	575.14	3.31	778.45	2%	−5.41%
Cr30x5a-2	200.00	626.01	0.04	826.05	1%	200.00	607.28	0.12	807.40	3%	−2.26%
Cr30x5a-3	200.00	507.99	17.56	725.55	2%	187.50	509.25	17.48	714.22	3%	−1.56%
Cr30x5b-1	200.00	682.97	0.32	883.29	2%	225.00	623.22	14.59	862.81	0%	−2.32%
Cr30x5b-2	200.00	625.32	0.00	825.32	2%	187.50	625.32	0.00	812.82	2%	−1.51%
Cr30x5b-3	200.00	684.58	5.95	890.53	1%	187.50	684.58	6.35	878.43	1%	−1.36%
Cr40x5a-1	200.00	733.47	3.22	936.70	0%	162.50	731.84	0.07	894.41	1%	−4.51%
Cr40x5a-2	200.00	691.47	11.15	902.63	1%	225.00	639.02	0.81	864.83	1%	−4.19%
Cr40x5a-3	200.00	748.64	9.88	958.52	1%	162.50	752.88	3.26	918.64	0%	−4.16%
Cr40x5b-1	200.00	858.58	1.94	1060.53	2%	162.50	852.04	12.24	1026.78	1%	−3.18%
Cr40x5b-2	300.00	690.57	0.65	991.22	2%	225.00	690.57	0.62	916.18	1%	−7.57%
Cr40x5b-3	200.00	780.62	0.07	980.69	2%	175.00	772.87	0.29	948.16	2%	−3.32%
Average	208.33	687.48	4.52	900.33	1.42%	191.67	672.00	4.93	868.59	1.42%	−3.45%
Medium variability
Cr30x5a-1	200.00	619.51	9.17	828.68	0%	200.00	575.14	9.67	784.81	2%	−5.29%
Cr30x5a-2	200.00	626.01	0.60	826.61	2%	200.00	607.28	1.94	809.22	3%	−2.10%
Cr30x5a-3	200.00	507.99	24.30	732.29	2%	187.50	509.25	24.10	720.85	3%	−1.56%
Cr30x5b-1	200.00	681.50	14.31	895.80	1%	225.00	623.22	18.32	866.53	3%	−3.27%
Cr30x5b-2	200.00	625.32	0.01	825.33	2%	187.50	625.32	0.00	812.82	2%	−1.52%
Cr30x5b-3	200.00	684.58	15.60	900.18	1%	187.50	684.58	12.79	884.87	1%	−1.70%
Cr40x5a-1	200.00	733.47	7.69	941.17	1%	162.50	739.24	0.01	901.75	3%	−4.19%
Cr40x5a-2	200.00	700.80	12.59	913.39	3%	225.00	642.02	0.24	867.26	3%	−5.05%
Cr40x5a-3	200.00	748.64	20.15	968.79	0%	162.50	752.88	8.57	923.95	1%	−4.63%
Cr40x5b-1	200.00	863.91	2.32	1066.23	3%	162.50	858.58	8.01	1029.09	2%	−3.48%
Cr40x5b-2	300.00	690.57	4.18	994.75	1%	225.00	690.57	3.77	919.33	0%	−7.58%
Cr40x5b-3	200.00	780.62	0.94	981.56	3%	175.00	772.87	2.53	950.40	2%	−3.17%
Average	208.33	688.58	9.32	906.23	1.58%	191.67	673.41	7.50	872.57	2.08%	−3.63%
High variability
Cr30x5a-1	200.00	619.51	20.69	840.20	0%	200.00	575.14	19.82	794.96	0%	−5.38%
Cr30x5a-2	200.00	621.45	5.66	827.12	3%	200.00	611.41	0.02	811.43	7%	−1.90%
Cr30x5a-3	200.00	507.99	30.16	738.15	4%	187.50	509.25	29.95	726.70	4%	−1.55%
Cr30x5b-1	200.00	681.50	18.85	900.35	0%	225.00	623.22	20.73	868.95	10%	−3.49%
Cr30x5b-2	200.00	625.32	0.14	825.46	5%	187.50	625.32	0.20	813.02	5%	−1.51%
Cr30x5b-3	200.00	684.58	30.23	914.81	1%	187.50	684.58	29.03	901.11	5%	−1.50%
Cr40x5a-1	200.00	737.94	5.78	943.73	2%	162.50	735.84	7.83	906.17	1%	−3.98%
Cr40x5a-2	200.00	700.80	15.98	916.78	3%	225.00	642.02	1.48	868.50	3%	−5.27%
Cr40x5a-3	200.00	748.64	32.89	981.54	0%	162.50	763.69	7.76	933.96	2%	−4.85%
Cr40x5b-1	200.00	858.58	22.53	1081.11	2%	237.50	792.36	2.84	1032.70	4%	−4.48%
Cr40x5b-2	300.00	693.03	12.66	1005.69	0%	225.00	690.57	12.59	928.15	2%	−7.71%
Cr40x5b-3	200.00	772.87	13.22	986.09	2%	175.00	780.62	4.90	960.52	3%	−2.59%
Average	208.33	687.68	17.40	913.42	1.83%	197.92	669.50	11.43	878.85	3.83%	−3.68%

compares two types of hybrid fuzzy-stochastic solutions. Firstly, we show our best solution with a single facility size alternative given by the original Akca’s instances—i.e., the solution is not flexible since only one size is available to select. Secondly, we show our best flexible solution, which corresponds to our best hybrid solution in Table 3. When comparing the total costs of both types of solutions, the negative gap obtained for all instances and under all variability levels shows the advantages of considering facility sizing decisions. For example, we reach a maximum absolute gap of $7.71\%$ in total cost savings for a single instance. In average, both opening and routing costs decrease whenever alternative depot sizes are available. Nevertheless, each instance shows different results regarding OC and RC. The most evident case is that in which opening costs decrease. Clearly, this is a direct result of having smaller facility size alternatives. Without loss of generality, all examples below take as reference the high variability scenario. For example, the instance Cr30x5b-3 has a total demand of 1620. Both flexible and nonflexible approaches design the same routes and yield equal routing costs. Nevertheless, the nonflexible approach locates two depots of size 1000 each. Conversely, our flexible approach locates one depot of size 1000 and one depot of size 750. Hence, the nonflexible solution assigns an extra capacity that is not necessary under the problem’s current conditions.

Some instances show an opposite behavior, i.e., opening costs either increase or remain the same while routing costs decrease. For example, the nonflexible solution of the instance Cr30x5a-1 opens two depots of size 1000 each. Alternatively, the flexible solution opens one depot of size 1500 and one depot of size 500, i.e., the total capacity is equal and, given our defined costs structure, also the opening costs. However, this slight change drives a redesign of routes that decreases RC. An additional example is given by the instance Cr40x5a-2. Figure 7 depicts the best solution found by both the nonflexible approach (a) and our flexible approach (b). The solution in Figure 7a locates two depots of size 1750 each, and the solution in Figure 7b locates three depots of size 875 each. The latter case has a total capacity that is smaller than the former’s; however, opening costs are higher since the fixed cost is clearly greater when 3 facilities are open instead of 2. This new configuration decreases considerably routing costs (Table 4), which shows that considering facility sizing decisions not only reduces total costs by decreasing depots capacity but also by increasing it, since shorter routes can be designed.

Managerial Insights

From a managerial perspective, we have shown a general algorithm useful to solve a flexible-size LRP where a subset of customers provides enough information to model stochastically their demand, while the complementary subset provides scarce data. In this case, decision makers may estimate a fuzzy demand. Our algorithm is general because scenarios where the demand of all customers is deterministic, stochastic or fuzzy represent particular cases of our described problem. Hence, decision makers can employ our approach more extensively than other algorithms. We analyze these scenarios through some numerical results and assess how the level of uncertainty influences opening, routing and route-failure costs. Clearly, more precise data decrease total costs. Furthermore, we have calculated the cost of assuming a deterministic demand when the real scenario is fuzzy or stochastic. It has been shown that our hybrid approach yields less average costs, which leads to a more competitive supply chain. Additionally, we have also shown that important cost savings are generated whenever a set of facility size alternatives are analyzed by decision makers, instead of considering a single alternative—as in most LRP studies. Finally, our algorithm is able to generate detailed information about the location-allocation-routing decisions that should be made.

7. Conclusions

This work presented a location routing problem where the facility size is an additional variable, instead of a known parameter as the traditional LRP assumes. Moreover, we consider a hybrid fuzzy-stochastic setting in which some customers’ demands are fuzzy and others are stochastic. Hence, a fuzzy simheuristic approach is proposed to solve this problem cost- and time- efficiently. Initially, our algorithm selects the best size for each open facility from a set of provided alternatives. We perform an iterative procedure in which a set of location-allocation-routing configurations are assessed in terms of opening and routing costs. Then a top list of complete LRP solutions is iteratively perturbed and simulated. The perturbation stage is performed by employing an iterated local search metaheuristic. The simulation stage is carried out by running a classic Monte Carlo simulation for the stochastic demands and a fuzzy simulation for the fuzzy demands. Failure costs are introduced as an additional performance indicator. Finally, a set of elite solutions is assessed through a refinement procedure where a larger number of simulation runs is executed.

Our fuzzy simheuristic approach has been proved to be flexible enough not only to combine efficiently stochastic and fuzzy demands in a single execution but also to address less general scenarios in which demands of all customers are either deterministic or fuzzy. Our approach has also been proved to be a cost-efficient algorithm when considering uncertainty scenarios. It decreases route failure costs when compared with the best deterministic solution tested in a hybrid fuzzy-stochastic environment. The use of a safety stock policy as a protection against uncertainty has also contributed to this decrease. In order to design a time-efficient algorithm, our current approach employs stochastic and fuzzy simulation only to assess the designed routes. Hence, our algorithm results can be enhanced by introducing fuzzy-stochastic aspects from the construction stage. However, this approach might also increase computational times.

To the best of our knowledge, this is the first time that a hybrid fuzzy-stochastic LRP with facility sizing decisions is addressed. Medium-sized benchmark instances considering three demand variability levels were used. Obtained results show that introducing such flexibility decreases total costs in two mutually nonexclusive ways: firstly, yielding savings in opening costs by locating facilities of smaller size; and secondly, yielding savings in routing costs by locating facilities of higher size, which drives a routes redesign that reduces the total traveled distance. We also have demonstrated that these savings are always incurred regardless of the demand variability level.

Multiple challenges remain open for future research. Since we are considering that only routes fail when demands are higher than expected, future work can include the simulation of facility failures, which would prompt a revision of location-allocation decisions. In addition, failure costs are currently measured only by considering the distances traveled to perform round-trips and detours. Still, real-life customers might not allow a delivery delay, e.g., because a time windows constraint must be met. This delay may drive lost sales or a goodwill reduction. Hence, this type of costs can be included in the computation of failure costs. Finally, large-sized instances can be used to assess the influence of the number of nodes in our approach performance.