Fuzzy Simheuristics for Optimizing Transportation Systems: Dealing with Stochastic and Fuzzy Uncertainty

Abstract

In the context of logistics and transportation, this paper discusses how simheuristics can be extended by adding a fuzzy layer that allows us to deal with complex optimization problems with both stochastic and fuzzy uncertainty. This hybrid approach combines simulation, metaheuristics, and fuzzy logic to generate near-optimal solutions to large scale NP-hard problems that typically arise in many transportation activities, including the vehicle routing problem, the arc routing problem, or the team orienteering problem. The methodology allows us to model different components–such as travel times, service times, or customers’ demands–as deterministic, stochastic, or fuzzy. A series of computational experiments contribute to validate our hybrid approach, which can also be extended to other optimization problems in areas such as manufacturing and production, smart cities, telecommunication networks, etc.

1. Introduction

Managers tend to rely on analytical methods that allow them to make informed decisions. This explains why optimization models play a key role in many industries and business, including the logistics and transportation sector. Whenever accurate information on the inputs and constraints of the optimization problem is available, the resulting deterministic models can be solved by using well-known methods, either of exact or approximate nature.

Many optimization problems in real-life transportation involve taking into account a large number of variables and rich constraints, which often makes them to be NP-hard [1]. When this is the case, the computational complexity makes it difficult to obtain optimal solutions in a short computational time. At this point, heuristic approaches can provide near-optimal solutions that, in turn, cover all the requirements of the problem [2]. When dealing with challenging optimization problems, there is a tendency to divide them into sub-problems, which simplifies the difficulty but might also lead to sub-optimal solutions [3,4]. Given the increase in computational power experienced during the last decade, and also the development of advanced metaheuristic algorithms, it is possible nowadays to solve rich and large-scale problems that were intractable in the past [5].

In the scientific literature on combinatorial optimization problems, it is often assumed that the input values are constant and known. However, in a real-world scenario this is rarely the case, since uncertainty is often present and affects these inputs. In the context of transportation and logistics, some examples of these inputs are: travel times, customer demands, service times, battery durability, etc. Whenever these inputs can be modeled by random variables, simheuristic algorithms—which combine heuristics with simulation—become a useful tool to address the associated optimization problem [6]. It should be noticed that simheuristics are designed to handle situations where uncertainty can be modeled by random variables, each of which follows a well-known probability distribution. When dealing with non-probabilistic uncertainty, fuzzy techniques might be a good choice. Therefore, fuzzy techniques can be particularly interesting for modeling uncertainty whenever it cannot be represented by random variables, for example: if not enough data are available, if the data cannot be fitted to a probability distribution, or if qualitative expert opinions must also be considered. Tordecilla et al. [7] illustrate with an example how to combine two types of uncertainty conditions using a fuzzy simheuristics, which hybridizes a metaheuristic with simulation and fuzzy logic. In their example, these authors assume that only some customer demands can be modeled by random variables, while others follow a fuzzy pattern.

A fuzzy system is based on fuzzy logic. Inputs enter the system, which computes fuzzy outputs on the basis of a set of rules established by a human expert [7]. In order to obtain solutions that mix information from different sources, the output of the fuzzy system includes different degrees of membership for different groups. This means that a fuzzy system can handlee decisions in a non-binary logic scenario, since the outputs have a partial degree of being ‘true’ or ‘false’. Therefore, the main contribution of this paper is to provide both conceptual and practical insights on how fuzzy simheuristics can be applied in the optimization of different transportation systems, which include the well-known vehicle routing problem (VRP) under uncertainty conditions, as well as the team orienteering problem (TOP) under uncertainty conditions. A comprehensive introduction to both problems can be found in Toth and Vigo [8] and Chao et al. [9], respectively. Therefore, we address and discuss the novel concept of fuzzy simheuristics, which has hardly been addressed in the literature. Accordingly, this new class of solution methodology is designed to solve the aforementioned transportation problems, whose performance and prospects have been duly analyzed and presented.

The remaining sections of this paper are organized as follows: Section 2 provides a description of the optimization problems discussed in this paper, the VRP and the TOP. Section 3 reviews related work on simheuristics and fuzzy sets in solving the aforementioned problems. The fuzzy simheuristic methodology is explained in Section 4. Section 5 describes how the proposed fuzzy simheuristic has been implemented, as well as the process of converting deterministic benchmarks into stochastic-fuzzy ones. A series of numerical experiments are included in Section 6. Finally, Section 7 summarizes the conclusions and main results of this work.

2. Popular Optimization Problems in Transportation

This section provides an overview of the two transportation problems considered in this paper, the VRP and the TOP.

2.1. The Vehicle Routing Problem

The VRP is a well-known combinatorial optimization problem with a vast number of applications in the transportation sector [10]. Solving the VRP aims to design cargo vehicle routes with minimum transportation costs to distribute goods between depots and a set of consumers. Since the capacity of the cargo vehicles is usually taken into account, the VRP is often referred to as capacitated VRP.

In its basic version, the distribution network of the VRP conists of a single depot and a set of customers, geographically distributed around a coverage area. A set of cargo vehicles, initially available at a central depot, visits customers to meet their demands. Once all customers assigned to a vehicle have been served, the vehicle returns to the central depot. The typical goal is to minimize the cost of distribution, serving all customers and without exceeding the loading capacity of the vehicles (which may or may not be homogeneous). This distribution network can be defined as a directed graph $G=(N,E)$, where: (i) $N=\{0,1,\dots ,|C\left|\right\}$ is the set of vertices, with node 0 being the central depot and C being the set of customers; and (ii) $E=\left\{\right(i,j)\mid i,j\in N,i<j\}$ is the set of edges connecting pairs of nodes. Each customer $i\in C$ requires a demand $d_i>0$, which affects the of the vehicle. The objective, in solving this problem, is to minimize the total cost of serving all customers, subject to: (i) each route starts and ends at the central depot; (ii) each customer is visited only once and by exactly one vehicle; and (iii) the total demand required by the costumers on a route does not exceed the vehicle capacity. Apart from this basic version, multiple extensions of the problem can be found in the literature, to name a few: heterogeneous fleet of vehicles [11,12], time-windows [13,14], multiple depots [15,16], multiple delivery levels [17,18] simultaneous pick-up and deliveries [19,20], or combination of the former [21,22,23]. Many real-life versions of this problem combine some of the aforementioned constraints and others, making them difficult to solve [24].

