# Vehicle Routing Problem with Deadline and Stochastic Service Times: Case of the Ice Cream Industry in Santiago City of Chile

Sebastián Dávila, Miguel Alfaro, Guillermo Fuertes, Manuel Vargas, Mauricio Camargo

► To cite this version:

Sebastián Dávila, Miguel Alfaro, Guillermo Fuertes, Manuel Vargas, Mauricio Camargo. Vehicle Routing Problem with Deadline and Stochastic Service Times: Case of the Ice Cream Industry in Santiago City of Chile. Mathematics, 2021, 9 (21), pp.2750. 10.3390/math9212750. hal-03411558

HAL Id: hal-03411558

https://hal.univ-lorraine.fr/hal-03411558

Submitted on 2 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Vehicle Routing Problem with Deadline and Stochastic Service Times: Case of the Ice Cream Industry in Santiago City of Chile

Sebastián Dávila<fnref n="1" /> ID, Miguel Alfaro<fnref n="1" /> ID, Guillermo Fuertes<fnref n="1" />,<fnref n="2" />,* ID, Manuel Vargas<fnref n="1" /> ID and Mauricio Camargo<fnref n="3" /> ID

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Dávila, S.; Alfaro, M.; Fuertes, G.; Vargas, M.; Camargo, M. *Vehicle Routing Problem with Deadline and Stochastic Service Times: Case of the Ice Cream Industry in Santiago City of Chile.* Mathematics **2021**, *9*, 2750. https://doi.org/10.3390/math9212750

Academic Editor: Poul Hjorth

Received: 24 August 2021

Accepted: 19 October 2021

Published: 29 October 2021

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

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** The research evaluates the vehicular routing problem for distributing refrigerated products. The mathematical model corresponds to the vehicle routing problem with hard time windows and a stochastic service time (VRPTW-ST) model applied in Santiago de Chile. For model optimization, we used tabu search, chaotic search and general algebraic modeling. The model's objective function is to minimize the total distance traveled and the number of vehicles using stochastic waiting restrictions at the customers' facilities. The experiments were implemented in ten scenarios by modifying the number of customers. Experiments were established with several customers that can be solved using the general algebraic modeling technique in order to validate the tabu search and the chaotic search methods. The study considered two algorithms modified with Monte Carlo (tabu search and chaotic search). Additionally, two modified algorithms, TSv2 and CSv2, were proposed to reduce execution time. These algorithms were modified by delaying the Monte Carlo procedure until the first set of sub-optimal routes were found. The results validate the metaheuristic chaotic search to solve the VRPTW-ST. The chaotic search method obtained a superior performance than the tabu search method when solving a real problem in a large city. Finally, the experiments demonstrated a direct relationship between the percentage of customers with stochastic waiting time and the model resolution time.

**Keywords:** supply chain management; optimization; vehicle routing problem; chaotic search; tabu search

## 1. Introduction

The cold chain is the key element in the transport and logistics of perishables. The cold chain consists of steps that are part of the refrigeration and freezing process that perishable products need in order to reach their final consumer in optimal conditions [1]. Food products often go bad due to long travel times and frequent stops to serve customers during the delivery process. Therefore, it is difficult to effectively manage cold chain distribution and ensure maximum freshness [2].

Logistics operators must have vehicles with temperature control, which are generally more expensive and consume more fuel than normal vehicles. The late delivery of perishable food significantly affects the costs of the logistics operator and retailers [3,4]. Some contracts with retailers involve penalties for delays within the established delivery times. Additionally, energy is required to maintain adequate temperature during distribution, and traffic conditions are different throughout the day [5]. Therefore, the impact of delays on delivering perishable food is more critical than general goods [6].

Requirements to serve customers within specific delivery time windows can increase the complexity of vehicle scheduling and routing problems for logistics operators.

<footnote n="1">¹ Industrial Engineering Department, University of Santiago de Chile, Avenida Ecuador 3769, Santiago 9170124, Chile; sebastian.davila@usach.cl (S.D.); miguel.alfaro@usach.cl (M.A.); manuel.vargasg@usach.cl (M.V.)</footnote>

<footnote n="2">² Facultad de Ingeniería, Ciencia y Tecnología, Universidad Bernardo O'Higgins, Avenida Viel 1497, Ruta 5 Sur, Santiago 8370993, Chile</footnote>

<footnote n="3">³ Université de Lorraine, ERPI, F-54000 Nancy, France; mauricio.camargo@univ-lorraine.fr</footnote>

<footnote>* Correspondence: guillermo.fuertes@usach.cl</footnote>
Therefore, vehicle routing problems (VRP) related to the delivery of goods have been stud-
ied extensively. VRPs are combinatorial optimization problems of high complexity; these
problems are categorized as NP-hard [7]. NP-hard models are widely used to solve opti-
mization models in supply chains [8,9], transportation [10,11], telecommunications [12,13],
cost optimization [14–16], manufacturing [17,18], service production [19], electrical net-
works [20], and social networks [21], among others. In the NP-hard class, processing time
and capacity increase exponentially as the number of customers increases. The need to
make decisions on real route assignment problems restricts processing times, and inves-
tigations focus on finding the best feasible solution. According to Toth and Vigo [22],
heuristics have been designed to narrow down search space, reducing processing times to
find feasible solutions.

Other solution methods are metaheuristics. Authors such as Gendreau and Potvin [23]
define metaheuristics as solutions which orchestrate an interaction between local improve-
ment procedures and high-level strategies, with processes that avoid local optimum by
performing a robust search for a solution space. The main challenge of these metaheuristics
is to solve large combinatorial optimization problems in an adequate time period. Different
authors have used metaheuristic algorithms to solve VRP: local search [24], simulated
annealing [25], greedy randomized adaptive search procedure (GRASP) [26], swarm in-
telligence [27], tabu search (TS) [28,29], genetic algorithms [30], colony optimization [31],
reactive search [32], and maximum coverage [33].

The problem analysis requires that each vehicle delivers the merchandise to the
customers within a specific time interval, known as a vehicle routing problem with time
windows (VRPTW). Considering Hsu et al. [34] and the randomness of the perishable food
delivery process, we propose using VRPTW-ST to generate delivery routes by minimizing
the total distance traveled and complying with client demand. This study aims to validate
the TS and chaotic search (CS) methods to minimize the total distance traveled and the
number of vehicles using stochastic waiting restrictions at the client’s facilities. The second
objective is to demonstrate the capacity of these techniques to reduce the execution times
of the computational model when establishing daily capacities and to generate new routes.

VRPTW-ST sets new restrictions on the VRP model; time windows limits (TW) are
constraints, and the service times (ST) are random variables. Other features of the general
VRPTW-ST model are that the vehicles can arrive early; however, the goods will still
be received within the time window at a random reception time. The complexity of
the VRPTW-ST model lies within the probability of vehicle arrival generated by the TW
restrictions and the ST variables [35].

This study uses VRPTW-SW in a frozen products distribution center of a transnational
food company located in the Quilicura commune in the Metropolitan Region of Santiago.
The distribution routes are planned daily to meet the demand of customers located in the
same region. For this case study, two customer segments are considered: (1) supermarkets
and (2) road customers. The former establishes a deadline for receiving orders, where
the truck is not serviced when it arrives after the scheduled time. Another restriction of
the supermarkets is the random waiting time of the trucks before being served. In the
case of road customers, they can receive the order at any time, and once the trucks arrive,
they immediately begin receiving the merchandise.

1.1. *Literature Review*

This section analyzes studies on VRP with stochastic service and/or travel times
and studies evaluating different sources of uncertainty. VRP is a classic combinatorial
optimization problem originally introduced by Dantzig and Ramser [36]. Taş et al. [37]
proposed a VRP model with a weak time window and stochastic travel time, comparing
the results with TS and an adaptive large neighborhood search. These solutions are useful
for large-scale problems, susceptible to rush hour.

