## 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) metheuristic with Monte Carlo simulation and fuzzy inference systems (FIS) to deal with stochastic and fuzzy variables, respectively. As discussed in [PERSON] 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.