Figure 1 depicts the network topology of a basic VRP, in which a fleet of three vehicles is used to serve, from a central depot, 12 customers. The vehicles depart from the depot with loaded demand, which is delivered at each customer location. In this case, the routes are constrained by the capacity of cargo vehicles.

2.2. The Team Orienteering Problem

One of the main differences between the TOP and the VRP is that in the former it is not mandatory to visit all customers. In other words, some nodes can be omitted during the generation of the routing plan. This is due to restrictions on the fleet size and on the maximum length that can be covered by any route. In a typical TOP, rewards are collected the first time a node is visited and therefore the objective is to maximize the total reward collected by a fixed fleet of vehicles. Vehicles depart from an origin node and have to reach a destination node. Cargo constraints are not usually considered in the basic version of the TOP. However, as in the case of the VRP, many variants can be found in the literature, e.g.,: capacity constraints [25,26,27], maximum driving range [28], stochastic travel times [29], etc.

As in the VRP case, the TOP network is composed of an origin depot, a destination depot, and a set of customers. It can be defined as a directed graph $G=(N,E)$, where (i) $N=\{0,1,\dots ,|C|,n\}$ is the set of nodes, including the origin depot (node 0) and the destination depot (node n); and (ii) $E=\left\{\left\{\right(i,j)\mid i,j\in C,i<j\right\}$ represents the set of edges connecting the nodes. Each customer $i\in C$ offers a reward $u_i>0$ the first time it is visited. A fleet of vehicles is available at the origin depot. Each vehicle visits a selected subset of nodes in C to collect the associated rewards, and moves to the destination depot. This process is presented in Figure 2. Note that, unlike the network topology depicted in Figure 1, this time a subset of nodes is not visited. Typically, these are customers located far away from each other or offering low rewards.

3. Related Work

Real-life transportation problems deal with uncertain parameters, such as customer demands, travel and service times, resource availability, etc. Modeling and solving problems that address these types of parameters are difficult challenges, given the high complexity and large scale that even the deterministic version of these problems tend to have in most practical applications. To cope with this uncertainty, Oliva et al. [30] introduces the concept of “fuzzy simheuristics”, which is a general approach that considers both stochastic and fuzzy uncertainty. However, the vast majority of the literature on optimization of transportation systems does not consider the combination of probabilistic and non-probabilistic uncertainty in the same problem.

3.1. Simheuristics in Transportation Problems

Combining metaheuristics with simulation has proven to be a successful approach when dealing with combinatorial optimization problems (COPs) involving probabilistic uncertainty.

Simheuristics can be considered a very efficient approach to deal with stochastic COPs. Such efficiency is measured both in terms of computing time and solution quality, i.e.,: (i) simheuristics consume relatively low computing times given the inclusion of fast metaheuristics in the solution search procedure. Furthermore, promising initial solutions are evaluated by an initial simulation using a small number of runs (i.e., more intensive simulations are reserved only for a small group of elite solutions); and (ii) the inclusion of metaheuristics also has an impact on enhancing the solution quality. Moreover, since considering stochastic inputs implies that the outputs are also stochastic, simheuristics not only assess quality of the solution in terms of traditional indicators, such as costs or profits, but also in terms of risk and reliability values [31].

Simheuristics have been successfully employed to solve problems related to different application fields, such as flow shop scheduling [32,33], job shop scheduling [34], waste collection [35,36,37], hazardous waste management [38], facility location [39], military applications [40], healthcare [41], finance [42], telecommunication networks [43], or disaster management [44]. Nevertheless, simheuristics have been mainly applied to the optimization of transportation systems. Different variants can be found in the literature. For instance, Latorre-Biel et al. [45] combine simheuristics with machine learning and Petri nets to solve a single-depot VRP with stochastic and correlated demands. The proposed algorithm is capable of forecasting both customer demands and their correlations. Stochastic demands have also been considered in the VRP with multiple depots [15]. Travel times have also been considered stochastic in the literature about VRPs. Different types of problems address this parameter, e.g., VRP in the context of the so-called omnichannel retailing mode with pick up and delivery [46], two-dimensional VRP [47], or routing of electric vehicles [48].

A natural and realistic extension of the VRP is achieved by including inventory management in transportation decisions. For instance, Gruler et al. [49] uses a simheuristic to solve a single-period inventory routing problem with stochastic demands. Stochastic demands have also been considered in the context of agri-food supply chains. In this case, the proposed approaches are tested by addressing a real-world case [50] or by using benchmark instances [51]. The latter work also includes perishable products. Finally, simheuristics have been used less frequently in other transportation and routing problems, such as the arc routing problem [52], the location routing problem [53], or the team orienteering problem [29].

3.2. Fuzzy Sets in the VRP and the TOP

Fuzzy techniques are typically used to deal with non-probabilistic uncertainty. This approach is especially useful when insufficient historical data are available as to model uncertainty with probability distributions. Early work shows the potential of using fuzzy sets in the VRP by considering travel times between nodes [54] or customer demands [55] as fuzzy. On the one hand, Teodorovic and Kikuchi [54] propose the introduction of fuzzy arithmetic operations in the savings computation when solving the VRP. On the other hand, Teodorović and Pavković [55] present a procedure to decide whether or not to visit the next customer on a route. To make this decision, these authors propose a preference index, which indicates the strength of the inclination to make this visit. The preference index is computed by comparing the fuzzy demand of the next customer and the fuzzy current load of the vehicle.

Uncertainty in customer demand is a recurrent topic in the VRP literature using fuzzy techniques. For instance, Erbao and Mingyong [56] propose a hybrid differential evolution algorithm to solve a VRP with uncertain demands, which is additionally formulated as a fuzzy chance constrained program. Different preference index thresholds are tested with the objective of minimizing the total distance traveled. A similar approach is employed by Cao and Lai [57] to solve an open vehicle routing problem with fuzzy demands. This parameter is also considered by Shi et al. [58], who address a home healthcare open VRP with time windows. A fuzzy chance constraint model is proposed, as well as a hybrid genetic algorithm. A set of benchmark instances is used to test their approach. Fuzzy customer demands are considered as well by Kuo et al. [59] and Werners and Drawe [60].

Time-windows constraints are also quite frequently considered uncertain in the fuzzy VRP literature. The idea is that the information about the earliest and latest times at which customers must be visited is imprecise or vague. For instance, Ghannadpour et al. [61] address a realistic multi-objective dynamic VRP with time windows. In this case, the time windows are related to the level of customer satisfaction. This satisfaction is sought to be maximized. In turn, the aim is to minimize the number of vehicles used, the total distance traveled, and the waiting time of the vehicles. A solving approach based on a genetic algorithm is proposed. The relation between the level of customer satisfaction and the fuzzy time windows is also examined by Tang et al. [62]. It is proposed a multi-objective model that seeks to both minimize the distance traveled and maximize the level of the customer service. Fuzzy time windows are also considered by Xu et al. [63], López-Castro and Montoya-Torres [64], and Brito et al. [65]. The latter authors also consider the vehicle capacity as a fuzzy parameter. Finally, fuzzy sets are additionally used to model parameters such as service times [66,67] and travel times [68].

Fuzzy techniques have hardly been used in the TOP. The TOP is similar to the VRP, but in the former a fixed fleet of vehicles needs to collect rewards by visiting customers, and since there is a maximum time or distance that each vehicle can cover, it is often the case that not all customers can be visited [9]. Hence, the main objective of the TOP is to maximize the collected reward without exceeding the route length threshold. The orienteering problem (OP) refers to the single-vehicle (and less challenging) version of the TOP. Verma and Shukla [69] and Ni et al. [70] consider OPs in which both the collected rewards and the travel times are fuzzy. The former authors propose a parallel algorithm as a solving approach, whereas the latter employ a genetic algorithm. Regarding the TOP, Brito et al. [71] propose a greedy randomized adaptive search procedure (GRASP) to solve this problem considering fuzzy rewards and fuzzy travel times. A fuzzy linear program is formulated to model the addressed problem, with the objective of maximizing the total collected reward. Oliva et al. [30] introduce the concept of “fuzzy simheuristics” to deal with the general case where both stochastic and fuzzy uncertainty is present, e.g., when the parameter(s) related to a subset of customers are stochastic, whereas the parameter(s) related to another subset of customers are fuzzy. Hence, considering all parameters of the problem as stochastic, fuzzy or deterministic are particular cases. These authors consider that the customer rewards are uncertain. Finally, fuzzy simheuristics are also employed by Tordecilla et al. [7] to solve a location routing problem in which the size of the depots is an additional variable to consider. Different variability levels are taken into account in order to test the behavior of the proposed solution approach when customer demands are uncertain. To the best of our knowledge, the works by Oliva et al. [30] and Tordecilla et al. [7] are the only papers dealing simultaneously with stochastic and fuzzy scenarios in NP-hard transportation or location problems. Hence, significant challenges open up for future research in this field.

4. A General Solving Methodology

Both the VRP and TOP are NP-hard problems [72,73]. Therefore, due to the combinatorial nature of these problems, the use of exact solving approaches is often limited by the size of the problem instances. When dealing which real-life instances, which are typically large-scale instances, the use of heuristics and metaheuristics has proven to be a good alternative [74]. Although metaheuristics are capable of finding near-optimal (or even optimal) solutions to many different combinatorial optimization problems in reasonable computing times, these approaches have been mainly designed to solve deterministic versions of these problems. Consequently, metaheuristics are not able to cope adequately with stochastic components, being their application constrained against uncertain scenarios as the ones proposed in this paper. In order to tackle uncertainty, metaheuristics have been combined with simulation methods in recent years. The resulting simheuristic approaches, apart from finding cost-effective solutions for the deterministic problem—through the optimization component—are also able to provide efficient solutions for the stochastic scenario [6].

In many real-life situations, large-scale and complex optimization problems assume different degrees of uncertainty, not only of a stochastic but also of a fuzzy (non-stochastic) nature. The latter might occur, for instance, when the volume of observations is low or the available data are of insufficient quality. In this case, the accurate modeling of the uncertainty sources simply does not follow the natural pattern of modeling them only as stochastic variables following a probability distribution. Instead, fuzzy inference systems (FIS) are also considered to achieve this goal.

In this paper, we propose a fuzzy simheuristic approach to solve both the VRP and the TOP under general uncertainty scenarios (i.e., those including both probabilistic as well as non-probabilistic uncertainty). This hybrid solution approach combines a multi-start (MS) metaheuristic with MCS and FIS to deal with stochastic and fuzzy variables, respectively. Specifically, this solution method is composed of three stages. The first stage refers to the construction of an initial feasible solution through a savings-based constructive heuristic, which is designed considering the characteristics of each problem. The second stage consists of enriching this heuristic with biased-randomized (BR) decisions [75], which are then incorporated into a MS framework–in order to generate multiple solutions. This stage, in addition to exploring different regions of the solution space, conducts a short number of simulations on a set of promising solutions in order to evaluate their efficiency under stochastic conditions. Finally, the third stage is a refinement one, in which a larger number of simulation runs are applied to a set of elite solutions. This procedure allows to obtain a more accurate estimation of the different solution properties. Since the number of solutions generated during the search can be large and the simulation process is time consuming, we have limited the number of short simulation runs to 100. Regarding the number of simulations in the refinement step, it has been set to 1000.

Figure 3 depicts a high-level description of the proposed methodology. As explained, this process starts from solving the deterministic problem, whose corresponding solution is submitted to a short simulation procedure, i.e., the exploratory stage. Consequently, new solutionsare generated for both the stochastic and the fuzzy environment. These steps are repeated until a stopping criterion is met. Finally, the best-found solutions (or a set of elite solutions) are submitted to a large number of simulation runs—the intensive stage—in order to obtain a more accurate summary of output variables, such as total cost/reward and risk/reliability values.

In order to facilitate reproducibility, the low-level details of each of the stages of the described methodology are provided below:

• The constructive heuristics for solving the VRP and TOP are based on the savings concept [76]. Despite being structurally similar for both problems, their particularities are introduced to adquately cope with each respective case, as follows:

  • Firstly, a dummy solution is constructed. This hypothetical solution is composed of a set of routes, each of them being designed to serve one customer. The vehicle departs from the origin depot, visits the customer and continues the trip towards the destination depot. In the case of the TOP, this stage takes into account the maximum tour length when designing these dummy routes. That is, those dummy routes whose total travel time is greater than this limit are automatically discarded. Similarly, in the case of the VRP, this stage takes into account the maximum loading capacity of each vehicle (i.e., if the demand of a customer is higher than this capacity, this customer is discared).

  • Secondly, a savings list ($\mathit{SL}$) is created, which includes all the edges connecting two different locations. For each edge $(i,j)\in \mathit{SL}$, a savings value is computed according to Equation (1), for the VRP, and (2), for the TOP. In both cases, $t_{\mathit{ij}}$ represents the time- or distance-based cost associated with traveling from node i to node j, 0 is the origin depot. In the case of the TOP, n represents the destination depot, while $u_i$ and $u_j$ represent the rewards obtained when customers i and j are visited for the first time. In the case of the TOP, considering a linear combination of both travel time and reward allows us to define an ‘enriched savings’ concept that reflects not only the desire of maximizing the total collected reward, but also takes into account how far a customer is from the rest of nodes on the emerging route [29]. Once computed, the $\mathit{SL}$ list is sorted in descending order of savings value, which implies that edges with the highest savings are placed at the top of the list.

    (1)$s_{\mathit{ij}}=t_{0i}+t_{j0}-t_{\mathit{ij}}$

    (2)$s_{\mathit{ij}}=\alpha (t_{\mathit{in}}+t_{0j}-t_{\mathit{ij}})+(1-\alpha )(u_i+u_j)$

  • The sorted $\mathit{SL}$ reflects the most promising movements to reduce the corresponding costs. In this way, the savings edge at the top of the $\mathit{SL}$ is selected to perform the merging of the associated routes. This procedure is performed only if a feasible combined route is generated. For the VRP, two routes can only be merged if the vehicle capacity is not exceeded. Alternatively, for the TOP, two routes can only be merged if the total travel distance does not exceed the maximum tour length. Once the selected savings edge is checked, it is deleted from the $\mathit{SL}$. Then, the new edge at the top is selected to continue this procedure, which is repeated until the $\mathit{SL}$ is empty. At the end of this process, a feasible solution is generated.

• The described heuristics are deterministic, which implies that the same decisions are made whenever they start from the same configuration. To change this behavior, these deterministic heuristics are transformed into a probabilistic algorithms by ‘smoothing’ the selection of candidates from the $\mathit{SL}$ using a probability distribution. This concept is called biased-randomization (BR), and is described in Dominguez et al. [77]. In our case, the geometric probability distribution was adopted, as suggested in Ferone et al. [78]. Introducing BR decisions in our heuristics requires dealing with additional parameters, such as the $\beta \in (0,1)$, which defines the geometric distribution. The value of this parameter was set after a quick tuning process over a random sample, establishing a good performance for both algorithms whenever $\beta$ falls in the interval $(0.3,0.4)$. Algorithm 1 describes the heuristic operational structure. Note that the difference between the two algorithms, designed to solve their respective problems, consists in how the $\mathit{SL}$ is constructed (line 2). Finally, the resulting BR algorithms are embedded into a multi-start (MS) framework in order to generate many alternative solutions. Then, the best-found solution is updated and returned at the end of this procedure.

  Algorithm 1: Biased-Randomized Algorithm.

  [Figure omitted. See PDF]

• The last stage refers to the incorporation of both simulation and fuzzy components into the BR-MS framework, so that promising solutions are processed to estimate their expected costs (Algorithm 2). For the VRP and TOP, the uncertain variables represent the customer demands and the travel times, respectively.

  • For stochastic variables, a new value is assigned to each random element based on its probability distribution. For stochastic variables, the MCS is used to estimate them.

  • For fuzzy variables, the new value of each element is based on its fuzzy function. Accordingly, fuzzy variables are estimated through the FIS. This procedure is explained more thoroughly in Section 5.

  Once the deterministic version of each problem is solved, their respective solutions are submitted to a exploratory stage, in which only a low number of simulation ($q_{\mathit{short}}$) runs are performed to avoid jeopardizing the time of the metaheuristic component [79]. These ‘short’ simulation runs are applied only to solutions that meet an acceptance criterion (line 8). Altering these stochastic and fuzzy values involves a re-evaluation of both the objective function and the constraints, so that the expected cost/reward of each promising solution can be computed. These short simulation runs allow multiple elite solutions to be found (line 11). In this way, once the BR-MS main loop is completed, a larger number of simulation ($q_{\mathit{long}}$) runs are executed for each elite solution (line 17). Consequently, the algorithm is able to obtain more accurate values of the output variables. Finally, a reduced set of best-found solutions is returned. From this set, managers can assess not only the expected costs/rewards but also the risk or reliability values associated with each solution, as described in Chica et al. [6].

Algorithm 2: Fuzzy Simheuristic.

[Figure omitted. See PDF]

5. Computational Experiments

The proposed fuzzy simheuristic has been implemented using Python $3.8$ and tested on a common PC with a multi-core processor Intel $i7$ and using 8 GB of RAM. The algorithm was executed five times with different seeds for a maximum time of 100 s for each instance. To the best of our knowledge, there are no instances in the literature for the stochastic-and-fuzzy problems described above. Accordingly, we have extended the well-known deterministic benchmarks proposed by Chao et al. [9] and Augerat et al. [80] for the TOP and the VRP problems, respectively. The following subsections describe in detail the process used to transform these deterministic benchmarks into stochastic-fuzzy ones.

5.1. A Fuzzy-Stochastic Approach for the VRP

In order to check the performance of our algorithm, we compare it with some benchmark instances that can be found in the literature. From Augerat et al. [80], we have chosen 29 of the classical instances that can be suitable for our study. The nomenclature of the instances is as follows: ‘$L-nXX-kY$’ where $L\in \{A,B,E\}$ is the set identification, $XX$ denotes the number of customers and Y establishes the number of vehicles. For carrying out the experiments, we assumed that the demand $d_i$ of each customer i is uncertain and, therefore, we have modeled it either as a stochastic or as a fuzzy variable.

Regarding the stochastic scenario, the instances were extended by considering that the stochastic demand $D_i$ follows a log-normal probability distribution. The parameters of this distribution were adjusted according to the mean $E\left[D_i\right]=d_i\forall i\in N$, where $d_i$ is the deterministic demand, and the variance $Var\left[D_i\right]=c·d_i$. The parameter c is a design parameter that allows us to set up the level of uncertainty. It is expected that, as c converges to zero, the results of the stochastic version will converge to those obtained in the deterministic scenario. In our experiments, we have utilized the value $c=0.25$, which introduces a medium level of uncertainty. Concerning to the fuzzy scenario, we consider the demand $D_i$, for each customer i, as a fuzzy variable. This demand can be estimated as low, medium, or high ($L,M,H$). Likewise, we assume that the vehicle remaining capacity, $RC$, is an input variable of the fuzzy system. Besides, each of the aforementioned demand levels is defined by a triangular fuzzy number $D_i=(d_{1i},d_{2i},d_{3i})$. Figure 4 shows the membership functions of these fuzzy sets. Similarly, the remaining vehicle capacity $RC$ is represented by a triangular fuzzy number $RC=(r{c}_1,r{c}_2,r{c}_3)$, which takes the values low (L), medium (M) or high (H) capacity. Figure 5 displays the membership function of the capacity fuzzy sets. Note that both the demands and the remaining capacities are expressed as a percentage of the total vehicle capacity, i.e., $0\le D_i\le 1$ and $0\le RC\le 1$.

For each node i, we define a preference index, $p_i$, as the output of the fuzzy system, such that $0\le p_i\le 1$. When this index takes the maximum value ($p_i=1$) then the next node of a route will be visited for sure as the remaining capacity $RC$ of the vehicle can meet the demand $D_{i+1}$. Moreover, if $p_i=0$, then we are sure that $D_{i+1}>RC$ and, consequently, the vehicle needs a replenishment at the depot. The preference index is classified into very low (VL), low (L), medium (M), high (H) and very high (VH) levels. The membership function related to each of these categories can be seen in Figure 6. The reasoning rules that determine the preference to travel to the next node–depending on the levels of the demand and the remaining capacity–are featured in Table 1. After performing a set of fine-tuning experiments, we established the threshold value to visit the next node to $p=0.25$. This means that whenever the calculated $p_i$ is greater than $0.25$, the next node will be visited; otherwise, the vehicle will return to the depot for a replenishment. The calculation of a specific value for $p_i$ requires converting the input variables into a crisp value. Hence, the estimated crisp values of the demand and the remaining capacity, the membership functions and the reasoning rules are employed in a fuzzification-defuzzification process to obtain the preference index. In our case, the defuzzification method applied was the well-known center-of-gravity method to obtain the output crisp value.

5.2. A Fuzzy-Stochastic Approach for the TOP

The deterministic benchmark used contains a total of 320 instances that are distributed in 7 subsets. The instances are identified following the nomenclature ‘$pa.b.c$’, where a represents the subset, b defines the number of available vehicles, and c identifies the specific instance under study. For experimentation purposes, we have considered that the uncertainty is located in the travel times between two pairs of nodes. To extend the instances to be employed in the stochastic scenarios, we have assumed that travel times, $T_{ij}$, follow a log-normal probability distribution. In setting up the stochastic instances, we assume that $E\left[T_{ij}\right]=t_{ij},$$\forall (i,j)\in N$, where $t_{ij}$ is the travel time for the corresponding deterministic instance. We set the variability in the travel times with reference to the deterministic travel time such that $Var\left[T_{ij}\right]=c·t_{ij}$, and $c\ge 0$. As in the VRP problem, we have utilized the value $c=0.25$ to induce a medium level of uncertainty in the travel times. With the aim of extending the instances to be used also in fuzzy scenarios, we have considered the travel times for each pair of nodes, $t_{ij}$, as a fuzzy variable. This variable has been modeled using a fuzzy inference system. We have assumed the case of electric vehicles and used their battery levels, as well as the reward of each node, as the input variables of the fuzzy system. The battery level (Q) of each vehicle can be estimated as low (L), medium (M), or high (H). The low and high levels are represented by a triangular fuzzy number $Q=(q_1,q_2,q_3)$, while the medium level follows a trapezoidal fuzzy number $Q=(q_1,q_2,q_3,q_4)$. All battery values are expressed as a proportion of the total battery level, i.e., $0\le Q\le 1$. The membership function of this fuzzy set is displayed in Figure 7. Similarly to the battery level, the reward of each node has been categorized using three fuzzy sets: low (L), medium (M), or high (H), where each of them follows a triangular distribution. The reward values have been represented as a proportion of the maximum reward that can be collected at any node of all the possible nodes to be visited.

Finally, the output of the fuzzy system gives a preference index, p, which indicates the inclination to visit the next node in the route. This index depends on both the reward of the next node and the remaining battery of the vehicle. This preference index has been defined between 0 and 1, i.e., $0\le p\le 1$. When $p=1$, the vehicle will definitely visit the next node in the route, since the vehicle will reach the node. On the contrary, when $p=0$, we are sure that the vehicle will not reach the next node, and the vehicle will remain in the current node. In this case, the route will present a failure, because the vehicle fails to reach the final depot, and consequently, the total reward of the route has been onlyt partially collected. We have classified the preference as: very low (VL), low (L), medium (M), high (H), or very high (VH). Each of these categories is represented by a fuzzy set. Finally, we have established a set of reasoning rules (Table 2), which describe the knowledge needed to determine the preference to visit the next node. After a quick fine-tuning process, we have set the threshold value for visiting the next node to $p\ge 0.45$. Note that this is a sensitive value, as a larger value could lead to generating overly conservative routes, while a value close to 0 could lead to risky decisions. In order to transform the input variables into a crisp value, the contribution of each membership function is combined on the inference, while a union operator is used to determine the output distribution. Subsequently, the center-of-gravity method is applied in order to obtain a crisp output value corresponding to the preference value.

6. Results and Discussion

Table 3 and Table 4 display the results of selected instances with different characteristics for the VRP and the TOP problems, respectively. In the case of the VRP and the TOP, the results—with the exception of the gap column—are measured in terms of distance and reward units, respectively. The first column of the tables identifies the instances. We have divided the remaining columns into three different parts. First, our best-found deterministic solutions (OBD) are presented (these solutions do not consider stochastic or fuzzy variables, they refer to the deterministic version of the problem). We compare the gap of our solutions (column 2) with respect to the best-known solutions (column 1). In the second part of the table, we present the obtained solutions for the stochastic scenario. Column 3 displays the expected cost when the OBD is evaluated under a stochastic scenario, with the corresponding level of uncertainty. To compute the expected cost, a simulation process has been applied to the OBD solution. Similarly, the next column shows the expected cost obtained using our simheuristic approach for the stochastic version of the problem. The last part of the table reports the results obtained considering fuzzy scenarios. Thus, column 5 reports the best hybrid fuzzy-stochastic solutions. To compute these solutions, we assume that half of the nodes follow a log-normal distribution, and the remaining half are considered to be fuzzy. In the case of the TOP, where the uncertainty is related to the edges, we have considered the origin node to evaluate the type of uncertainty. Finally, the last column of the table reports the solutions obtained in a scenario with a high level of uncertainty, where all the uncertain variables are considered as fuzzy.

Notice that, although the goal of this paper is not to solve the deterministic version of the problem, the results show that our approaches are highly competitive for the deterministic version of both problems. For the VRP problem, we obtain an average gap of $0.39\%$, with a maximum gap of $1.27$%. Furthermore, the obtained gap is $0.0\%$ for the TOP problem. These results highlight the quality of our base algorithms, which constitutes the optimization component in our fuzzy simheuristic, validating their potential to be used in uncertainty scenarios.

Regarding the uncertainty scenarios—which represent the main contribution of this paper—Figure 8a,b depict an overview of Table 3 and Table 4, respectively. In these box plots, the vertical axis represents the gap that was obtained in the stochastic and fuzzy scenarios with respect to the OBD solution—which is used here as a reference value. The latter can be considered as an ideal scenario with perfect information, which is not the case in scenarios with stochastic or fuzzy components. Concerning the stochastic solutions, the results show that those provided by the simheuristic clearly outperform the solutions of the deterministic version of the problem when these are simulated for all the considered problems, i.e., using the OBD solutions for the stochastic scenario is not a good idea, since it leads to sub-optimal solutions. On average, an improvement of about $7.91\%$ is observed for the VRP problem, while an improvement of about $1.72\%$ is reported for the TOP problem. These results justify the importance of using hybrid simulation-fuzzy methods when dealing with optimization problems under uncertainty.

Regarding the fuzzy scenarios, Figure 9 displays a summary of the presented results for different problems. The vertical axis represents the gap obtained for the different optimization methods with respect to the OBD solution. This figure shows that the solution quality worsens as the uncertainty level increases. This is due to route failures occurring during the execution stage, which penalize the entire route and, therefore, cause an extra cost. Figure 10, Figure 11, Figure 12 and Figure 13 illustrate a numerical example for the VRP instance A-n80-k10. Figure 10 depicts the configuration of the deterministic solution and its associated cost ($1797.05$). This cost can be seen as a lower-bound reference value in a scenario with perfect information. Figure 11, Figure 12 and Figure 13 show a representation of the obtained solution considering different levels of uncertainty. As the level of uncertainty increases, the total cost also increases, i.e., the highest cost ($1860.94$) is reached in the completely fuzzy scenario, where the solution cost has increased up to $3.43\%$ with respect to the deterministic solution. This extra cost is mainly caused by the rise in failure costs, since a greater number of detours and round-trips are expected when the uncertainty in the demand at each node is higher. Note also that the solutions have different configurations in each scenario, since the optimization algorithm attempts to generate routes that reduce the risk of failure.

7. Conclusions

This work has introduced the “fuzzy simheuristic” methodology to deal with NP-hard transportation problems under uncertainty scenarios, both probabilistic and fuzzy in nature. This uncertainty is tackled in a general way, since we consider that both stochastic and fuzzy uncertainty are present in many real-life transportation systems. Hence, pure deterministic, pure stochastic, and pure fuzzy scenarios represent particular cases that can also be addressed by employing our fuzzy simheuristic methodology. Since our methodology combines metaheuristics with stochastic and fuzzy simulation, it takes the best characteristics of both worlds, i.e., (i) the metaheuristics component provides the efficiency necessary to explore the solution space in order to find near-optimal solutions in short computational times. This characteristic becomes highly relevant when dealing with transportation problems, which are usually NP-hard; and (ii) the stochastic/fuzzy simulation component provides suitable tools to cope with different types of uncertainty, in order to provide high-quality solutions in terms of expected costs, expected profits, or risk/reliability indicators. A set of numerical instances of two well-known transportation topics—the vehicle routing problem (VRP) and the team orienteering problem (TOP)—demonstrates these advantages.

The simultaneous consideration of stochastic and fuzzy uncertainty arises whenever a subset of elements in a transportation problem–e.g., customers, roads, or vehicles–allows us to model some uncertainty aspects using probability distributions, while others require fuzzy techniques due to their vagueness or to the lack of enough historical data. The well-known VRP and TOP have been useful in testing our approach. For the VRP, we have studied a numerical example in which demands associated with a group of customers are stochastic, while a different group of customers presents fuzzy demands. Regarding the TOP, we have analyzed a case study in which travel times between customers are stochastic for a group of edges, and fuzzy for another group. The obtained results show that employing our approach leads to improve the solution quality—in terms of total cost for the VRP, and total collected reward for the TOP—when uncertainty is considered. All in all, these numerical examples illustrate the efficiency of the proposed methodology to solve transportation problems combining, at the same time, deterministic, stochastic, and fuzzy elements, something that has been rarely explored in the existing literature despite its relevance in real-life applications.

The challenges for future work are huge, given the increasing necessity of providing agile and very good solutions to real-world transportation problems considering uncertainty. In this sense, richer versions of the VRP and the TOP can employ and adapt our fuzzy simheuristic approach. This can be done, for instance, by including decisions about inventory management, multiple depots, facility location, or time windows, as well as by incorporating external dynamic conditions that force us to constantly update the routing plans.

Author Contributions

Conceptualization, A.A.J. and J.P.; methodology, A.A.J., J.P., P.J.C. and R.D.T.; validation, P.J.C., J.P. and R.D.T.; formal analysis, R.D.T., E.P.-B. and L.d.C.M.; writing—original draft preparation, R.D.T., L.d.C.M. and E.P.-B.; writing—review and editing, A.A.J., R.D.T. and E.P.-B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

This work has been partially supported by the Spanish Ministry of Science (PID2019-111100RB-C21/AEI/10.13039/501100011033, RED2018-102642-T), and the Erasmus+ program (2019-I-ES01-KA103-062602).

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Network topology of a basic VRP.

Figure 2. Network topology of a basic TOP.

Figure 3. High-level flowchart of the proposed solution method.

Figure 4. Fuzzy sets for the customer i demand.

Figure 5. Fuzzy sets for the remaining capacity after visiting customer i.

Figure 6. Fuzzy sets for the preference strength to travel to customer i.

Figure 7. Fuzzy sets for the battery of each vehicle.

Figure 8. Gaps of different optimization methods with respect to the OBD solution. (a) Results forthe VRP dataset. (b) Results for the TOP dataset.

Figure 9. Gaps of different optimisation methods with respect to the OBD solution.

Figure 10. Best solution for VRP-Deterministic scenario.

Figure 11. Best solution for VRP-Stochastic scenario.

Figure 12. Best solution for VRP-Stochastic and Fuzzy scenario.

Figure 13. Best solution for VRP-Fuzzy scenario.

Table 1

The rules used in the Fuzzy system for the VRP.

Demand	Remaining Capacity
	L	M	H
L	M	H	VH
M	L	M	H
H	VL	L	M

Table 2

The rules used in the Fuzzy system for the TOP.

Battery	Reward
	L	M	H
L	VL	L	M
M	L	M	H
H	M	H	VH

Table 3

Comparison of results, in terms of traveled distance, for the different VRP scenarios.

Instance	Deterministic Scenario			Stochastic Scenario		Fuzzy Scenario
	BKS	OBD	GAP (%)	Det	Stochastic	Stoch-Fuzzy	Fuzzy
	(1)	Sol. (2)	(1)–(2)	Sol. (3)	Sol. (4)	Sol. (5)	Sol. (6)
A-n32-k5	787.1	787.2	0.01%	1117.0	797.5	1119.0	1245.4
A-n33-k5	662.1	662.1	0.00%	895.4	667.8	811.6	988.8
A-n33-k6	742.7	742.7	0.00%	860.6	755.0	868.6	990.8
A-n37-k5	672.5	674.2	0.26%	974.8	682.4	926.8	1044.3
A-n38-k5	733.9	739.7	0.79%	768.4	761.5	945.0	1118.5
A-n39-k6	833.2	835.2	0.24%	835.7	835.3	1118.1	1236.4
A-n45-k6	952.2	957.1	0.51%	1074.2	1030.0	1264.9	1510.4
A-n45-k7	1147.4	1156.4	0.79%	1251.6	1161.2	1438.8	1754.7
A-n55-k9	1074.5	1085.9	1.06%	1128.6	1126.9	1292.8	1584.8
A-n60-k9	1360.6	1365.8	0.38%	1380.4	1371.5	1791.1	2135.4
A-n61-k9	1040.3	1049.0	0.83%	1128.4	1118.6	1361.2	1622.5
A-n63-k9	1633.7	1641.0	0.45%	1720.9	1689.4	2248.6	2847.1
A-n65-k9	1184.7	1195.2	0.89%	1332.2	1271.5	1594.6	1936.4
A-n80-k10	1776.2	1792.7	0.93%	2013.0	1838.5	2761.4	3334.5
B-n31-k5	676.1	676.5	0.06%	697.5	677.3	693.2	1068.4
B-n35-k5	958.9	959.4	0.05%	1053.4	1033.2	1215.4	1588.0
B-n39-k5	553.2	553.7	0.08%	585.6	584.6	855.3	1037.3
B-n41-k6	835.8	840.8	0.60%	924.9	862.8	1152.2	1344.5
B-n45-k5	754.0	754.7	0.10%	776.4	768.3	1085.0	1202.7
B-n50-k7	744.2	744.2	0.00%	789.4	778.1	991.7	1235.8
B-n52-k7	754.5	756.8	0.31%	852.1	763.8	1084.5	1355.6
B-n56-k7	716.4	719.4	0.42%	802.6	735.5	1023.0	1340.3
B-n57-k9	1602.3	1603.8	0.09%	1915.8	1700.0	2047.6	2723.9
B-n68-k9	1300.2	1306.5	0.48%	1491.7	1359.1	1776.7	2348.9
B-n78-k10	1250.6	1256.6	0.48%	1413.0	1383.9	1759.8	2146.6
E-n22-k4	375.3	375.3	0.00%	375.3	375.3	376.4	502.3
E-n30-k3	505.0	505.0	0.00%	505.0	507.9	742.4	838.0
E-n33-k4	837.7	839.4	0.21%	839.7	839.7	1183.0	1506.9
E-n76-k10	841.3	852.0	1.27%	926.0	902.9	1194.6	1336.8
Average	941.6	945.8	0.39%	1049.3	978.6	1266.3	1549.2

Table 4

Comparison of results, in terms of collected reward, for the different TOP scenarios.

Instance	Deterministic Scenario			Stochastic Scenario		Fuzzy Scenario
	BKS	OBD	GAP (%)	Det	Stochastic	Stoch-Fuzzy	Fuzzy
	(1)	Sol. (2)	(1)–(2)	Sol. (3)	Sol. (4)	Sol (5)	Sol. (6)
p1.2.f	80	80	0.0	78.9	79.3	77.4	76.8
p1.2.i	135	135	0.0	127.6	129.4	125.2	117.9
p1.2.k	175	175	0.0	169.3	174.4	164.6	143.2
p1.2.n	235	235	0.0	228.8	232.7	216.2	189.7
p1.3.n	190	190	0.0	182.5	189.6	180.0	174.1
p1.4.j	75	75	0.0	59.3	63.3	51.3	45.1
p1.4.k	100	100	0.0	98.3	99.9	99.9	95.8
p1.4.l	120	120	0.0	118.9	119.2	118.6	117.1
p1.4.m	130	130	0.0	98.2	102.9	91.5	86.3
p1.4.n	155	155	0.0	99.9	104.0	107.5	98.2
p1.4.o	165	165	0.0	159.4	164.2	148.6	143.6
p1.4.p	175	175	0.0	171.3	174.4	163.7	155.5
p2.2.d	160	160	0.0	150.6	150.6	150.0	137.5
p2.2.i	230	230	0.0	223.4	226.3	228.5	204.0
p2.3.i	200	200	0.0	191.5	195.2	186.7	177.9
p3.2.c	180	180	0.0	179.1	179.2	178.9	160.2
p3.2.d	220	220	0.0	212.3	217.5	197.7	179.8
p3.2.g	360	360	0.0	358.3	358.8	308.2	297.9
p3.2.q	760	760	0.0	748.5	755.2	663.4	630.6
p3.2.r	790	790	0.0	768.3	774.9	656.1	638.5
p3.3.e	200	200	0.0	198.2	199.0	195.7	187.6
p3.4.g	220	220	0.0	212.6	217.3	205.0	191.3
p5.2.d	80	80	0.0	75.5	77.4	73.7	70.7
p5.2.k	670	670	0.0	643.3	662.1	646.0	612.5
p5.2.p	1150	1150	0.0	1135.4	1138.1	1105.1	1073.0
p5.3.f	110	110	0.0	107.4	109.1	107.6	103.0
p5.3.o	870	870	0.0	856.2	865.1	836.9	806.9
p5.4.g	140	140	0.0	135.3	137.9	134.3	129.2
p5.4.t	1160	1160	0.0	1139.5	1148.4	1107.4	1068.1
p5.4.u	1300	1300	0.0	1279.5	1286.3	1239.2	1198.2
p6.2.d	192	192	0.0	185.4	188.1	177.8	164.2
p6.2.e	360	360	0.0	276.4	297.2	285.4	277.9
p6.2.f	588	588	0.0	577.4	580.0	519.5	501.2
p6.2.g	660	660	0.0	648.3	650.5	584.9	569.2
Average	362.8	362.8	0.0	349.9	354.3	333.3	318.3