Ehmke et al. [38] proposed using programming with probability restrictions in VRP to
guarantee a certain level of service for all customers. Stochastic programming must be used
to solve this vehicle routing problem. Xu et al. [39] used an improved hybrid ant colony optimization algorithm, K-means, 2-Opt, and crossover. The experimental results showed that the ant colony optimization algorithm is capable of finding high-quality solutions. Furthermore, Tao et al. [32] proposed a metaheuristic based on the hybrid topological graph, genetic algorithm, and TS to minimize travel and waiting times. The algorithm considers time windows, load capacity, and the origin of tasks. The results showed the efficiency and effectiveness of the proposed algorithm, while Urzúa-Morales et al. [33] used a VRP model for a merchandise distribution system in the historic center of the city of Santiago de Chile. The last mile modeling considered a maximum coverage optimization model, k-nearest neighbors, and the analytic hierarchy process. The results reduced 53 tons of carbon dioxide within the square kilometer and reduced 1103 h of interruptions per year in vehicle congestion.

Yu et al. [40] proposed two lower limit models that determine optimal quantity for the bottleneck process in production and distribution logistics. The distribution problem is modeled as a VRPTW. A Lagrangian relaxation model was designed to optimize the lower limit, and an improved subgradient algorithm was proposed. The results show that the suggested algorithm can calculate an adjusted lower limit. Laporte et al. [41] used programming with probability restrictions to assign a certain probability of failure or several restrictions to VRP problems.

Taking into account the third type of uncertainty (service times), Errico et al. [35] developed a VRPTW-ST model for route planning. The model includes a branch and cut-price algorithm and an adaptation of the dynamic programming algorithm to determine the probabilistic consumption of resources. Han et al. [42] presented a single-depot vehicle-routing problem with the time windows (SDVRPTW) model for attended home delivery. This study considers two types of behaviors: customer unavailability and random response time. The solution method proposes a hybrid algorithm which connects dynamic programming and TS heuristics. The proposed approach generates high-quality solutions within a reasonable execution time.

Zhang et al. [43] consider a stochastic VRP problem with uncertain deadlines (VRPD) to find profitable routes complying with the service level, subject to the uncertainty of the distribution. They considered a robust distribution optimization framework. The computational difficulties caused by the probability constraints are approximated using value-at-risk constraints. The results show that the routing policies developed from this framework are robust in all the scenarios considered. Authors such as Shukla et al. [44] use a portfolio of evolutionary algorithms to reduce computer time when searching for local solutions to resolve a vehicle routing problem with stochastic demand (VRPSD). Furthermore, Hand and Chu [45] proposed a multi-start variable neighborhood search approach for solving the split-delivery vehicle routing problem with minimum delivery amounts (SDVRP-MDA), in which the demand of a customer can be split and delivered to the customer by several vehicles on different routes, but the customers each require a minimum delivery amount every time they are visited. The results suggest that a relatively simple multi-start heuristic with good search operators can be as effective as sophisticated metaheuristics. Braaten et al. [46] considered a robust vehicle routing problem with time windows (RVRPTW) with applications in maritime transportation. Computational results show a considerable benefit from the main novel heuristic components.

## 1.2. *Research Gap and Contribution of This Study*

The random nature of the variables implicit in route planning (demand, travel times, service times, and physical availability of the client) makes decision making difficult. According to Ehmke et al. [38], there are two main problems in the planning of routes that cause a loss of customers and cost overruns: (1) dissatisfaction in service due to non-compliance with time windows and (2) cost overruns in the final sale price of the product due to overtime pay to drivers. Since 1992, random variables have been incorporated into
the VRP models, giving rise to the stochastic vehicle routing problem (SVRP) incorporating
random elements.

The main variants studied in recent years are the stochastic demand, when there is no certainty of the demand of each client [47], stochastic travel time [37], and stochastic clients with probability *p* that the client is found and a probability of 1 − *p* that the client is not found [37]. Other jobs that incorporate stochastic service time are referenced in [48].

The literature considers two types of time windows: (1) strict time window—if the vehicle arrives before the time window it must wait, in general, without an associated cost; (2) flexible time window, which allows violating the restriction assuming a penalty [22]. For Li et al. [49], high-value penalties greatly limit the probability of violating the restriction and result in adapting behaviors similar to the strict time window. If there is only an upper limit in the time window, the model is called a vehicle routing problem with a deadline (VRPD).

Finally, hybrid methods solve VRPTW-ST. The two-phase method solves a model
without the probability of failure to obtain a preliminary result and subsequently resolves a
second model that considers the costs associated with the failure of any of the restrictions.

The main contributions of this work are:

1. In studies [34,35], a waiting time is established for vehicles only when they arrive before the time window. Our research generalizes random wait times for all delivery vehicles arriving before or during the time window.

2. In the article, we develop an extension of the CS and TS method for the VRPTW-ST problem. This extension solves problems with reduced calculation times for the resolution of the model.

3. The proposed model is made up of clients with and without reception restrictions. The study presents different combinations of the confidence level $(\alpha, \beta, \gamma)$, maintaining the percentage of clients with restrictions. The different variations of the confidence levels modify the best value of $Z(x)$. Very high confidence levels decrease the solution quality.

1.3. *Problem Statement and Insights from This Study*

In this section, the problem described above was represented as a VRP with a dead-
line and stochastic waiting time (VRPTW-ST) for a real distribution problem of frozen
food for a transnational company in the city of Santiago. The vehicle is serviced within
the attention intervals defined by the client, and the waiting time is defined. Based on
previous works [50], we propose a mathematical model with probabilistic restrictions that
considers the random nature of the problem when establishing a level of service for each
customer [38].

To address the complexity of the VRPTW-ST, Li et al. [49] developed a metaheuristic
TS-based model, and Ikeguchi et al. [50] used a CS method, obtaining better results than
the TS method.

To validate the proposed methodology, models were constructed to compare both
methods with the empirical results of the industry. The solution methods were programmed
in the Python language.

The data used were (i) customer demand, (ii) average travel times and distances between customers, and (iii) average and variance of waiting time for each customer. The proposed modifications to the TSv2 and CSv2 algorithms free up memory and reduce numerical computation. These new algorithms do not use the Monte Carlo (MC) method to choose the new best move. Subsequently, the algorithm recalculates the best value of Z(x) using the MC method. These modifications allow an efficient optimization design of the vehicle route graph. The program configuration was designed to run all four codes (TS, CS, TSv2, and CSv2). For each code, we created ten experiments, and the results per code correspond to the average of these ten experiments.

The novelty of this research is the development of a metaheuristic that solves real
VRPTW-ST to distribute merchandise in large cities. These problems cannot be solved
quickly by the general algebraic modeling systems (GAMS) technique. Here, we proposed the use of CS and TS heuristics to solve the problem of time and computational capacity.

The insights of the study correspond (1) to the design of an optimization strategy to solve VRPTW-ST problems associated with the supply chain of frozen products, (2) the techniques TS and CS modified with MC reduce execution times and program the new vehicle routing daily, (3) the optimization strategy can adapt to new VRPTW problems associated with large cities, (4) the modified CS technique provides adequate resolution times for customer volumes higher than 80, and (5) the implementation of the algorithms reduces overtime for drivers, reducing costs. This reduction in labor costs and the reduced processing times of the different scenarios allows commercial implementation of the solution.

The development of this study is divided into four sections. Section 1 reviews the literature concerning the main concepts and heuristics that have served as solution methods to different VRPs with time windows and stochastic service. Section 2 describes the methodology. Section 3 presents the results of the proposed VRPTW-ST. Finally, Section 4 exposes the main conclusions of the proposed investigation.

## 2. Materials and Methods

### 2.1. Previous Mathematical Approach

The VRP is defined according to the graph theory as a unidirectional graph represented by $G = (V, E)$. Where $V = (0,1,2,...,n+1)$ is the set of vertices with $n + 1$ representing a depot, $E = \{(i,j) : i, j \in V_{i \neq j}\}$ is the set of arcs, $V_0 = V \setminus \{0, n + 1\}$ denotes the set of customers, 0 represents the depot where a fleet of vehicles with identical capacity $Q$ to serve customers. Each customer $i$ has a fixed demand $q_i$, which is delivered from the depot in a service time $\delta_i$. Each arc $(i, j)$ is associated with a travel distance or a cost $c_{ij}$. In general, it is assumed that $C = (c_{ij})$ satisfies the triangular inequality, that is, $c_{ij} \le c_{ik} + c_{kj}$ for every $i, j, k \in V$.

The VRPTW-ST model probability restrictions are based on the author's work [51] and assign service levels to comply with the probability restrictions.

