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

: The location routing problem integrates both a facility location and a vehicle routing problem. Each of these problems are NP-hard in nature, which justiﬁes 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.


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 flexiblesize 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 fuzzystochastic 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.

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.

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 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 fuzzystochastic 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.

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.

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   Formally speaking, the LRP can be defined on a complete, weighted, and undirected graph G(V, E, C), in which V = J ∪ 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 ∈ 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 ∪ 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 ∈ L D i = Uncertain demand of customer i ∈ I f j = Fixed opening cost of depot j ∈ J o jl = Variable opening cost of depot j ∈ J with size of type l ∈ L c e = Cost of traversing arc e ∈ E q = Capacity of each vehicle %SS = Safety stock percentage Decision variables y jl = Binary variable that indicates whether the depot j ∈ J is open with size l ∈ L or not.
x ij = Binary variable that indicates whether customer i ∈ I is assigned to the depot j ∈ J or not w ek = Binary variable that indicates whether arc e ∈ E is used in the route performed by vehicle k ∈ 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).
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:

1.
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.

2.
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 ∅ = S ⊂ V be a subset of nodes, δ + (S) the set of edges leaving S, δ − (S) the set of edges entering S, and A(S) 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: subject to: 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.

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 (nbDepots 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 nbDepots 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 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, α, β, λ, Inc, T 0 , K min , K max , I 0 , t max ) 1 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.

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:

1.
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 ∈ J, and L is the set of available sizes, our approach' alternative sizes are s jl ∈ {(1 − 2r)s j , (1 − r)s j , s j , (1 + r)s j , (1 + 2r)s j }, where l ∈ 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.

2.
Traditional LRP instances consider a fixed cost ( f j ) incurred whenever a depot j ∈ 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 = (s jl − s j ) 2s j ∑ j f j |J| . This formula preserves o jl in the same order as f j for each depot j ∈ 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.

3.
An uncertain demand D i for each customer i ∈ I is considered. The demand of half of the customers is assumed to follow a log-normal probability distribution. If φ i is the deterministic demand in the Akca's set, then E[D i ] = φ i . In addition, three different values of variance are considered: low, medium and high, i.e., for λ ∈ {0.05, 0.10, 0.20}, 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 ≤ D i ≤ 1. The membership function of these fuzzy sets are displayed in Figure 2.

A 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 ∈ I. C i is expressed as a proportion of q, i.e., 0 ≤ C i ≤ 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 ≤ p i ≤ 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) 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 ∈ I. This procedure is described as follows:

1.
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:

2.
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.

3.
Estimate the fuzzy demand and the fuzzy available capacity according to the categories previously defined: low, medium or high.

4.
Determine the membership function of the preference index using the reasoning rules defined in Table 2.

5.
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 ≤ p * ≤ 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 ≥ 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 β 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 β 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. Table 3 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 (λ) are shown. Clearly, our best deterministic solutions are the same regardless of the variability level, given the total absence of uncertainty.   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 in
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 fuzzystochastic 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 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.

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. Funding: This work has been partially supported by the Spanish Ministry of Science (PID2019-111100RB-C21/AEI/10.13039/501100011033). In addition, it has received the support of the Doctoral School at the Universitat Oberta de Catalunya (Spain) and the Universidad de La Sabana (INGPhD-12-2020).