### 2.2. Set, Parameters and Variable

The model is composed of the following sets, parameters, and variables:

$m$ : number of vehicles used, can be a constant or a decision variable.

$n_k$ : number of customers served by the vehicle $k$.

$K$ : set of vehicles, which is defined as $K = \{1,2,...,m\}$.

$Q$ : vehicle capacity; assumed that all the vehicles are homogeneous and each one has limited capacity $Q = 3000$.

$B$ : established working time of the driver. In this case, 480 min. The value of $B$ is composed of travel time, service time, and waiting time.

$W$ : the maximum time that can be spent waiting by every vehicle.

$q_i$ : customer demand $i \in V_0$.

$d_{ij}$ : distance between vertex $i$ and $j$ with $i, j \in V$.

$t_{ij}$ : travel time between vertex $i$ and $j$.

$\delta_i$ : customer service time $i$.

$w_i$ : customer timeout $i$; in this work, $w_i$ is assumed as continuous random variable with normal distribution $w_i \sim N(\mu_i, \sigma_i^2)$.

$s_{ik}$ : time for vehicle $k$ to visit the customer $i$; $s_{ik}$ is a continuous random variable because it depends on customer timeout.

$A$ : maximum waiting time for vehicle $k$.

Model variables:

$$x_{ijk} = \begin{cases} 1 & \text{if the edge } (i,j) \text{ is routed to truck } k, \ 0 & \text{otherwise} \\ k \in K, y, j \in V_0 \end{cases}$$
## 2.3. Objective Function

The probability constrained programming model of the vehicle routing problem with stochastic lead time, and the delivery deadline is shown below.

$$ \min F_0(x) = M \sum_{j \in V_0} \sum_{k \in K} x_{0jk} + \sum_{i \in V} \sum_{j \in V} \sum_{k \in K} d_{ij} x_{ijk} \quad (1) $$

Subject to:

$$ \sum_{k \in K} \sum_{j \in V} x_{ijk} = 1 \qquad \forall i \in V \setminus \{0, n+1\} \quad (2) $$

$$ \sum_{j \in V} x_{0jk} = 1 \qquad \forall k \in K \quad (3) $$

$$ \sum_{i \in V} x_{in+1k} = 1 \qquad \forall k \in K \quad (4) $$

$$ \sum_{i \in V_0} x_{ijk} - \sum_{i \in V_0} x_{jik} = 0 \qquad k \in K, j \in V_0 \quad (5) $$

$$ \sum_{i \in V_0} q_i \sum_{j \in V} x_{ijk} \le Q \qquad \forall k \in K \quad (6) $$

$$ s_{jk} - s_{ik} \ge t_{ijk} + \delta_i + w_i - M(1 - x_{ijk}) \qquad \forall k \in V\{0, n+1\}, k \in K \quad (7) $$

$$ P\left(\sum_{i \in V} \sum_{j \in V} t_{ij} x_{ijk} + \sum_{i \in V_0} (\delta_i + w_i) \sum_{j \in V} x_{ijk} \le B\right) \ge \beta \qquad \forall k \in K \quad (8) $$

$$ \left\{ \begin{array}{l} P(s_{ik} \le l_i) \ge \alpha, \\ P(s_{n+1k} \le l_{n+1}) = 1 \end{array} \right. \qquad i \in V_0, \qquad \forall k \in K \quad (9) $$

$$ P\left(\sum_{i \in V_0} \sum_{j \in V_0} w_i x_{ijk} \le W\right) \ge \gamma \qquad \forall k \in K \quad (10) $$

$$ x_{ijk} \in \{0, 1\} \qquad \forall (i, j) \in V, k \in K \quad (11) $$

$$ s_{ik} \ge 0 \qquad \forall i \in V, k \in K \quad (12) $$

Equation (1) represents the objective function, where parameter M is a large value to minimize the number of vehicles and the total distance traveled. Equation (2) represents how each customer can only be visited once by a vehicle. Equations (3) and (4) indicate that vehicle k must leave and return to the depot. Restriction (5) forces each vehicle k to leave the node. Equation (6) represents the capacity constraint. Equations (7) and (8) are the stochastic restrictions of driving time and waiting time of vehicles using the decision variable $x_{ijk}$. Equation (9) represents the probability of arriving before the deadline ($l_{n+1}$) with a level of confidence $\alpha$. Equation (10) represents the maximum waiting time of each vehicle. The model prioritized the customer to reduce the difficulty of finding a solution. Finally, constraint (11) does not allow partial services and ensures non-negativity of variables.

## 2.4. Stochastic Restrictions

Equation (8) establishes the duration of a driver's route with a B value under a probability P and a confidence level $\beta$. Excess in driver hours has a probability of $1 - \beta$. For very low $\beta$ values the probability of overtime is increased.

Equation (9) establishes the probability P with confidence level $\alpha$ of the vehicle arriving before the deadline. The random $s_{ik}$ variable depends on the travel time, the service time, the random waiting time, and the deadline. The confidence level $\alpha$ adjusts the problem's difficulty level; no feasible solutions may be found for very high values.

Given the high randomness of the waiting time, this constrains the decision-maker to select adequate confidence levels. Therefore, it is often an unnecessary or high cost to hold all restrictions. Firms can establish lower values for the constraint with less relevance than a higher significance value for the limitations.

## 2.5. Improvement Stochastic Restrictions

Equation (10) is a new restriction to the original model proposed by [49]. It restricts the extra times of each vehicle by adding all the waiting times of each vehicle W under a $\gamma$ probability.
Equation (8) is transformed into several restrictions where a certain number of possible scenarios are evaluated. The purpose of the transformation is to meet these restrictions a minimum number of times, and the compliance number will depend on the level of trust in the original restriction.

The set of scenarios is defined as $s = (s_1, \dots, s_{scen})$, where $s_{scen}$ is the number of scenarios. The $w_i$ variable must be redefined by $w_{is}$ to represent the customer's wait time $i$. The variable $y_s$ is a binary variable with a value of 1 when the constraint (8) is satisfied and 0 otherwise. Meanwhile, $\pi_s$ is the probability of occurrence of the scenarios; it is assumed that all the scenarios have the same probability. Finally, $cc_1$ is the probability of failure of the constraint (8), which must be greater than zero and less than probability P.

Equation (8) can be written as the following set of restrictions.

$$ \sum_{i \in V_0} \sum_{j \in V_0} w_{is} x_{ijk} \leq W + M(1 - y_s) \quad \forall k \in K, \forall s \in S \quad (13) $$

$$ \sum_{i \in V_0} \sum_{j \in V_0} w_{is} x_{ijk} \leq W - My_s \quad \forall k \in K, \forall s \in S \quad (14) $$

$$ cc_1 = 1 - \sum_{s \in S} \pi^s y_s \quad \pi^s = \frac{1}{card(S)} \quad (15) $$

$$ 0 \leq cc_1 \leq (1 - \gamma) \quad (16) $$

$$ y_s \in (0,1) \quad \forall s \in S \quad (17) $$

Finally, the model plans a set of routes subject to probability restrictions. Such scenarios allow routes not satisfied under a certain threshold. However, it does not consider the costs associated with the failure of routes. In this case, the objective is to minimize the number of vehicles and the total distance traveled; for Sawik [52], solutions with fewer vehicles are better than solutions with less total distance traveled.

## 2.6. Programming of Tabu Search

An initial solution for VRPWT-ST must be available for TS programming. The shortest path from each node must be adapted: the minimum distance between the customers and depot. This solution is built considering a greedy algorithm. Initially, we searched the minimum distance between the depot and customer, which was then repeated with each node until violating any restrictions. To build the solution, we used the upper wait time confidence interval with an await time probability of 0.9999. This level of confidence guarantees a workable initial solution.

The restrictions are addressed as follows: (1) when the next customer exceeds the maximum tuck capacity, then a new route is generated; (2) when the truck arrives after the deadline for delivery, then the delivery is postponed until the next route is generated and another node is analyzed, and (3) if the route time exceeds the working time $B (=480)$ then a new route is generated.

Each visited customer is assigned a large-value datum point to avoid being reassigned, eliminating the probability of partial deliveries. Subsequently, the possible violations of the restrictions of the original problem are checked, penalizing the objective function $F(x)$ of the model. This generates a new $F(x)$ represented in the following way (Equation (18)):

$$ Z(x) = F(x) + \sum_{k=1} (P_1 L_k + P_2 G_k + P_3 H_k + P_4 J_k) \quad (18) $$

where $F(x)$ is the original objective function defined in Equation (1). $L_k$ represents the feasibility of the capacity constraint. $G_k$, $H_k$, and $J_k$ represent the feasibility of the probability constraint: deadline, overtime and total wait time, respectively. $L_k$, $G_k$, $H_k$, and $J_k$ take a value of 1 in case the respective restriction is violated and 0 in the other case. We carried out three movements to build a new neighbor (i) exchange, (ii) relocation, and (iii) the 2-opt movement in each iteration. The first movement (i) switches the position between two random clients, the second movement (ii) transfers a random client to another route, and the third 2-opt movement fixes two random clients of the same route and reverses the order of the intermediate customers. Algorithm 1 summarizes the proposed TS heuristic.
**Algorithm 1.** TS heuristic.

Set *x* as the given initial solution

Update statistical data: set *x*<sup>*</sup>: *x* and update *F*<sup>*</sup>(*x*) and *Z*<sup>*</sup>(*x*). Allow that every neuron shot in the first iteration *X* = 1

while *t* ≤ *t*<sub>max</sub> and *t*2 ≤ *t*1

Set *n* = 0

while *n* ≤ *n*<sub>max</sub>

Set Z(x<sub>best</sub>) = ∞

Randomly relocated, exchange or two-opt operator, and two random customers (i, j) such constraints for each operator (e.g., the route must differ for relocated customers). Move current solution x<sub>crnt</sub> to its neighboring x<sub>iter</sub>, Set F<sub>-</sub>(x<sub>iter</sub>) and Z(x<sub>iter</sub>).

If and Z(x<sub>best</sub>) > Z(x<sub>iter</sub>) and (i, j) does not tabu then x<sub>best</sub> = x<sub>iter</sub>, F(x<sub>best</sub>) = F(x<sub>iter</sub>), Z(x<sub>best</sub>) = Z(x<sub>iter</sub>),
(i, j)<sub>best</sub> = (i, j) otherwise if x is feasible and F<sup>*</sup>(x<sup>*</sup>) > F(x<sub>iter</sub>) then x<sub>best</sub> = x<sub>iter</sub>, F(x<sub>best</sub>) = F(x<sub>iter</sub>),
Z(x<sub>best</sub>) = Z(x<sub>iter</sub>), (i, j)<sub>best</sub> = (i, j)

Update list tabu: add (i, j)_{best} for random value ~ U(5, 15) iterations, and reduce the rest in one.

If and Z<sub>best</sub> > Z(x<sub>i,k</sub>) then Z<sub>best</sub> = Z(x<sub>i,k</sub>) and x = x<sub>i,k</sub> If Z<sup>*</sup>(x<sup>*</sup>) > Z<sub>best</sub>(x) then Z<sup>*</sup>(x<sup>*</sup>) = Z<sub>best</sub>(x), Set t2 = 0 and t3 = 0

If *x* is feasible and *F*<sup>*</sup>(*x*<sup>*</sup>) > *F*(*x*), *x*<sup>*</sup> = *x*, and *F*<sup>*</sup>(*x*<sup>*</sup>) = *F*(*x*). Set *t*2 = 0 and *t*3 = 0. If *t*2 = *t*2 + 1 and *t*3 = *t*3 + 1

if *x*<sup>*</sup> is not updated for √*t*1 iterations then *x*<sub>best</sub> = *x*<sup>*</sup> and *t*3 = 0

If Z<sup>*</sup>(x<sup>*</sup>) > Z<sub>best</sub>(x) then Z<sup>*</sup>(x<sup>*</sup>) = Z<sub>best</sub>(x)

*t* = *t* + 1

2.7. Programming of Chaotic Search

Based on the research of Kato et al. [53], the algorithm uses two local searches:
(1) customer change and (2) customer relocation. The method requires 2n neurons to
solve a problem with n customers. When a neuron is triggered, the effect over the search is
a change of customer or customer relocation. The variables and the parameters to use in
the method are the following:

*α* : scalar parameter of tabu effect.

β : scalar parameter of the effect gained.

θ : bias positive.

n : number of customers.

m : number of local searches. In this case, there are only 2.

k<sub>r</sub> : decay parameters for tabu effects, k<sub>r</sub> ∈ {0, 1}.

$\varphi_{ij}(t)$ : internal state of the (i, j)-th neuron in the period $t$ corresponding to the effect gained, where $i = 1, \dots, n$ and $m = 1 \text{ or } 2$.

$\omega_{ij}(t)$ : internal state of the (i, j)-th neuron in the period $t$ corresponding to the refractory effect.

Each local search previously selects the best possible movement, that is, the movement which decreases the value of F(x). The movement is performed when the most efficient neuron is fired. Each neuron has a gained effect and a refractory effect described below:

$$
\text{Effect gained:} \qquad \varphi_{ij}(t+1) = \beta \max_s \{ \Delta_{ijs}(t) \} \qquad (19)
$$

The effect gained refers to the profit obtained when performing the local search; in
this case, it corresponds to the greatest savings. Variable Δ<sub>ijs</sub>(t) is obtained by the change
or relocation of the customer *i* with the customer *s* due to the local search *j*.

$$
\text{Refractory effect:} \qquad \omega_{ij}(t+1) = -\alpha \sum_{d=0}^{t} k_r^d x_{ij}(t-d) + \theta \qquad (20)
$$

The refractory effect is used to avoid falling into local optimum. This element performs
a job similar to memory in the TS method. The last shot is memorized as a previous state
when the intensity of the refractory effect of the current state is decided. The value increases
only after the neuron is fired and recovers exponentially over time, depending on the value
of *k*<sub>*r*</sub> and *θ*. This does not escape possible local minimums.

The trigger functions used in the CS are Equations (21)–(23). Equation (22) is f(y)
with sigmoid function and ε is a very small value. The sigmoid function generates more
complex shots.
For the condition $x_{ij}(t) > \frac{1}{2}$, the ij-th neuron is triggered in the time t, and the corresponding local search is performed; that is, if the n1-th neuron is fired, the n-th customer is changed, and if the n2-th neuron is fired, then the mn-th neuron is relocated.

Neurons update asynchronously, and the iteration ends when all neurons are updated. As in the TS, if iterations do not improve the solution, we propose intensifying the search by returning to the best solution found so far. Algorithm 2 summarizes the proposed CS heuristic.

$$x_{ij}(t+1) = f(y_{ij}(t+1)) \qquad (21)$$

$$f(y) = \frac{1}{1 + e^{-\frac{y}{\bar{\epsilon}}}} \qquad (22)$$

$$y_{ij}(t+1) = q_{ij}(t+1) + \omega_{ij}(t+1) \qquad (23)$$

**Algorithm 2.** CS heuristic.

Set x as the given initial solution

Update statistical data: set $x^*$: x and update $F^*(x)$ and $Z^*(x)$. Allow every neuron shot in the first iteration $X = 1$ while $t \le t_{max}$ and $t2 \le t1$

$$Z(x_{best}) = \infty$$

for (i,k) in neuron:

Move current solution $x_{best}$ to its neighboring $x_{i,k}$, Set $F_-(x_-(i,k))$ and $Z(x_{i,k})$ for CS calculation $Z(x_{i,k})$ used MC.

Calculate Refractory and Gain. Set y. Update $X_{i,k}(t+1) = f(y)$

If $X_{i,k}(t) \ge$ shot the neuron (i,k) and $Z_{best} > Z(x_{i,k})$ then $Z_{best} = Z(x_{i,k})$ and $x = x_{i,k}$ (for CS v2 update $Z_{best}$ with MC). If $Z^*(x^*) > Z_{best}(x)$ then $Z^*(x^*) = Z_{best}(x)$

Else If x is feasible and $F^*(x^*) > F(x)$, $x^* = x$, and $F^*(x^*) = F(x)$. Set $t2 = 0$ and $t3 = 0$. Else $t2 = t2 + 1$ and $t3 = t3 + 1$. Update of values of $\tau_i$ ($i = 1,2,3,4$)

If $x^*$ is not updated for $\sqrt{t1}$ iterations then $x_{best} = x^*$ and $t3 = 0$

If $\text{mod}(t,15) = 0$, $P_i$ ($i = 1,2,3,4$) is updated as follows: $P_i \leftarrow P_i \cdot 2^{(1 - \frac{2\tau_i}{T})}$, set $\tau_i = 0$ ($i = 1,2,3,4$) $t = t + 1$

## 2.8. Monte Carlo and Penalization Procedure

Both TS and CS were initially built for VRP deterministic cases. The authors of [49] adapted the TS to apply in stochastic cases. They implemented an MC simulation to validate the feasibility of chance constraints, providing a maximum number of runs ($N_0$), any chance constraint (e.g., (10)), and confidence level (e.g., $\gamma$). In each simulation, waiting time w is set randomly from a normal distribution $N(\bar{w}, \sigma_w)$ where $\bar{w}$ and $\sigma_w$ are the mean and standard deviation of the waiting time, respectively. Later, it is verified whether it holds the constraint. If the percentage of simulation it is greater or equal to the level of confidence, then the constraint is not violated; otherwise, a penalty is added to Z. Additionally, we updated the penalty value $P_i$, $i = 1,2,3,4$ after a certain number of the iteration T (in this case 15) to allow for diversification. Hence, the $P_i$ is multiplied for $2(1 - \frac{2\tau_i}{T})$, where $\tau_i$ is the number of instances that satisfy the constraint i. Notice that if the $\tau_i$ is small (the constraint has been violated in the last iterations), then the penalty increases with the proposed force to satisfy the successive 15 iterations, otherwise it is reduced to acceded new solutions. These procedures were considered for both TS and CS. It embeds these procedures into CS, extending its usefulness to stochastic problems. However, notice that CS for every (i, j) neuron must search the best move possible, which impacts its performance. For this reason, we considered only two operator exchanges for CS, and relocated.

Two modifications of the TS shown in [49] were proposed. They proposed MC to prove the feasibility of the chance constraints, calculate the value of $Z(x)$, and repeat it for every neighboring candidate generated. However, the runtime depended drastically on the number of candidates and the number of simulations established. Therefore, we proposed an alternative not to run MC proceedings while generating the candidates, but after selecting the best move to reduce the runtime. We used the average to random variable
w to calculate the best value of Z(x), and we updated its value using MC when finalizing this procedure. Additionally, for the neighborhood structure, we chose a pair of customers (i, j) randomly to generate an admissible solution. After, we selected an operator randomly to relocate the customers in a different route. The heuristics with this consideration are called TSv2 and CSv2.

## 3. Results

### 3.1. Results of the Previous Mathematical Approach

The investigation corresponds to the resolution of a real VRPTW-ST product transport problem in the City of Santiago. To validate the implementation of the TS and CS methods, they were compared with the authors' works [51]. The problem with probability restrictions could only be resolved by a small number with the GAMS mathematical programming module using exact methods. The programming with probability restrictions is resolved by evaluating the same model in different scenarios for computational downtime of 6 h, obtaining a maximum resolved size of 20 nodes.

Table 1 shows the results obtained for the instances with parameters $\alpha, \beta$, of 0.75. The effectiveness of the TS method for the solution of the proposed problem was evaluated. The results obtained in the instances that could be resolved through GAMS were compared with those obtained by TS. Table 1 shows the comparative results between GAMS-TS and GAMS-CS.

<table><caption>**Table 1.** Results GAMS and TS-CS.</caption><thead><tr><th>Customers</th><th>%SM</th><th>GAMS</th><th>Global Optimum Value</th><th>TS</th><th>Difference between GAMS-TS</th><th>CS</th><th>Difference between GAMS-CS</th></tr></thead><tbody><tr><td>10</td><td>10</td><td>10,081.42</td><td>0.01%</td><td>10,081.42</td><td>0.00%</td><td>10,081.65</td><td>0.00%</td></tr><tr><td>10</td><td>20</td><td>10,081.42</td><td>0.03%</td><td>10,084.82</td><td>0.03%</td><td>10,087.37</td><td>0.06%</td></tr><tr><td>10</td><td>30</td><td>20,118.34</td><td>0.04%</td><td>20,132.63</td><td>0.07%</td><td>20,130.64</td><td>0.06%</td></tr><tr><td>20</td><td>10</td><td>20,128.23</td><td>0.02%</td><td>20,155.18</td><td>0.13%</td><td>20,134.30</td><td>0.03%</td></tr><tr><td>20</td><td>20</td><td>20,133.47</td><td>0.04%</td><td>20,141.66</td><td>0.04%</td><td>20,130.00</td><td>0.02%</td></tr><tr><td>20</td><td>30</td><td>30,179.48</td><td>0.05%</td><td>30,180.95</td><td>0.00%</td><td>30,172.53</td><td>0.02%</td></tr></tbody></table>

The results have the expected behavior, the increase in customers, and/or an increase in the stochastic elements of the problem increases the objective function F(x). The results of the GAMS module failed to obtain the optimal solution; however, the difference concerning the optimal value is negligible. The resolution model is validated to solve this problem.

The difference obtained between both solutions GAMS and TS is similar. With TS, feasible results are obtained from the mathematical model proposed in Section 2.

The results of the metaheuristic CS confirm the validity of the technique by comparing the results with the GAMS technique for the problem posed in Section 2. In this case study, the main challenge was delivery to supermarkets because they only received orders until midday. Additionally, Table 1 shows that the mean waiting time was between 52 and 63 min, with a standard deviation of 33 and 53 min, depending on the brand.

Figures 1–3 show the average methods. The results showed that there is no clear advantage between heuristics. For instance, CSv2 obtained better results when alpha, beta, and gamma values were lower and worse when they were higher. CS and TS are the best results, with CS being the best in most instances.

Figure 4 shows the waiting time distribution for all brands (Table 2). Figure 5 shows the waiting time distribution as a function of $\sqrt{x}$. The waiting time distribution is almost normally distributed, validating the probability distribution suggested in Section 2. Historically, 10% to 60% of a supermarket's total customers visit every day. We recommend a deterministic mean service time between 15 and 25 min for road customers and supermarkets, respectively.
<figure><img src="image_4.png" /><figcaption>Figure 1. Total distance by methods. $\alpha = 0.5$ and supermarket = 60%.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2. Total distance by methods. $\alpha = 0.75$ and supermarket = 60%.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 3. Total distance by methods. A = 0.99 and supermarket = 60%.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4. Histogram of supermarket waiting time.</figcaption></figure>
<table><caption>**Table 2.** Statistical data of supermarket-type customers.</caption><thead><tr><th>Supermarket</th><th>Mean</th><th>StDev</th></tr></thead><tbody><tr><td>Brand 1</td><td>57.14</td><td>36.83</td></tr><tr><td>Brand 2</td><td>75.85</td><td>53.19</td></tr><tr><td>Brand 3</td><td>62.78</td><td>39.81</td></tr><tr><td>Brand 4</td><td>55.75</td><td>45.88</td></tr><tr><td>Brand 5</td><td>57.42</td><td>38.75</td></tr><tr><td>Brand 6</td><td>52.99</td><td>37.15</td></tr><tr><td>Total</td><td>57.72</td><td>37.75</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Histogram applying transformation $\sqrt{x}$ to supermarket waiting time.</figcaption></figure>

Meanwhile, we used historical data to establish the 80 customers' demand, mean wait time, travel times, and distances randomly selected, 10, 20, 30, ... 80, and then set the supermarket variable at 20%, 40% and 60% and randomly selected them again. We only selected supermarkets with random waiting times and set a deadline of half the working time (0.5·B). We only used the average waiting time for road customers. To study how the confidence level (α, β, γ) influenced chance constraints, we used instances of 0.5, 0.75, and 0.99 for each confidence level, i.e., α = 0.5, 0.75 and 0.99, β = 0.5, 0.75, and 0.99, and γ = 0.5, 0.75, and 0.99. Additionally, we used the highest value for the standard deviation of the random variable w. We set the rest of the parameters as B = 480, W = 90, and Q = 1000. We set the parameters of the TS $t_{max}$ = 2000, $t1$ = 3·|V₀|, the CS to $t_{max}$ = 1000, $t1$ = 3·|V₀|, and established 100 simulations for every iteration of the MC proceedings. We ran 10 instances with different seeds for each set of parameters. Likewise, we incorporated the maximum runtime as a stop criterion set in 3600 min.

The resolution was made with TS and CS with their respective version 2. Tables 3, 4 and S1-S4 show the result for the average distance, vehicles, and execution time; the results establish a direct relationship between the increase in the percentage of stochastic customers and the number of customers with the increased distances. Likewise, the number of trucks and distances increased with higher confidence levels, alpha (associated with the deadline constraint) being the one which harms the objective function the most.
<table><caption>**Table 3.** Summary of the results obtained by the company data using GAMS.</caption><thead><tr><th></th><th>Truck</th><th>Distance</th><th>Deadline Exceed</th><th>Working Time Exceed</th><th>Runtime</th></tr></thead><tbody><tr><td>Center</td><td>3</td><td>516.35</td><td>316.99</td><td>00.00</td><td>3020.46</td></tr><tr><td>North</td><td>4</td><td>260.83</td><td>95.56</td><td>00.00</td><td>1915.77</td></tr><tr><td>East</td><td>4</td><td>406.26</td><td>693.62</td><td>00.00</td><td>1263.43</td></tr><tr><td>West</td><td>5</td><td>407.86</td><td>48.78</td><td>00.00</td><td>1085.05</td></tr><tr><td>South</td><td>5</td><td>571.12</td><td>161.81</td><td>00.00</td><td>32.59</td></tr><tr><td>Total</td><td>21</td><td>2162.42</td><td>1316.76</td><td>00.00</td><td>7317.3</td></tr></tbody></table>

<table><caption>**Table 4.** TS and CS results for the real case.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Truck</th><th colspan="2">Distance</th><th colspan="2">Deadline Exceed</th><th colspan="2">Working Time Exceed</th><th colspan="2">Runtime</th></tr><tr><th>TS</th><th>CS</th><th>TS</th><th>CS</th><th>TS</th><th>CS</th><th>TS</th><th>CS</th><th>TS</th><th>CS</th></tr></thead><tbody><tr><td>Center</td><td>4</td><td>4</td><td>452.96</td><td>417.74</td><td>00.00</td><td>00.00</td><td>00.00</td><td>00.00</td><td>54.80</td><td>17.40</td></tr><tr><td>North</td><td>4</td><td>4</td><td>262.14</td><td>240.92</td><td>00.00</td><td>00.00</td><td>00.00</td><td>00.00</td><td>42.85</td><td>16.25</td></tr><tr><td>East</td><td>5</td><td>5</td><td>356.12</td><td>344.86</td><td>00.00</td><td>00.00</td><td>00.00</td><td>00.00</td><td>25.07</td><td>18.26</td></tr><tr><td>West</td><td>5</td><td>5</td><td>411.05</td><td>390.04</td><td>00.00</td><td>00.00</td><td>00.00</td><td>00.00</td><td>72.81</td><td>21.54</td></tr><tr><td>South</td><td>5</td><td>5</td><td>673.66</td><td>618.46</td><td>00.00</td><td>00.00</td><td>00.00</td><td>00.00</td><td>20.01</td><td>15.96</td></tr><tr><td>Total</td><td>23</td><td>23</td><td>2155.93</td><td>2012.02</td><td>00.00</td><td>00.00</td><td>00.00</td><td>00.00</td><td>215.54</td><td>89.41</td></tr></tbody></table>

The results show that version 2 (i.e., when the feasibility of only the best candidate is tested) has reduced runtime; however, the average distances and amount of trucks are not necessarily reduced.

The above is explained because CS and TS test the feasibility of each candidate and select the best of them; in contrast, the best candidate selected for version 2 could have worse values after testing the feasibility. However, version 2 reduces the runtime by increasing the accuracy (i.e., increasing the number of MC simulations) and the number of candidates of each iteration. Meanwhile, CS is executed in a greater time than TS because each candidate is chosen from the best move; in contrast, in TS, the candidates are selected randomly, performing fewer calculations in each iteration.

### 3.2. Case Ice Cream Industry Results

The city of Santiago, Chile is in a valley of 15,000 km². The big customers of the Ice Cream Industry are 157 establishments with 43 customers with hard time windows and stochastic service divided into five zones: center, north, east, west and south. Table 3 shows the distance, the number of trucks, deadline exceed, working time exceed, and runtime obtained by GAMS. The large number of customers generates a low-quality solution using the GAMS method. GAMS programming was carried out by creating subgroups throughout the different areas of Santiago. The total post-resolution time exceeds 20 h. A daily excess of a 22 h deadline is incurred with the truck fleet, and it is possible not to exceed the normal working hours of the drivers.

The results obtained with the TS and CS algorithms for the ice cream industry problem show a (1) clear decrease in processing time, (2) reduction in the distance traveled, (3) do not exceed the deadline, and (4) do not they generate extra work. The results delivered by TS and CS increase the overall quality of service by reducing the deadlines (Table 4).

## 4. Conclusions

In this work, the VRPTW-ST model is used to solve a routing problem in Santiago de Chile, specifically in the distribution center of the ice cream area of the company Unilever. The TS and CS optimization methods were validated by comparison with the GAMS algorithm. The results validated the use of CS for VRPTW-ST, and the resolution capacity
of the CS algorithm for problems with stochastic time-out variables was established. The
CS algorithm obtained better results compared with TS.

This study demonstrates that CS and TS can solve stochastic variables using MC
simulations at each stage.

A greater number of MC simulations in each algorithm iteration increases the robustness of the results but increases computation time. The effect of using the MC simulation is not similar for both heuristics; CS is subjected to more candidates than TS in each interaction.

The best candidate is selected to improve the model, and then the MC simulation is
run. These modifications (TSv2 and CSv2) reduce the processing time without increasing
the total distance traveled.

There is a slight advantage in the CS2 vs. TS2, showing the efficiency of incorporating chaotic neurons when choosing candidates. Chaotic neurons generate greater diversification in the search space for possible solutions.

A model with probability restrictions concerning a model with deterministic variables (i) requires more time for resolution, (ii) reduces the maximum number of customers in the model, and (iii) increases the distances and vehicles required for resolution. The uncertainty of the waiting times generates an increase in the calculation times of the optimal solution; if the number of customers increases with stochastic waiting times, the complexity of the model increases. This model uncertainty must be controlled using different levels of services in F(x) to categorize each customer.

To validate the implementation of TS and CS methods, the results were compared with the authors' works [32,39]. The adjustment of the TS and CS parameters improved convergence and achieved better results compared with the GAMS algorithm. The improved performance in the distribution of refrigerated elements of the Unilever Company is achieved by the best search features of the CS method for VRPTW-ST models, (i) compliance with restrictions without decreasing service levels; (ii) the ability to plan routes in shorter times; (iii) shorter total routes, and (iv) reducing cost overruns by controlling overtime.

The experiments were performed comparing real models with a size limit of 20 cus-
tomers for the three optimization techniques. Using the GAMS technique with sizes greater
than 20 customers, the optimization model developed in the research cannot be solved
quickly enough for operational decision making. The applicability limits of the experiment
were restricted to 80 customers using the CS optimization techniques. The growth of
processing times and computer memory generates operating restrictions for the model
with a client size equal to or less than 100.

The high randomness of the waiting time generates a very costly restriction.
The decision-maker must select the product delivery policy with the confidence levels that
best suit their needs by modifying the constraint values of the proposed model.

The use of efficient mathematical models positively influences the correct daily decision making of a company. This type of optimization tool facilitates the work of administrators, reducing costs and above all, time, allowing adjustments to be quickly made before possible eventualities that huddle in the logistic chain.

For future work, we propose to consider (i) the randomness of travel times. This random time variable would allow us to adjust the model's results to the real conditions in a city with high traffic flow.

The model presented considers deterministic customer demand. Therefore, another extension of the model corresponds (ii) to the uncertain customer demand. Finally, (iii) the randomness presented in the future extended model requires a study on execution times. This study will allow us to establish the real computational cost of each random variable.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3
390/math9212750/s1, Table S1: Averages of the execution time, amount of trucks, and distance
obtained for the TS method with $\alpha = 0.5$ according to the confidence levels $\beta$ and $\gamma$, and percentage
of supermarkets, Table S2: Averages of the execution time, amount of trucks, and distance obtained
for the TSv2 method with $\alpha = 0.75$ according to the confidence levels $\beta$ and $\gamma$, and percentage of
supermarkets, Table S3: Averages of the execution time, amount of trucks, and distance obtained
for the CS method with $\alpha = 0.5$ according to the confidence levels $\beta$ and $\gamma$, and percentage of supermarkets, Table S4: Averages of the execution time, amount of trucks, and distance obtained for the CSv2 method with $\alpha = 0.5$ according to the confidence levels $\beta$ and $\gamma$, and percentage of supermarkets.

**Author Contributions:** Conceptualization, S.D., M.A. and G.F.; methodology, M.V., G.F. and S.D.; validation, M.A. and M.C.; formal analysis, M.C. and M.A.; Investigation, S.D., M.A. and M.V.; writing—original draft preparation, M.C. and M.V.; writing—review and editing, G.F., S.D. and M.V.; project administration, M.A., M.C. and G.F.; funding acquisition, S.D., M.C. and M.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been supported by DICYT (Scientific and Technological Research Bureau) of the University of Santiago of Chile (USACH) and Department of Industrial Engineering. This work was supported in part by Fondecyt (Chile) Grant No. 11200993 (M.V.).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** This research has been supported by DICYT (Scientific and Technological Research Bureau) of the University of Santiago of Chile (USACH) and Department of Industrial Engineering.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Notations**

<table><tr><td>TS</td><td>Tabu search</td></tr><tr><td>CS</td><td>Chaotic search</td></tr><tr><td>TW</td><td>Time windows</td></tr><tr><td>ST</td><td>Service times</td></tr><tr><td>MC</td><td>Monte Carlo</td></tr><tr><td>TSv2</td><td>Modified tabu search algorithm</td></tr><tr><td>CSv2</td><td>Modified chaotic search algorithm</td></tr><tr><td>VRP</td><td>Vehicle routing problems</td></tr><tr><td>SVRP</td><td>Stochastic vehicle routing problem</td></tr><tr><td>VRPD</td><td>Vehicle routing problem with deadline</td></tr><tr><td>VRPSD</td><td>Vehicle routing problem with stochastic demand</td></tr><tr><td>GAMS</td><td>General algebraic modeling systems</td></tr><tr><td>VRPTW</td><td>Vehicle routing problem with time windows</td></tr><tr><td>GRASP</td><td>Greedy randomized adaptive search procedure</td></tr><tr><td>RVRPTW</td><td>Robust vehicle routing problem with time windows</td></tr><tr><td>SDVRPTW</td><td>Single-depot vehicle-routing problem with time windows</td></tr><tr><td>VRPTW-ST</td><td>Vehicle routing problem with hard time windows and stochastic service time</td></tr><tr><td>SDVRP-MDA</td><td>Split-delivery vehicle routing problem with minimum delivery amounts</td></tr></table>

**References**

1. Theophilus, O.; Dulebenets, M.A.; Pasha, J.; Lau, Y.Y.; Fathollahi-Fard, A.M.; Mazaheri, A. Truck scheduling optimization at a cold-chain cross-docking terminal with product perishability considerations. *Comput. Ind. Eng.* **2021**, *156*, 107240. [CrossRef]

2. Yu, Y.; Xiao, T. Analysis of cold-chain service outsourcing modes in a fresh agri-product supply chain. *Transp. Res. Part E Logist. Transp. Rev.* **2021**, *148*, 102264. [CrossRef]

3. Zhang, S.; Chen, N.; She, N.; Li, K. Location optimization of a competitive distribution center for urban cold chain logistics in terms of low-carbon emissions. *Comput. Ind. Eng.* **2021**, *154*, 107120. [CrossRef]

4. He, B.; Yin, L. Prediction Modelling of Cold Chain Logistics Demand Based on Data Mining Algorithm. *Math. Probl. Eng.* **2021**, *3421478*. [CrossRef]

5. Hu, B.; Huang, B.; Liu, Z.; Guo, H.; Chen, Z.; Shi, L. Optimization Model of Carbon Footprint of Fresh Products in Cold Chain from the Energy Conservation and Emission Reduction Perspective. *Math. Probl. Eng.* **2021**, *2021*, 5559021. [CrossRef]
6. Li, X.; Zhou, K. Multi-objective cold chain logistic distribution center location based on carbon emission. *Environ. Sci. Pollut. Res.* **2021**, *28*, 32396–32404. [CrossRef]

7. Papadimitriou, C.H.; Steiglitz, K. *Combinatorial Optimization: Algorithms and Complexity*; Prentice-Hall: Hoboken, NJ, USA, 1982.

8. Demirel, E.; Demirel, N.; Gökçen, H. A mixed integer linear programming model to optimize reverse logistics activities of end-of-life vehicles in Turkey. *J. Clean. Prod.* **2016**, *112*, 2101–2113. [CrossRef]

9. Gupta, S.; Haq, A.; Ali, I.; Sarkar, B. Significance of multi-objective optimization in logistics problem for multi-product supply chain network under the intuitionistic fuzzy environment. *Complex Intell. Syst.* **2021**, *7*, 2119–2139. [CrossRef]

10. López-Ospina, H.; Cortés, C.E.; Pérez, J.; Peña, R.; Figueroa-García, J.C.; Urrutia-Mosquera, J. A maximum entropy optimization model for origin-destination trip matrix estimation with fuzzy entropic parameters. *Transp. A Transp. Sci.* **2021**, *1–39*. [CrossRef]

11. Gmira, M.; Gendreau, M.; Lodi, A.; Potvin, J.Y. Tabu search for the time-dependent vehicle routing problem with time windows on a road network. *Eur. J. Oper. Res.* **2021**, *288*, 129–140. [CrossRef]

12. Heidari, A.; Dong, Z.Y.; Zhang, D.; Siano, P.; Aghaei, J. Mixed-integer nonlinear programming formulation for distribution networks reliability optimization. *IEEE Trans. Ind. Inform.* **2018**, *14*, 1952–1961. [CrossRef]

13. Fuertes, G.; Navarrete, R.; Alfaro, M.; Millan, G.; Vargas, M.; Lagos, C. Adaptive equalization using artificial neural networks for a visible light communication system. In Proceedings of the IEEE International Conference on Automation, Concepcion, Chile, 17–19 October 2018; pp. 1–6.

14. Kumar, R.; Chandrawat, R.K.; Sarkar, B.; Joshi, V.; Majumder, A. An advanced optimization technique for smart production using α-cut based quadrilateral fuzzy number. *Int. J. Fuzzy Syst.* **2021**, *23*, 107–127. [CrossRef]

15. Garai, A.; Chowdhury, S.; Sarkar, B.; Roy, T.K. Cost-effective subsidy policy for growers and biofuels-plants in closed-loop supply chain of herbs and herbal medicines: An interactive bi-objective optimization in T-environment. *Appl. Soft Comput.* **2021**, *100*, 106949. [CrossRef]

16. Sarkar, B.; Mondal, S.P.; Hur, S.; Ahmadian, A.; Salahshour, S.; Guchhait, R.; Iqbal, M.W. An optimization technique for national income determination model with stability analysis of differential equation in discrete and continuous process under the uncertain environment. *RAIRO Oper. Res.* **2019**, *53*, 1649–1674. [CrossRef]

17. Kłosowski, G.; Kozłowski, E.; Gola, A. Integer linear programming in optimization of waste after cutting in the furniture manufacturing. In *Intelligent Systems in Production Engineering and Maintenance*; Springer International Publishing: Wroclaw, Poland, 2018; pp. 260–270.

18. Güiza, J.; Luque, R.; Murillo, J.; Romero, R.; Barrera, D.; López-Ospina, H. Integrating pricing and coordinated inventory decisions between one warehouse and multiple retailers. *J. Ind. Prod. Eng.* **2021**, *38*, 1–11. [CrossRef]

19. Kang, C.W.; Imran, M.; Omair, M.; Ahmed, W.; Ullah, M.; Sarkar, B. Stochastic-petri net modeling and optimization for outdoor patients in building sustainable healthcare system considering staff absenteeism. *Mathematics* **2019**, *7*, 499. [CrossRef]

20. Sabattin, J.; Fuertes, G.; Alfaro, M.; Quezada, L.; Vargas, M. Optimization of large electric power distribution using a parallel genetic algorithm with dandelion strategy. *Turk. J. Electr. Eng. Comput. Sci.* **2018**, *26*, 1–13. [CrossRef]

21. Samanta, S.; Dubey, V.K.; Sarkar, B. Measure of influences in social networks. *Appl. Soft Comput.* **2021**, *99*, 106858. [CrossRef]

22. Toth, P.; Vigo, D. *Vehicle Routing: Problems, Methods, and Applications*; SIAM: Philadelphia, PA, USA, 2014. ISBN 1611973589.

23. Gendreau, M.; Potvin, J.-Y. *Handbook of Metaheuristics*; Springer: Melbourne, Australia, 2010. ISBN 9781441916655.

24. Linfati, R.; Escobar, J.W. Reoptimization heuristic for the capacitated vehicle routing problem. *J. Adv. Transp.* **2018**, *2018*, 8. [CrossRef]

25. Feng, Y.; Zhang, R.-Q.; Jia, G. Vehicle routing problems with fuel consumption and stochastic travel speeds. *Math. Probl. Eng.* **2017**, *2017*, 16. [CrossRef]

26. Zeng, Z.; Xu, W.; Xu, Z.; Shao, W. A hybrid GRASP+VND heuristic for the two-echelon vehicle routing problem arising in city logistics. *Math. Probl. Eng.* **2014**, *2014*, 517467. [CrossRef]

27. Kaiwartya, O.; Kumar, S.; Lobiyal, D.K.; Tiwari, P.K.; Abdullah, A.H.; Hassan, A.N. Multiobjective dynamic vehicle routing problem and time seed based solution using particle swarm optimization. *J. Sens.* **2015**, *2015*, 14. [CrossRef]

28. Wang, Q.; Ji, Q.; Chiu, C.-H. Optimal routing for heterogeneous fixed fleets of multicompartment vehicles. *Math. Probl. Eng.* **2014**, *2014*, 847630. [CrossRef]

29. Zhang, X.; Zhong, S.; Liu, Y.; Wang, X. A framing link based tabu search algorithm for large-scale multidepot vehicle routing problems. *Math. Probl. Eng.* **2014**, *2014*, 13. [CrossRef]

30. Yue, Y.; Zhang, T.; Yue, Q. Improved fractal space filling curves hybrid optimization algorithm for vehicle routing problem. *Comput. Intell. Neurosci.* **2015**, *2015*, 9. [CrossRef]

31. Xu, H.; Pu, P.; Duan, F. A hybrid ant colony optimization for dynamic multidepot vehicle routing problem. *Discret. Dyn. Nat. Soc.* **2018**, *2018*, 3624728. [CrossRef]

32. Tao, N.-R.; Jiang, Z.-H.; Liu, J.-F.; Xia, B.-X.; Li, B.-H. A metaheuristic algorithm to transporter scheduling for assembly blocks in a shipyard considering precedence and cooperating constraints. *Discret. Dyn. Nat. Soc.* **2019**, *2019*, 14. [CrossRef]

33. Urzúa-Morales, J.G.; Sepulveda-Rojas, J.P.; Alfaro, M.; Fuertes, G.; Ternero, R.; Vargas, M. Logistic modeling of the last mile: Case study Santiago, Chile. *Sustainability* **2020**, *12*, 648. [CrossRef]

34. Hsu, C.I.; Hung, S.F.; Li, H.C. Vehicle routing problem with time-windows for perishable food delivery. *J. Food Eng.* **2007**, *80*, 465–475. [CrossRef]
35. Errico, F.; Desaulniers, G.; Gendreau, M.; Rei, W.; Rousseau, L.M. The vehicle routing problem with hard time windows and stochastic service times. *EURO J. Transp. Logist.* **2018**, *7*, 223–251. [CrossRef]

36. Dantzig, G.B.; Ramser, J.H. The truck dispatching problem. *Manag. Sci.* **1959**, *6*, 80–91. [CrossRef]

37. Taş, D.; Dellaert, N.; van Woensel, T.; de Kok, T. Vehicle routing problem with stochastic travel times including soft time windows and service costs. *Comput. Oper. Res.* **2013**, *40*, 214–224. [CrossRef]

38. Ehmke, J.F.; Campbell, A.M.; Urban, T.L. Ensuring service levels in routing problems with time windows and stochastic travel times. *Eur. J. Oper. Res.* **2015**, *240*, 539–550. [CrossRef]

39. Xu, H.; Pu, P.; Duan, F. Dynamic vehicle routing problems with enhanced ant colony optimization. *Discret. Dyn. Nat. Soc.* **2018**, *2018*, 1295485. [CrossRef]

40. Yu, N.-K.; Hu, R.; Qian, B.; Wang, L. An improved lagrangian relaxation algorithm for solving the lower bound of production logistics. In *Proceedings of the International Conference Intelligent Computing Theories and Application, Shenzhen, China, 12–15 August 2021*; Springer: Shenzhen, China, 2021; pp. 652–662.

41. Laporte, G.; Louveaux, F.; Mercure, H. The vehicle routing problem with stochastic travel times. *Transp. Sci.* **1992**, *26*, 161–170. [CrossRef]

42. Han, S.; Zhao, L.; Chen, K.; Luo, Z.W.; Mishra, D. Appointment scheduling and routing optimization of attended home delivery system with random customer behavior. *Eur. J. Oper. Res.* **2017**, *262*, 966–980. [CrossRef]

43. Zhang, D.; Li, D.; Sun, H.; Hou, L. A vehicle routing problem with distribution uncertainty in deadlines. *Eur. J. Oper. Res.* **2021**, *292*, 311–326. [CrossRef]

44. Shukla, N.; Choudhary, A.K.; Prakash, P.K.S.; Fernandes, K.J.; Tiwari, M.K. Algorithm portfolios for logistics optimization considering stochastic demands and mobility allowance. *Int. J. Prod. Econ.* **2013**, *141*, 146–166. [CrossRef]

45. Han, A.F.W.; Chu, Y.C. A multi-start heuristic approach for the split-delivery vehicle routing problem with minimum delivery amounts. *Transp. Res. Part E Logist. Transp. Rev.* **2016**, *88*, 11–31. [CrossRef]

46. Braaten, S.; Gjønnes, O.; Hvattum, L.M.; Tirado, G. Heuristics for the robust vehicle routing problem with time windows. *Expert Syst. Appl.* **2017**, *77*, 136–147. [CrossRef]

47. Marinakis, Y.; Marinaki, M.; Migdalas, A. A hybrid clonal selection algorithm for the location routing problem with stochastic demands. *Ann. Math. Artif. Intell.* **2016**, *76*, 121–142. [CrossRef]

48. Binart, S.; Dejax, P.; Gendreau, M.; Semet, F. A 2-stage method for a field service routing problem with stochastic travel and service times. *Comput. Oper. Res.* **2016**, *65*, 64–75. [CrossRef]

49. Li, X.; Tian, P.; Leung, S.C.H. Vehicle routing problems with time windows and stochastic travel and service times: Models and algorithms. *Int. J. Prod. Econ.* **2010**, *125*, 137–145. [CrossRef]

50. Ikeguchi, T.; Hasegawa, M.; Kimura, T.; Matsuura, T.; Aihara, K. Theory and applications of chaotic optimization methods. In *Innovative Computing Methods and Their Applications to Engineering Problems*; Springer: Berlin/Heidelberg, Germany, 2011; pp. 131–161.

51. Berhan, E.; Beshah, B.; Kitaw, D.; Abraham, A. Stochastic vehicle routing problem: A literature survey. *J. Inf. Knowl. Manag.* **2014**, *13*, 1450022. [CrossRef]

52. Sawik, B. Weighted-sum approach for bi-objective optimization of fleet size with environmental aspects. In *Applications of Management Science*; Emerald Publishing Limited: Bingley, UK, 2018; pp. 101–116.

53. Kato, H.; Kimura, T.; Ikeguchi, T. Self-organized neural network structure depending on the STDP learning rules. In *Applications of Nonlinear Dynamics*; Springer: Berlin/Heidelberg, Germany, 2009; pp. 413–416.