# Metaheuristic Approaches for Solving Truck and Trailer Routing Problems with Stochastic Demands: A Case Study in Dairy Industry

Authors: Seyedmehdi Mirmohammadsadeghi; Shamsuddin Ahmed
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Manufacturers and service providers often encounter stochastic demand scenarios. Researchers have, thus far, considered the deterministic truck and trailer routing problem (TTRP) that cannot address ubiquitous demand uncertainties and/or other complexities. The purpose of this study is to model the TTRP with stochastic demand (TTRPSD) constraints to bring the TTRP model closer to a reality. The model is solved in a reasonable timeframe using data from a large dairy service by administering the multipoint simulated annealing (M-SA), memetic algorithm (MA), and tabu search (TS). A sizeable number of customers whose demands follow the Poisson probability distribution are considered to model and solve the problem. To make the solutions relevant, first, 21 special TTRPSD benchmark instances are modified for this case and then these benchmarks are used in order to increase the validity and efficiency of the aforementioned algorithms and to show the consistency of the results. Also, the solutions have been tested using sensitivity analysis to understand the effects of the parameters and to make a comparison between the best results obtained by three algorithms and sensitivity analysis. Since the differences between the results are insignificant, the algorithms are found to be appropriate and relevant for solving real-world TTRPSD problem.


---

## Body

## 1. Introduction

These days, complex customer demands are required to be satisfied by many companies. Therefore, a large number of companies are trying to achieve a high level of reliability, flexibility, and agility for different demands. As a result, supply chain management (SCM) has become a thought-provoking subject for various companies, seeking for a way out of efficiently improving their customer satisfaction. In a way, according to the position and role, supply chain is categorized into three classes; the outbound, intracompany, and inbound supply chain. As the network of supplies begins at the inbound supply chain, the role of this group is transporting the semifinished products or the raw materials to the site of manufacturing. The main concern of the intracompany supply chain, as the intermediary part, is with the flow of material in the site of manufacturing. Finally, the outbound supply chain is concerned with the delivery of final products to the customers [1]. The inventory allocation and transportation are strongly considered in the outbound supply chain for minimizing the cost and satisfying the customers’ demands. One significant part of the supply chain management is providing the services or/and goods from a supply point to different destinations, which are geographically distributed with significant implications of economics. Aside from the cost of purchasing the goods, on the average and compared to the other relative activities, a higher percentage of the costs of logistics are absorbed by transportation. Therefore, efficiency improvement through the maximum usage of the necessities of transportation and decreasing the costs of transportation along with the improvement of services for customers are the frequent and significant decision analysis problems [2].Customers, warehouses, manufacturers, and suppliers are the main elements of a supply chain (SC), carrying the goods from the upstream to the downstream links of the chain. In a supply chain, there are four main business functions to be performed: purchasing, manufacturing, inventory, and distribution. The function of distribution includes two activities: the shipment of finished products from the companies to the locations of demand, and transportation of parts or raw materials from the suppliers to the companies [3].In order to manage a supply chain, a large number of business processes need to be carried out and many decisions are required to be made. Particular design versions of these general supply systems and inventory planning problems have been studied for a long time. It is pretty obvious that these main supply chain problems are greatly related. As the time goes by, more companies are awakened about their supply chain performance and how important it is that they improve this performance. They also have become aware of the competitive advantage of distribution operations, inventory, integration, and coordination of supplies. One of the main problems in supply chain management and logistics is the routing of a series of vehicles, which are assigned to transport goods from a warehouse to the customers or/and retailers. Since goods are hardly ever produced and consumed in one particular place, transportation is considered as a significant factor in the supply chain management.In this study, the supply chain is considered in terms of transportation and distribution. Due to the complexity, a large number of realistic solutions are disregarded. Any company in the world currently faces with a number of challenges in serving their customers. Transportation cost is considered to be the largest logistics expense for a vast number of firms and companies. Transportation is the area where costs should be diminished. This is a very bearing question, how servicing and manufacturing companies could successfully diminish transportation expenses without overshadowing customers’ satisfaction [4].It is a widely accepted principle that firms aiming to serve the customers scattered in a vast area should possess a servicing plan if they do not want to waste time and money. One of the best approaches to work out the arising problematic issues is to apply a special and unique method under the title of vehicle routing problem (VRP). Also, the truck and trailer routing problem (TTRP) is one of the widely studied and most important combinatorial optimization problems in VRP domain and because of its natural complications and efficacies in a large number of real world and supply chain management applications, it is finding its place in the transportation logistic.This paper discusses a real case study under TTRP in which demands are stochastic (TTRPSD) in nature. This work is an advancement of the well-known deterministic TTRP. TTRP is a variant of the vehicle routing problem (VRP). It is related to transporting manufacturing goods within a plant or between factory floors and delivering products to intended markets or customers. Indeed, VRP has been known as one of the most studied combinatorial optimization problems in the past few decades. It covers certain areas in practice and considers complexities to a reasonable extent such as stochastic VRP, multi depot VRP and TTRP [5–8]. Dantzig and Ramser [8] defined VRP as a generalized solution based on Travelling Salesman Problem (TSP) in 1959. After that many researchers have extended this concept to make it useful in diverse areas. During the last two decades, some constraints such as time window, travel and service time were added to the VRP. Due to some other practical issues, such as narrow roads and bridges in village or government restrictions, maneuvering a complete vehicle sometime appears to be difficult - the traditional VRP approach is found inadequate and these issues are considered in TTRP model. In general, TTRP is more extensive than VRP and can cover more real life aspects since some limitations in VRP as mentioned above can be considered in TTRP.In TTRP, the customers may be serviced either by a single truck or complete vehicle (truck with a trailer). This feature is usually ignored in VRP. However, because of some obstacles that appear in real-life situations, such as road conditions, market locations, government regulations or limited space to manoeuver at the customer’s site, only a single truck may be needed to serve a few destinations and/or customers. These constraints are obvious in many practical situations [9–11]. Several researchers have so far contributed in this area. One instance is that of Gerdessen’s work on VRP with a trailer in 1996 [12]. He demonstrated two real applications pertaining to TTRP. In one case it was the distribution of dairy products in cities where the distributor faced heavy traffic. Due to the fact that maneuvering a complete vehicle (a truck with a trailer) was difficult in some areas, and that some streets had traffic restrictions not allowing trailers to enter, the trailers were often parked in parking spots from which customers were serviced. The second case was the distribution of animal feed components to farmers. Because most villages have narrow roads or small bridges, different kinds of vehicles were required to deliver the feed to farmers, one of which was called a double bottom, consisting of a truck and a trailer. While the truck made deliveries using subtours (some parts of the tour) the trailer parked in a parking area. Gerdessen showed the necessity of considering TTRP by demonstrating real applications.In practical situations, a dispatcher may not know the exact demands in advance. Therefore, the company may face the problem of delivering products to customers with random demands. Consequently, unexpected extra costs might be imposed on the company. These issues can be considered in vehicle routing problems with stochastic demands (VRPSD). Moreover, when facing the limitations and restrictions mentioned above, VRPSD cannot cover these issues and needs to consider TTRP with stochastic demands (TTRPSD). The literature survey revealed no TTRP solution with stochastic parameters. Few articles were published on TTRP with deterministic parameters. However, papers published on stochastic VRP are simply large in number. These concepts need to be considered together in order to formulate TTRPSD. Therefore, this paper classifies the relevant papers into two groups—standard TTRP papers and VRP with stochastic demands. To solve TTRPSD, it appears that its solution is computationally more difficult than VRPSD. Indeed, VRP itself is one of the most difficult combinatorial optimization modelling problems. It is generally framed and solved by heuristics approaches [13–15]. To formulate and solve TTRPSD, one can make an effort to reduce it to VRPSD. As VRPSD is also a complex type of combinatorial optimization problem, it can be solved by heuristics approaches [7, 10, 16]. The purpose of TTRPSD is to design a set of routes to cover all customers and optimize the total costs or distance that will satisfy the relevant constraints.Standard TTRP was first proposed by Chao [7]. Scheuerer applied 0-1 integer programming formulation to solve TTRP [16]. Chao and Scheuerer used a two-phase approach to solve TTRP [7, 16]. They used heuristics to construct the initial TTRP solution in the first phase. In the second phase, the Tabu search algorithm was used to improve the initial solution. Chao followed Fisher and Jaikumar’s construction to solve VRP [17]. Scheuerer [16] used Chao’s model [7] and improved it by using two-construction heuristics, T-Cluster and T-Sweep, and applied a new Tabu search improvement algorithm to solve TTRP. Lin et al. [9] introduced simulated annealing to solve TTRP and obtained 17 best solutions to the 21 benchmarked TTRPs. Then they applied time window constraints to the TTRP solution for the first time to bring the model closer to reality [18]. Villegas et al. [19] considered a single TTRP with satellite depots (STTRPSD). Variable neighbourhood descent (VND) and greedy randomized adaptive search procedures (GRASP) were proposed by them for solving TTRP. In addition, they applied the GRASP/VND algorithm for multidepot VRPs and improved the previous analysis. Villegas et al. [20] improved this solution by considering a hybrid algorithm based on the GRASP/VND and a path relinking (PR) algorithm and proved that the performance of this hybrid algorithm exceeds that of GRASP/VND. Finally, Villegas et al. [11] proposed a new hybrid algorithm by considering GRASP and an iterated local search (ILS) and found a new solution for benchmarking, which was considered by Lin et al. for the first time [9]. Derigs et al. proposed TTRP without transferring load between truck and trailer for the first time [10]. A hybrid algorithm was applied to solve the TTRP by considering the large neighbourhood search (LNS) and local search (LS). In addition, time window constraints were also considered for each customer to bring the model closer to reality. Tillman introduced stochastic demand in capacitated VRP (CVRP) for the first time [21, 22]. In his study, the multidepot VRP was considered in which the demands are stochastic with Poisson distribution. Bertsimas demonstrated critical evidence for the stochastic VRP and compared reoptimization anda priori strategies [23]. Gendreau et al. [24] presented a Tabu search algorithm for solving VRP with stochastic demands and customers’ demands following a known distribution and customers were presented with a probability. The integer L-shaped method for CVRP with normal and Poisson demands was improved by Laporte et al. [25]. Lei et al. [26] extended CVRP with stochastic demand (CVRPSD), and time window constraints were added to the vertices. They solved the vehicle routing problem with stochastic demands and time windows (VRPSDTW) by considering discrete and continuous distributions for stochastic demands.In this study, the multipoint simulated anneling (M-SA), memetic algorithm (MA) and tabu search (TS) algorithms are applied in solving the TTRPSD. Multipoint SA is a local search heuristic method, which uses a few initial solutions instead of one solution and generates better solutions in the process in order to find the best solutions. The MA is a hybrid algorithm which combines the population and LS procedures to find the best solutions. Tabu search is an iterative local search heuristics. A tabu mechanism is put in place to prevent the process from cycling over a sequence of solutions. This paper demonstrates these algorithms that are efficient for solving TTRPSD. The new real case study has been considered for this purpose. This problem has been solved by these algorithms and the results have been compared with each other. In addition, sensitivity analysis is conducted to validate the results.This paper is organized as follows. Section2 provides the formulation of the TTRPSD. The detailed implementation of metaheuristic algorithms for solving the TTRPSD is described in Section 3. The computational study is reported in Section 4. The paper is concluded in Section 5. In addition, some further researches are recommended in this section.


## 2. Formulation of Truck and Trailer Routing Problem with Stochastic Demand

TTRPSD is defined as an undirected graphG
=
(
V
,
E
), where the set of vertices is V
=
v
0
,
v
1
,
v
2
,
…
,
v
n and E
=
{
(
v
i
,
v
j
)
:
v
i
,
v
j
∈
V
,
i
<
j
} is the set of edges. The central depot is represented by “v
0” and the other vertices in V
∖
{
v
0
} correspond to customers. Each vertex v
i is associated with a stochastic and nonnegative demand ξ
i. They can be split and are unknown until the vehicle arrives at the vertex. A customer type t
i is available for all customers, where t
i
=
1 shows that customer v
i is a truck customer (TC) and can be serviced only by a single truck. If t
i
=
0, a customer i is a vehicle customer (VC) and it can be serviced by a single truck or a complete vehicle (truck pulling a trailer). C
=
(
c
(
v
i
,
v
j
)
) is a symmetric travel cost, which is defined on E. It is assumed that all vehicles are the same and have a unit speed, so the travel cost is equal to the Euclidean distance between v
i and v
j. Say a fleet of m
k and m
r, trucks and trailers, respectively, is available. However, some trucks and trailers may not be used at any time in the TTRPSD solution. Without losing generality, it is assumed that m
k
≥
m
r, as in Chao [7], Scheuerer [16] and Lin et al. [9]. All trucks have the same capacity Q
k, and all trailers also have the same capacity Q
r, where Q
k and Q
r are different.Three types of route are available in TTRPSD as follows:(
1
) a pure travel route (PTR) which can be travelled by only a single truck; (
2
) a pure vehicle route without any subtours (PVR), which can only be travelled by a complete vehicle; (
3
) complete vehicle route (CVR) which consists of a main tour and at least one subtour. A subtour starts and finished at the same VC (v
r) (the trailer is parked in a parking place which is called the root) and it can be travelled only by a single truck; however, it should be serviced by a complete vehicle in the main tour. It is assumed that shifting demand loads between the truck and trailer is possible at the parking places. In addition, it should be mentioned that there may be a lack of capacity in serving the customers since the customers’ demands are stochastic. Therefore, the vehicle must return to the depot or to the parking place (while the truck is delivering on the subtours) and refill to capacity to serve the remaining customers. This is the so-called failure in the route.The TTRPSD can be cast as stochastic programming. Two main solution concepts for solving the aforesaid types of TTRPSD can be discerned from stochastic programming [26]. The first is known as chance constrained programming (CCP). In CCP, the problem can be solved by imposing a constraint ensuring that the probability of route failure is bounded by some parameters, such as time limitation and service time [26, 27]. This concept attempts to convert stochastic parameters to equivalent deterministic values. For instance, the TTRPSD can be converted to an equivalent deterministic program. Stewart and Golden [28] and Laporte et al. [29] demonstrated this transformation by considering the statistical relationships among the parameters. The second concept is stochastic programming with recourse (SPR). Two main solution strategies are available under SPR. The first is known asa priori optimization [23, 30–32] and the second is reoptimization [33, 34]. In ana priori optimization solution, the set of tours and subtours is constructed in the first stage. Recourse actions considering random variables are then revealed. The common recourse policy in TTRPSD is as follows: if the cumulative demand exceeds the vehicle capacity on the main tour (it means that the cumulative demand exceeds the summation of capacities of truck and trailer on the main tour up to a vertex), the vehicle returns to the depot to unload, fills the capacity (truck and trailer) and then comes back to the last visited vertex. However, if the vehicle capacity is reached exactly for any vertex on the main tour (it means that the cumulative demand is equal to the summation of capacities of truck and trailer on the main tour up to a vertex), the vehicle should return to the depot, fills the capacity and then start servicing from the next vertex along its route. But if the cumulative demand exceeds the truck capacity on the subtour, the vehicle should return to its root to upload, use the remaining trailer capacity and return to the last visited vertex on the subtour. In addition, if the truck capacity exactly reached the customer demand for any vertex on the subtour (it means that the cumulative demand on the subtour is equal to the truck capacity up to a vertex), the truck returns to the root, using the remaining trailer capacity to fill the truck and comes back to the next vertex on the subtour.The TTRPSD consists of designing the first-stage routes, including the truck route, vehicle route and complete route, satisfying all constraints and if any failure occurs, a recourse action is applied. The purpose is to minimize the sum of the total distance of the first-stage routes and the expected recourse costs.


### 2.1. The Expected Cost of a Solution

Since all trucks and trailers in their respective groups are the same, each having its own features in terms of capacity and associated unit speed, the cost is associated with the distance travelled. Initially, the expected recourse cost (distance) is computed under some assumptions. Then the total cost is computed and the recourse costs will be computed based on these assumptions.It should be mentioned that all demands are independent random variables with known distributions. Also each demandξ
i is a nonnegative random variable that never exceeds the truck capacity. In addition, a maximum of one return to the root is possible on any subtour (ST) and a maximum of one return to the depot can be allowed on any main tour (in a complete route), truck or vehicle route. Because more than one failure is impossible on any kind of route r
k
(1)
P
∑
v
i
∈
r
k
ξ
i
<
2
Q
k
+
Q
r
≅
1
I
f
r
k
∈
P
V
R
,
C
V
R
,
P
∑
v
i
∈
r
k
ξ
i
<
2
Q
k
≅
1
I
f
r
k
∈
P
T
R
,
S
T
.This failure means that the demand of a customer cannot be satisfied while serving the customers since the vehicle does not have enough capacity to serve the respective customer and has to come back to the depot (or parking place) and fill the capacity and return to the respective customer to serve it completely. This is to ensure that this failure may occur a maximum once. Therefore, the cumulative demand of the customers must be less than twice the route capacity in a worst case.


### 2.2. Mathematical Estimation of the Total Expected Cost

It should be considered that TTRPSD has three different types of route. Each route is planned in the first stage of the solutionR
=
(
r
1
,
r
2
,
…
,
r
m
), where r
k
=
(
v
0
k
=
v
0
,
v
1
k
,
…
,
v
n
k
k
,
v
n
k
+
1
k
=
v
0
). The final objective function T
(
R
) is the sum of two costs: Φ
R and R
(
R
), where Φ
(
R
) is the total distance of the first-stage routes, and R
(
R
) is the expected recourse cost. So, the objective function is,(2)
T
R
=
Φ
R
+
R
R
.The computation of Φ
R is not difficult. To compute Φ
R, the cost of each route should be found and the costs (total distances) of all routes, which are planned in the first stage, should be computed. The estimations of θ
(
R
) is complicated and it is shown in the next section. First, the probability of failure should be computed. ξ
i
k and X
i
k are the demand from customer v
i and the cumulative demand up to customer v
i in a route r
k.The probability of the cumulative demand up to the customerv
i
k, on a route r
k, can be calculated by the following recursive equation [24], where P
ξ
0
k
(
t
) is a boundary condition for this equation and P
X
i
k
t
=
P
{
X
i
k
=
t
}
(3)
P
X
i
k
t
=
∑
d
=
0
t
P
X
i
-
1
k
t
-
d
P
ξ
i
d
.For instance, consider P
X
0
k
t
=
P
ξ
0
k
(
t
), P
X
4
k
3
=
P
X
3
k
(
3
)
P
ξ
i
(
0
)
+
P
X
3
k
2
P
ξ
i
1
+
P
X
3
k
1
P
ξ
i
2
+
P
X
3
k
(
0
)
P
ξ
i
(
3
). This means that the cumulative demand up to customer 4 can be three if the cumulative demand up to customer 3 is three and the demand of customer 4 is zero, or the cumulative demand up to customer 3 is two and the demand of customer 4 is one, or the cumulative demand up to customer 3 is one and the demand of customer 4 is two, or the cumulative demand up to customer 3 is zero and the demand of customer 4 is three. Therefore, the probability of these cases should be summed up to calculate the total cumulative demand up to customer 4.The failure occurred at vertexv
i
k on a route r
k as follows. If failure occurred in a pure truck route or ST, then X
i
-
1
k
<
Q
k and X
i
k
≥
Q
k. However, X
i
-
1
k
<
Q
k
+
Q
r and X
i
k
≥
Q
k
+
Q
r, if failure occurred in a pure vehicle route or a main tour. Then depending on the type of route, the probability P
i
k can be computed. The probability of route failure P
i
k at customer v
i
k on a route r
k can be computed [26] as(4)
P
i
k
=
∑
d
=
0
Q
k
-
1
P
X
i
-
1
k
t
-
∑
d
=
0
Q
k
-
1
P
X
i
k
t
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
T
R
o
r
s
u
b
t
o
u
r
,
∑
d
=
0
Q
k
+
Q
r
-
1
P
X
i
-
1
k
t
-
∑
d
=
0
Q
k
+
Q
r
-
1
P
X
i
k
t
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
V
R
,
m
a
i
n
t
o
u
r
,
0
o
t
h
e
r
w
i
s
e
.For example, if the failure occurs on the main tour of CVR the probability P
i
k can be computed as(5)
P
i
k
=
P
X
i
-
1
k
<
Q
k
+
Q
r
,
X
i
k
≥
Q
k
+
Q
r
=
1
-
P
X
i
-
1
k
≥
Q
k
+
Q
r
-
P
X
i
k
<
Q
k
+
Q
r
=
P
X
i
-
1
k
<
Q
k
+
Q
r
-
P
X
i
k
<
Q
k
+
Q
r
=
∑
l
=
0
Q
k
+
Q
r
-
1
P
X
i
-
1
k
t
-
∑
l
=
0
Q
k
+
Q
r
-
1
P
X
i
k
t
.


### 2.3. The Expected Recourse Cost

Failure can happen in different situations. First, failure may occur in the PTR, PVR or main-route on the CVR; second failure may occur in an ST on the CVR. The recourse cost can be computed according to four failure types as follows.(1)
Type 1:H
i
k
=
c
(
v
i
k
,
v
r
k
)
+
c
v
r
k
,
v
i
+
1
k
-
c
(
v
i
k
,
v
i
+
1
k
), if X
i
k
+
ξ
i
k
=
Q
k, if failure occurred in the PTR or ST.
(2)
Type 2:G
i
k
=
2
c
(
v
i
k
,
v
r
k
), if X
i
k
+
ξ
i
k
>
Q
k, if failure occurred in the PTR or ST.
(3)
Type 3:F
i
k
=
c
(
v
i
k
,
v
0
k
)
+
c
v
0
k
,
v
i
+
1
k
-
c
(
v
i
k
,
v
i
+
1
k
), if X
i
k
+
ξ
i
k
=
Q
k
+
Q
r, if failure occurred in the PVR, main tour.
(4)
Type 4:E
i
k
=
2
c
(
v
i
k
,
v
0
k
), if X
i
k
+
ξ
i
k
>
Q
k
+
Q
r, if failure occurred in the PVR, main tour.For instance, if the failure is of the first type, it means that after serving customerv
i
k, the vehicle has to return to the root of the route if the vehicle is in ST, otherwise it must return to the depot, then continue serving customers from the next vertex v
i
+
1
k. If the second failure type occurs, after serving customer v
i
k, the single truck must return to the depot or to the root of its route and again proceed to vertex v
i
k. The third and fourth types of the failure are the same as the first and second one.The expected cost of recourse for router
k can be computed as(6)
E
θ
r
k
=
∑
i
=
2
n
k
A
i
k
H
i
k
+
B
i
k
G
i
k
+
C
i
k
F
i
k
+
D
i
k
E
i
k
,where A
i
k, B
i
k, C
i
k, D
i
k is the probability that first, second, third, fourth types of failure are occurring respectively. Considering (3) and (4), the two terms can be computed as(7)
A
i
k
=
∑
l
=
1
Q
k
P
ξ
i
k
l
P
X
i
-
1
k
Q
k
-
l
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
T
R
o
r
s
u
b
t
o
u
r
,
0
o
t
h
e
r
w
i
s
e
,
B
i
k
=
∑
l
=
1
Q
k
-
1
P
X
i
-
1
k
t
-
∑
l
=
1
Q
k
-
1
P
X
i
k
t
-
A
i
k
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
T
R
o
r
s
u
b
t
o
u
r
,
0
o
t
h
e
r
w
i
s
e
,
C
i
k
=
∑
l
=
1
Q
k
+
Q
r
P
ξ
i
k
l
P
X
i
-
1
k
Q
k
+
Q
r
-
l
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
V
R
,
m
a
i
n
t
o
u
r
,
0
o
t
h
e
r
w
i
s
e
,
D
i
k
=
∑
l
=
1
Q
k
+
Q
r
-
1
P
X
i
-
1
k
t
-
∑
l
=
1
Q
k
+
Q
r
-
1
P
X
i
k
t
-
C
i
k
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
V
R
,
m
a
i
n
t
o
u
r
,
0
o
t
h
e
r
w
i
s
e
.To recognize the failure types and their probabilities, it should be mentioned that the expected cost of a recourse is an extra travel cost. Therefore, four failure types and four recourse actions have been considered for this model. For instance, H
i
k is the cost of recourse action type 1 and should be multiplied by the relevant failure type which is A
i
k to compute the extra travel cost type 1. Also, G
i
k should be multiplied by B
i
k and so on.Finally, the total expected cost of recourse action in terms of distance can be computed as(8)
E
θ
R
=
∑
k
=
1
m
E
θ
r
ksubject to: (9)
∑
j
∈
V
∑
k
∈
K
X
i
j
k
+
Y
i
j
k
=
1
,
∀
i
∈
V
,
(10)
∑
j
∈
V
X
0
j
k
+
Y
0
j
k
=
1
,
∀
k
∈
K
,
(11)
∑
i
∈
V
X
i
0
k
+
Y
i
0
k
=
1
,
∀
k
∈
K
,
(12)
∑
i
∈
V
X
i
j
k
-
∑
i
∈
V
X
j
i
k
=
0
,
∀
j
∈
V
,
k
∈
K
,
(13)
∑
i
∈
V
Y
i
j
k
-
∑
i
∈
V
Y
j
i
k
=
0
,
∀
j
∈
V
,
k
∈
K
,
(14)
P
∑
i
∈
v
ξ
i
∑
j
∈
v
X
i
j
k
<
2
Q
k
+
Q
r
≅
1
I
f
v
i
∈
P
V
R
o
r
C
V
R
(15)
P
∑
i
∈
v
ξ
i
∑
j
∈
v
Y
i
j
k
<
2
Q
k
≅
1
I
f
r
k
∈
P
T
R
o
r
S
T
(16)
X
i
j
k
∈
0,1
,
Y
i
j
k
∈
0,1
,
∀
i
,
j
∈
V
,
k
∈
K
,where, r
k is the kth route and m is the number of routes including PTR, PVR and CVR. X
i
j
k is equal to 1 if customers i and j (edge (
i
,
j
)
∈
E) are serviced by a complete vehicle (the kth truck with a trailer), and 0 otherwise. Y
i
j
k is equal to 1 if edge (
i
,
j
)
∈
E is used by the kth truck (without trailer). Equation (9) gives information about each customer that must be serviced by a single truck or complete vehicle. Equation (10) is concerned with starting each route from the depot and going to exactly one customer using vehicle k. Equation (11) is similar to (10); however, it shows the route terminates by leaving one customer. Equations (12) and (13) indicate the flow on the route to be followed by vehicle k. Equations (14) and (15) are the capacity constraints for the vehicle route and the truck route, respectively, and these ensure the feasibility of servicing customers, and q
i is the expected value of the stochastic customers’ demands.


## 2.1. The Expected Cost of a Solution

Since all trucks and trailers in their respective groups are the same, each having its own features in terms of capacity and associated unit speed, the cost is associated with the distance travelled. Initially, the expected recourse cost (distance) is computed under some assumptions. Then the total cost is computed and the recourse costs will be computed based on these assumptions.It should be mentioned that all demands are independent random variables with known distributions. Also each demandξ
i is a nonnegative random variable that never exceeds the truck capacity. In addition, a maximum of one return to the root is possible on any subtour (ST) and a maximum of one return to the depot can be allowed on any main tour (in a complete route), truck or vehicle route. Because more than one failure is impossible on any kind of route r
k
(1)
P
∑
v
i
∈
r
k
ξ
i
<
2
Q
k
+
Q
r
≅
1
I
f
r
k
∈
P
V
R
,
C
V
R
,
P
∑
v
i
∈
r
k
ξ
i
<
2
Q
k
≅
1
I
f
r
k
∈
P
T
R
,
S
T
.This failure means that the demand of a customer cannot be satisfied while serving the customers since the vehicle does not have enough capacity to serve the respective customer and has to come back to the depot (or parking place) and fill the capacity and return to the respective customer to serve it completely. This is to ensure that this failure may occur a maximum once. Therefore, the cumulative demand of the customers must be less than twice the route capacity in a worst case.


## 2.2. Mathematical Estimation of the Total Expected Cost

It should be considered that TTRPSD has three different types of route. Each route is planned in the first stage of the solutionR
=
(
r
1
,
r
2
,
…
,
r
m
), where r
k
=
(
v
0
k
=
v
0
,
v
1
k
,
…
,
v
n
k
k
,
v
n
k
+
1
k
=
v
0
). The final objective function T
(
R
) is the sum of two costs: Φ
R and R
(
R
), where Φ
(
R
) is the total distance of the first-stage routes, and R
(
R
) is the expected recourse cost. So, the objective function is,(2)
T
R
=
Φ
R
+
R
R
.The computation of Φ
R is not difficult. To compute Φ
R, the cost of each route should be found and the costs (total distances) of all routes, which are planned in the first stage, should be computed. The estimations of θ
(
R
) is complicated and it is shown in the next section. First, the probability of failure should be computed. ξ
i
k and X
i
k are the demand from customer v
i and the cumulative demand up to customer v
i in a route r
k.The probability of the cumulative demand up to the customerv
i
k, on a route r
k, can be calculated by the following recursive equation [24], where P
ξ
0
k
(
t
) is a boundary condition for this equation and P
X
i
k
t
=
P
{
X
i
k
=
t
}
(3)
P
X
i
k
t
=
∑
d
=
0
t
P
X
i
-
1
k
t
-
d
P
ξ
i
d
.For instance, consider P
X
0
k
t
=
P
ξ
0
k
(
t
), P
X
4
k
3
=
P
X
3
k
(
3
)
P
ξ
i
(
0
)
+
P
X
3
k
2
P
ξ
i
1
+
P
X
3
k
1
P
ξ
i
2
+
P
X
3
k
(
0
)
P
ξ
i
(
3
). This means that the cumulative demand up to customer 4 can be three if the cumulative demand up to customer 3 is three and the demand of customer 4 is zero, or the cumulative demand up to customer 3 is two and the demand of customer 4 is one, or the cumulative demand up to customer 3 is one and the demand of customer 4 is two, or the cumulative demand up to customer 3 is zero and the demand of customer 4 is three. Therefore, the probability of these cases should be summed up to calculate the total cumulative demand up to customer 4.The failure occurred at vertexv
i
k on a route r
k as follows. If failure occurred in a pure truck route or ST, then X
i
-
1
k
<
Q
k and X
i
k
≥
Q
k. However, X
i
-
1
k
<
Q
k
+
Q
r and X
i
k
≥
Q
k
+
Q
r, if failure occurred in a pure vehicle route or a main tour. Then depending on the type of route, the probability P
i
k can be computed. The probability of route failure P
i
k at customer v
i
k on a route r
k can be computed [26] as(4)
P
i
k
=
∑
d
=
0
Q
k
-
1
P
X
i
-
1
k
t
-
∑
d
=
0
Q
k
-
1
P
X
i
k
t
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
T
R
o
r
s
u
b
t
o
u
r
,
∑
d
=
0
Q
k
+
Q
r
-
1
P
X
i
-
1
k
t
-
∑
d
=
0
Q
k
+
Q
r
-
1
P
X
i
k
t
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
V
R
,
m
a
i
n
t
o
u
r
,
0
o
t
h
e
r
w
i
s
e
.For example, if the failure occurs on the main tour of CVR the probability P
i
k can be computed as(5)
P
i
k
=
P
X
i
-
1
k
<
Q
k
+
Q
r
,
X
i
k
≥
Q
k
+
Q
r
=
1
-
P
X
i
-
1
k
≥
Q
k
+
Q
r
-
P
X
i
k
<
Q
k
+
Q
r
=
P
X
i
-
1
k
<
Q
k
+
Q
r
-
P
X
i
k
<
Q
k
+
Q
r
=
∑
l
=
0
Q
k
+
Q
r
-
1
P
X
i
-
1
k
t
-
∑
l
=
0
Q
k
+
Q
r
-
1
P
X
i
k
t
.


## 2.3. The Expected Recourse Cost

Failure can happen in different situations. First, failure may occur in the PTR, PVR or main-route on the CVR; second failure may occur in an ST on the CVR. The recourse cost can be computed according to four failure types as follows.(1)
Type 1:H
i
k
=
c
(
v
i
k
,
v
r
k
)
+
c
v
r
k
,
v
i
+
1
k
-
c
(
v
i
k
,
v
i
+
1
k
), if X
i
k
+
ξ
i
k
=
Q
k, if failure occurred in the PTR or ST.
(2)
Type 2:G
i
k
=
2
c
(
v
i
k
,
v
r
k
), if X
i
k
+
ξ
i
k
>
Q
k, if failure occurred in the PTR or ST.
(3)
Type 3:F
i
k
=
c
(
v
i
k
,
v
0
k
)
+
c
v
0
k
,
v
i
+
1
k
-
c
(
v
i
k
,
v
i
+
1
k
), if X
i
k
+
ξ
i
k
=
Q
k
+
Q
r, if failure occurred in the PVR, main tour.
(4)
Type 4:E
i
k
=
2
c
(
v
i
k
,
v
0
k
), if X
i
k
+
ξ
i
k
>
Q
k
+
Q
r, if failure occurred in the PVR, main tour.For instance, if the failure is of the first type, it means that after serving customerv
i
k, the vehicle has to return to the root of the route if the vehicle is in ST, otherwise it must return to the depot, then continue serving customers from the next vertex v
i
+
1
k. If the second failure type occurs, after serving customer v
i
k, the single truck must return to the depot or to the root of its route and again proceed to vertex v
i
k. The third and fourth types of the failure are the same as the first and second one.The expected cost of recourse for router
k can be computed as(6)
E
θ
r
k
=
∑
i
=
2
n
k
A
i
k
H
i
k
+
B
i
k
G
i
k
+
C
i
k
F
i
k
+
D
i
k
E
i
k
,where A
i
k, B
i
k, C
i
k, D
i
k is the probability that first, second, third, fourth types of failure are occurring respectively. Considering (3) and (4), the two terms can be computed as(7)
A
i
k
=
∑
l
=
1
Q
k
P
ξ
i
k
l
P
X
i
-
1
k
Q
k
-
l
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
T
R
o
r
s
u
b
t
o
u
r
,
0
o
t
h
e
r
w
i
s
e
,
B
i
k
=
∑
l
=
1
Q
k
-
1
P
X
i
-
1
k
t
-
∑
l
=
1
Q
k
-
1
P
X
i
k
t
-
A
i
k
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
T
R
o
r
s
u
b
t
o
u
r
,
0
o
t
h
e
r
w
i
s
e
,
C
i
k
=
∑
l
=
1
Q
k
+
Q
r
P
ξ
i
k
l
P
X
i
-
1
k
Q
k
+
Q
r
-
l
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
V
R
,
m
a
i
n
t
o
u
r
,
0
o
t
h
e
r
w
i
s
e
,
D
i
k
=
∑
l
=
1
Q
k
+
Q
r
-
1
P
X
i
-
1
k
t
-
∑
l
=
1
Q
k
+
Q
r
-
1
P
X
i
k
t
-
C
i
k
0
i
f
f
a
i
l
u
r
e
o
c
c
u
r
r
e
d
i
n
P
V
R
,
m
a
i
n
t
o
u
r
,
0
o
t
h
e
r
w
i
s
e
.To recognize the failure types and their probabilities, it should be mentioned that the expected cost of a recourse is an extra travel cost. Therefore, four failure types and four recourse actions have been considered for this model. For instance, H
i
k is the cost of recourse action type 1 and should be multiplied by the relevant failure type which is A
i
k to compute the extra travel cost type 1. Also, G
i
k should be multiplied by B
i
k and so on.Finally, the total expected cost of recourse action in terms of distance can be computed as(8)
E
θ
R
=
∑
k
=
1
m
E
θ
r
ksubject to: (9)
∑
j
∈
V
∑
k
∈
K
X
i
j
k
+
Y
i
j
k
=
1
,
∀
i
∈
V
,
(10)
∑
j
∈
V
X
0
j
k
+
Y
0
j
k
=
1
,
∀
k
∈
K
,
(11)
∑
i
∈
V
X
i
0
k
+
Y
i
0
k
=
1
,
∀
k
∈
K
,
(12)
∑
i
∈
V
X
i
j
k
-
∑
i
∈
V
X
j
i
k
=
0
,
∀
j
∈
V
,
k
∈
K
,
(13)
∑
i
∈
V
Y
i
j
k
-
∑
i
∈
V
Y
j
i
k
=
0
,
∀
j
∈
V
,
k
∈
K
,
(14)
P
∑
i
∈
v
ξ
i
∑
j
∈
v
X
i
j
k
<
2
Q
k
+
Q
r
≅
1
I
f
v
i
∈
P
V
R
o
r
C
V
R
(15)
P
∑
i
∈
v
ξ
i
∑
j
∈
v
Y
i
j
k
<
2
Q
k
≅
1
I
f
r
k
∈
P
T
R
o
r
S
T
(16)
X
i
j
k
∈
0,1
,
Y
i
j
k
∈
0,1
,
∀
i
,
j
∈
V
,
k
∈
K
,where, r
k is the kth route and m is the number of routes including PTR, PVR and CVR. X
i
j
k is equal to 1 if customers i and j (edge (
i
,
j
)
∈
E) are serviced by a complete vehicle (the kth truck with a trailer), and 0 otherwise. Y
i
j
k is equal to 1 if edge (
i
,
j
)
∈
E is used by the kth truck (without trailer). Equation (9) gives information about each customer that must be serviced by a single truck or complete vehicle. Equation (10) is concerned with starting each route from the depot and going to exactly one customer using vehicle k. Equation (11) is similar to (10); however, it shows the route terminates by leaving one customer. Equations (12) and (13) indicate the flow on the route to be followed by vehicle k. Equations (14) and (15) are the capacity constraints for the vehicle route and the truck route, respectively, and these ensure the feasibility of servicing customers, and q
i is the expected value of the stochastic customers’ demands.


## 3. Metaheuristic Algorithms to Solve TTRPSD

Metaheuristics use two principal methods to improve the solution from that of heuristic in terms of different performance criteria. These methods are known as local search method [26] and population search method [35]. When using the local search method, one should know that a thorough search of the solution space is implemented by moving the existing solution to another likely solution in its neighborhood at each step. Tabu search (TS) and simulated annealing (SA) are the most well-known algorithms in this area [36]. The population search includes upholding a pool of good parent solutions and then reassociating them to create new offsprings. Genetic algorithm (GA) is the main example in this principle. To solve TTRPSD, multipoint simulated annealing (M-SA), tabu search and memetic algorithms (MA) have been used. MAs fit into the category of the evolutionary algorithms (EAs) where LS procedures are used to distribute the search area and enhance the search to refine the individual. In fact, MA is a hybrid algorithm which combines the population and LS procedures to improve the initial solutions.In MAs, the populationP consists of a set of individuals generated randomly. Each individual is called a “chromosome,” which is simply a permutation of the set of n nodes (customers) 1,2
,
…
,
n and N dummy zeroes (artificial stores or the root of ST). This idea was initially proposed by Beasley [37] as part of a route-first cluster-second heuristic, and was then used by Prins [38]. The presentation of the solution for TTRPSD is explained more as follows. Initially, parents are selected based on tournament selection. Then, in each stage, offspring are randomly generated from the initial population. First, two chromosomes (parents) are randomly selected and two new chromosomes (offspring) are produced by crossover operation. The functional value is computed by the chromosomes. The new solutions are compared with the existing solutions. If a new solution is better than an existing solution, the algorithm replaces the existing solution with the new one. Then the new solution can be improved by mutations and LS procedures. This procedure is continued until the termination condition occurs.In M-SA, firstly, the number of predetermined initial solutions should be produced and it is known as the best solutions. In each iteration, the algorithm produces some new solutions from the neighborhood of the current solutions and finds the best one between them and chooses this as a new solution. If a new solution appears to be better than the current solution, this new solution is termed as a best solution and the procedure is continued. In this process, the number of predetermined iterations is applied in each temperature level to improve the possibility of a set of appropriate solutions. However, sometimes the algorithm occurs in local optima. The procedure may escape from this issue by accepting worse solution with some conditions. This new solution with a worse objective function value can be accepted as the current solution with a small probability determined by the Boltzmann function,e
x
p
(
-
Δ
/
K
T
), where K is a predetermined constant and T is the current temperature in Boltzmann function [18].The presentation of the TS algorithm to solve TTRPSD is explained as follows. At each iteration, the algorithm explores the solution space and produces a new solution from the neighbourhood of the current solution. Even if, the objective function value becomes worse with this new solution. A tabu mechanism is put in place to prevent the process from cycling over a sequence of solutions. An intuitive way to prevent cycles is to forbid the process from returning to previously encountered solutions, which is achieved by exploiting excessive bookkeeping. Some attributes of the past solutions are registered and any solution possessing these attributes may not be considered, and temporarily declared tabu for number of iterations (it is called tabu tenure).


### 3.1. Initial Solution

In this paper, this method is used to produce an initial solution considering the following assumptions. The initial solution can be used for M-SA, TS and MA. Consider(17)
P
∑
v
i
∈
r
k
ξ
i
<
2
Q
k
+
Q
r
≅
1
,
I
f
r
k
∈
P
V
R
,
C
V
R
,
P
∑
v
i
∈
r
k
ξ
i
<
2
Q
k
≅
1
,
I
f
r
k
∈
P
T
R
,
S
T
.All customers are classified as TCs and VCs. A string of numbers represents an initial solution, which is denoted by the set 1,2
,
…
,
n and N dummy zeroes. The parameter N is utilized to terminate the ST or the route. This parameter is only used to predict the number of routes or STs. If demand was deterministic and failure could not occur, this parameter could be ∑
i
q
i
/
(
Q
k
), where q
i is customer’s demand and Q
k is the truck capacity. However, since demand is stochastic and failure can occur a maximum of once, the parameter N can be computed by ∑
i
q
i
/
(
2
Q
k
) where · denotes the largest integer, which is equal to or smaller than the enclosed number and q
i is the expected value of customers’ demand. In the first n
+
N positions, the ith nonzero number denotes the ith customer to be serviced. VC can be serviced either by a complete vehicle or a single truck and the service vehicle type can be of 1 or 0. Therefore, the vehicle type for VC is 1 if it is serviced by a single truck. Otherwise, its vehicle type is 0 and must be serviced by a complete vehicle. Hence, TCs must be serviced only by a single truck; therefore there is no need to mention it in the solution [18].The presentation of the solution is explained more as follows. In this solution, the first number addresses the first customer that is to be served on the first route. A PTR is decided for the route, if a single truck is to service the first customer on a route. Without violating the following assumptionsP
{
∑
v
i
∈
r
k
ξ
i
<
2
(
Q
k
+
Q
r
)
}
≅
1 if r
k
∈ PVR and CVR or P
{
∑
v
i
∈
r
k
ξ
i
<
2
(
Q
k
)
}
≅
1 if r
k
∈ PTR or ST to represent the servicing sequence, from left to right, one by one, other customers are added to the route. Note that if it is a complete vehicle on the CVR main tour or on a PVR, the vehicle capacity may be (
Q
k
+
Q
r
) or on a CVR ST or on a PTR, if it is a single truck, it may be Q
k and it all depends on the vehicle type in use for the service. If, in the solution representation, the next customer to be served is zero, the vehicle will come back to the depot or the parking place. If it is on a CVR ST, the ST will be terminated, because the vehicle returns to the parking place. If not, it is on a CVR main tour, on a PVR or on a PTR. Consequently, the vehicle goes back to the depot and the route is terminated. It is worth mentioning that when recourse cost occurs, it should be considered and computed.In the solution representation, generation of a new route will be occurred, if the termination of previous route has been taken place and there are still customers yet to be serviced. Therefore, the next customer will be selected for a new route. A TTRP solution is given by this solution representation, without violating the assumption and it can be verified. On the other hand, by using this solution representation, the number of the vehicles that are used might exceed the available vehicles. As a result, after the solution representation has generated the routes, in order to decrease the number of the vehicles that are used a procedure for route combination is considered necessary. The procedure of the route combination, checks the possibility of combining the two existing routes. However, this combination process must not violate the constraint of the capacity of the vehicle and if there are routes that can be combined together without the violation of this constraint they will be merged without any modification. This process goes on until the number of the used vehicles is not greater than the number of the vehicles available or it will stop if there are not any more routes that can be combined together without violating the assumptions. If the outcome solution still persists on using more vehicles than available, for each additional trailer or truck that is used, a cost ofP, as a penalty, is added to the objective function in order to make the solutions of this type undesirable [9].


### 3.2. Simulated Annealing to Solve TTRPSD

Simulated annealing uses some parameters for improving the solutions, such asT
0, T
f, I
i
t
e
r, K, N
n
o
n, α, P, n
p
o
p, n
m
o
v
e, I
i
t
e
r
p
e
r
t
e
m
p, where T indicates the thermodynamic temperature, which is used to simulate the system at different temperatures, whereas T gradually cools down from an initial high temperature (T
0) to a final low temperature (T
f). It means that the procedure will be stopped when the temperature reaches lower than T
f. I
i
t
e
r represents the number of iterations during the procedure. K is the Boltzmann constant which is used in the probability function. N
n
o
n denotes the maximum number of allowable iterations in temperature. The algorithm will be terminated if the number of reductions exceed the N
n
o
n without improving the best cost. α is the coefficient for controlling the cooling scheme. P is the penalty cost associated with the number of extra vehicles used. Hence it is not allowed to use more than the number of available vehicles in 21 benchmark instance problems that are used for TTRPSD, the penalty cost is assigned to be very large (e.g., P
=
10
5 is considered for this model) to prevent from this issue. n
p
o
p is the number of initial solution that should be considered for producing new solutions. n
m
o
v
e is the number of move from current solution to generate new solutions. It means that each current solution can produce n
m
o
v
e new solutions and then the best one should be chosen as new solution. I
i
t
e
r
p
e
r
t
e
m
p is the number of iterations in each temperature T.In M-SA, it is important to define an appropriate method for finding the effective neighbors to improve an existing solution. A random neighborhood structure including swap, reversion, insertion, and change of service vehicle type of VCs is used for generating an appropriate neighborhood from a current solution. In each iteration, the algorithm generates a new solutionY from the current solution X by using swap, reversion, insertion and change of service vehicle type randomly.The swap is carried out by selecting two customers randomly and swapping them to generate a new solution from the current solution. The reversion is performed by selecting two numbers (customers) from the string of numbers representing the current solution then reversing the route from bigger number to smaller one. The insertion is performed by selecting two customers randomly and inserting the first customer immediately after the second one. The change of service vehicle type of VCs is fulfilled by selecting a vehicle customer randomly. If it is serviced by a single truck in the current solution, its service vehicle type will be changed to a complete vehicle and vice versa. It means that the vehicle service type of a selected VC will be changed from 1 to 0 or 0 to 1. For example, VC customer servicing in the main-tour by a complete vehicle will be serviced on a subtour by a single truck and vice versa.The probability of selecting swap, reversion, insertion and change service vehicle type methods will be set to 0.25. For each of them as it is assumed that each of these methods has an equal chance to improve the solution. In addition, failure may occur. Therefore, the failure should be diagnosed in the route and recourse action cost should be calculated for each failure.The M-SA procedure is as follows. At first, the current temperature is set to be the initial temperature and the algorithm generates the initial solutionsX
i randomly. In addition, the best solution X
i
b
e
s
t and the best objective function considering recourse cost, if occurred, are set to be X
i and o
b
j
(
X
i
,
P
), respectively. In each iteration, the algorithm generates new solutions Y
i from its neighborhood and evaluates its objective function value. Let Δ
i represents the differences between objective functional values from X
i and Y
i. Therefore, Δ
i
=
o
b
j
Y
i
,
P
-
o
b
j
(
X
i
,
P
). The probability of replacement of the current solution with the next solution is as follows(18)
1
,
Δ
≤
0
,
exp
⁡
-
Δ
i
K
T
,
Δ
>
0
.If Δ
i
≤
0, it means that the next solution is better than the current solution and should be replaced for the current solution. If Δ
i
>
0, it means that the next solution is not appropriate. However, as it mentioned, since escaping from trapping in local search, these kinds of solutions will be accepted with exp
⁡
-
Δ
i
/
K
T probability by generating a random number r
∈
[
0,1
] and replacing the solution X
i with Y
i if r
<
exp
⁡
-
Δ
i
/
K
T.After accomplishing the number of iterations, the current temperature will be decreased according to the formulaT
=
α
T, where 0
<
α
<
1. As there are more chances to find a better solution, the algorithm uses a local search procedure, which sequentially performs 2-Opt, swap, reversion, insertion and change of service vehicle types in every three temperature reductions [9].If the current temperatureT becomes lower than T
f or the number of reductions exceeds the N
n
o
n without improving the best cost, the algorithm will be terminated.


### 3.3. Tabu Search Algorithm to Solve TTRPSD

Similar to M-SA, TS also needs to start its procedure from an initial solution in the solution space. This initial solution can be infeasible. For this purpose, the initial solution method which is explained in previous section can be used for TS. TS algorithm tries to encourage the procedure to explore that part of the solution space which has not been visited yet. This purpose can be attained by preventing the reverse moves. The reversal of previous moves should become tabu for prohibiting a return to the previous solutions [39]. Cycle avoidance can be obtained over the tabu tenure. At each iteration, the inverse modification is identified as a tabu and occurred in the tabu list. In this paper, four operators (swap, reversion, insertion, and change of service vehicle type of VCs) have been considered for finding new solutions. In addition, tabu tenure ɳ is generated randomly with integer uniform distribution from [
r
,
s
]. Four tabu lists T
A
B
U
i
i
=
1,2
,
3,4 are considered to store the relevant inverse modifications for four operators. If the tabu status T
A
B
U
i
(
j
,
k
) is less or equal 0, it means that arc(
j
,
k
) is not tabu and vice versa.Tabu moves can be overridden if the aspiration criterion is satisfied. If the tabu solution has less objective function value, it can be overridden.The TS algorithm is explained as follows. At first, the relevant symbols are mentioned as follows:(i)
S: set of candidate solutions;
(ii)
t: iteration counter;
(iii)
t
n
o
n: current number of iteration without improvement;
(iv)
t
n: maximum number of allowable iterations;
(v)
t
m
a
x: maximum number of iterations;
(vi)
n
c: current number of generated candidate solutions;
(vii)
n
m
a
x: maximum number of candidate moves;
(viii)
x
∗: best-found solution for complete vehicle;
(ix)
y
∗: best-found solution for single truck;
(x)
T
x
,
y: the objective value;
(xi)
P
x
,
y: the penalized objective value.Tabu Search Heuristic
(1)
Generate an initial solution.
(2)
Initialize the parameters:A
1
=
A
2
=
1, A
3
=
A
4
=
1000, S
=
Ø, t
=
0, n
c
=
0, t
n
o
n
=
0, δ
i
=
0, n
m
a
x, t
m
a
x, t
n, r, s.
(3)
Update the data:x
→
x
∗, y
→
y
∗, T
x
,
y and P
x
,
y.
(4)
Ifn
c
≥
n
m
a
x, go to (3).
(5)
Generate random neighborhoodN
(
x
,
y
) and moves from current solution to new solution.
(6)
Updatet
+
1
→
t, n
c
+
1
→
n
c, S
=
S
∪
x
n
c
∪
{
y
n
c
}.
(7)
Evaluate all solutions inS by procedure 1 and set i
=
1, x
b
e
s
t
=
Ø, y
b
e
s
t
=
Ø and P
x
,
y
=
∞ and insert all nontabu solutions in N
’
(
x
,
y
).
(8)
Ifi
>
n
c, go to (10).
(9)
If(
x
i
,
y
i
)
∈
N
’
(
x
,
y
) and P
x
i
,
y
i
<
P
x
b
e
s
t
,
y
b
e
s
t, x
b
e
s
t
,
y
b
e
s
t
=
(
x
i
,
y
i
) and P
x
b
e
s
t
,
y
b
e
s
t
=
P
x
i
,
y
i; else if the solution is feasible, and T
x
i
,
y
i
<
T
∗
(
x
∗
,
y
∗
), some parameters should be updated as follows. x
b
e
s
t
,
y
b
e
s
t
=
(
x
i
,
y
i
), P
x
b
e
s
t
,
y
b
e
s
t
=
P
x
i
,
y
i and i
+
1
→
i and go to (8).
(10)
IfP
x
b
e
s
t
,
y
b
e
s
t
≥
P
∗
(
x
,
y
), set t
n
o
n
+
1
→
t
n
o
n, else update: x
b
e
s
t
,
y
b
e
s
t
→
(
x
,
y
), x
b
e
s
t
,
y
b
e
s
t
→
P
∗
(
x
,
y
).
(11)
Ifx
b
e
s
t
,
y
b
e
s
t is feasible and T
x
b
e
s
t
,
y
b
e
s
t
<
T
∗
(
x
∗
,
y
∗
), then set x
b
e
s
t
,
y
b
e
s
t
→
(
x
∗
,
y
∗
). In addition, T
x
b
e
s
t
,
y
b
e
s
t
→
T
∗
(
x
∗
,
y
∗
) and update the values of E
k, F
k, G
k, H
k and δ
i.
(12)
SetT
A
B
U
i
j
,
k
-
1
→
T
A
B
U
i
(
j
,
k
). It should be mentioned that the corresponding inverse modifications are tabooed for a tabu tenure ɳ, which is generated randomly from [
r
,
s
].
(13)
Ifm
o
d
t
,
10
=
0, then A
i should be updated. In addition, δ
i
=
0 and t
+
1
→
t.Ift
≤
t
m
a
x and t
n
o
n
≤
t
n, then go to (5); else the algorithm is terminated and the best solution found is (
x
∗
,
y
∗
).


### 3.4. Memetic Algorithm to Solve TTRPSD

To employ the MA, various crossovers, mutations and local search approaches were applied. The order crossover (OX) and partial-mapped crossover (PMX) were considered as the representations of permutation. Goldberg and Holland (1988) were the first ones to propose the PMX [40]. The operator of the PMX is an extended version of the two-point crossover. The two-point crossover is applied for a binary string and by using this operator, the occurrence of some infeasible solutions might be expected. However, some procedures are used by PMX for repairing these infeasible solutions by fixing the solutions resulted from the two point crossover [41, 42]. First, a random selection of two positions in the chromosome is performed and then the subchromosomes, which have been located in between these positions, are replaced. After that, the relations of the two sections of mapping are decided. Finally, the rest of the customers are randomly arranged according to the discovered relations in the remaining positions [40]. For the TTRPSD solution, this operator of PMX has been considered. Figure 1 depicts a PMX sample representation. This figure indicates that two samples, which are considered as the chosen parent chromosomes and according to PMX, two new offspring are generated. At first, the subchromosomes (the shadow customers) are selected. After that, all of the customers in this subchromosome are substituted and at the end, the rest of the customers are arbitrarily allocated consistent with the developed relations.Figure 1
Partial-mapped crossover for TTRPSD.The OX operator uses different methods for repairing the procedure. This operator has been successfully applied for the TSP and VRP [37]. First, a subchromosome from one parent is randomly selected. Then the new offspring is made by dropping the subchromosome into the same position. Finally, the remaining customers are arranged according to their positions in the other parent. A sample representation of an OX is illustrated in Figure 2.Figure 2
Order crossover for TTRPSD.Various mutation types such as swaps, insertions, inversions, displacements and changes in service vehicle type for the VCs have been used by MA in the production of offspring. By random selection of subchromosomes and putting them in a random situation, the implementation of the displacement mutation is performed. Swap is implemented by random selection of two customers and changing them for creation of new solutions out of the current solution. Reversion is performed by selection of two chromosome customers and overturning the route from the greater number to the smaller number. The achievement of insertion is obtained by selecting two random customers and inserting the first customer exactly after the second one. The fulfilment of the change in service vehicle type for the VCs is achieved by choosing a random VC. If in the current solution, it is serviced by a single truck, then its complete vehicle type will be changed to a service vehicle type and vice versa. This means that the selected VC vehicle service type will be changed from 0 to 1 or 1 to 0. For instance, a VC that has been serviced by a complete vehicle in the main tour will be serviced by a single truck on a ST and vice versa. Different types of mutations applied in the MA are illustrated in Figures3–6.Figure 3
Displacement mutation for TTRPSD.Figure 4
Swap mutation for TTRPSD.Figure 5
Reversion mutation for TTRPSD.Figure 6
Insertion mutation for TTRPSD.After applying the crossovers and mutation procedures, the approaches of LS are used for the improvement of the chromosomes in the candidates’ pool. Three different procedures of LS, change ST root, two-point exchange (TPE), and one-point movement (OMP) are applied in the MA. Each one of the approaches can be selected randomly with probability of equal value.A customer is moved in OPM from one particular route to one another. If this movement decreases the cost, then this movement can be accepted. In addition, only one customer is considered at a time. In movement execution, it is forbidden to move a TC customer to the main tour of the CVR or two PVRs. At first, the customers on the main tour and the PTR are examined by the algorithm and then the examination of the ST customers is started.Two customers of two different routes need to be replaced in a TPE. However, it is forbidden to move a TC customer into a PVR or into the main tour in the CVR. Considering the capacity of the truck, when customers are switched between two routes, all of the subroute continues to be feasible. In the aforementioned procedures, the position of the root node is never changed. This is due to an improvement in the solution when a number of root nodes are substituted. In this case, customer resequencing or root reselection is considered.


## 3.1. Initial Solution

In this paper, this method is used to produce an initial solution considering the following assumptions. The initial solution can be used for M-SA, TS and MA. Consider(17)
P
∑
v
i
∈
r
k
ξ
i
<
2
Q
k
+
Q
r
≅
1
,
I
f
r
k
∈
P
V
R
,
C
V
R
,
P
∑
v
i
∈
r
k
ξ
i
<
2
Q
k
≅
1
,
I
f
r
k
∈
P
T
R
,
S
T
.All customers are classified as TCs and VCs. A string of numbers represents an initial solution, which is denoted by the set 1,2
,
…
,
n and N dummy zeroes. The parameter N is utilized to terminate the ST or the route. This parameter is only used to predict the number of routes or STs. If demand was deterministic and failure could not occur, this parameter could be ∑
i
q
i
/
(
Q
k
), where q
i is customer’s demand and Q
k is the truck capacity. However, since demand is stochastic and failure can occur a maximum of once, the parameter N can be computed by ∑
i
q
i
/
(
2
Q
k
) where · denotes the largest integer, which is equal to or smaller than the enclosed number and q
i is the expected value of customers’ demand. In the first n
+
N positions, the ith nonzero number denotes the ith customer to be serviced. VC can be serviced either by a complete vehicle or a single truck and the service vehicle type can be of 1 or 0. Therefore, the vehicle type for VC is 1 if it is serviced by a single truck. Otherwise, its vehicle type is 0 and must be serviced by a complete vehicle. Hence, TCs must be serviced only by a single truck; therefore there is no need to mention it in the solution [18].The presentation of the solution is explained more as follows. In this solution, the first number addresses the first customer that is to be served on the first route. A PTR is decided for the route, if a single truck is to service the first customer on a route. Without violating the following assumptionsP
{
∑
v
i
∈
r
k
ξ
i
<
2
(
Q
k
+
Q
r
)
}
≅
1 if r
k
∈ PVR and CVR or P
{
∑
v
i
∈
r
k
ξ
i
<
2
(
Q
k
)
}
≅
1 if r
k
∈ PTR or ST to represent the servicing sequence, from left to right, one by one, other customers are added to the route. Note that if it is a complete vehicle on the CVR main tour or on a PVR, the vehicle capacity may be (
Q
k
+
Q
r
) or on a CVR ST or on a PTR, if it is a single truck, it may be Q
k and it all depends on the vehicle type in use for the service. If, in the solution representation, the next customer to be served is zero, the vehicle will come back to the depot or the parking place. If it is on a CVR ST, the ST will be terminated, because the vehicle returns to the parking place. If not, it is on a CVR main tour, on a PVR or on a PTR. Consequently, the vehicle goes back to the depot and the route is terminated. It is worth mentioning that when recourse cost occurs, it should be considered and computed.In the solution representation, generation of a new route will be occurred, if the termination of previous route has been taken place and there are still customers yet to be serviced. Therefore, the next customer will be selected for a new route. A TTRP solution is given by this solution representation, without violating the assumption and it can be verified. On the other hand, by using this solution representation, the number of the vehicles that are used might exceed the available vehicles. As a result, after the solution representation has generated the routes, in order to decrease the number of the vehicles that are used a procedure for route combination is considered necessary. The procedure of the route combination, checks the possibility of combining the two existing routes. However, this combination process must not violate the constraint of the capacity of the vehicle and if there are routes that can be combined together without the violation of this constraint they will be merged without any modification. This process goes on until the number of the used vehicles is not greater than the number of the vehicles available or it will stop if there are not any more routes that can be combined together without violating the assumptions. If the outcome solution still persists on using more vehicles than available, for each additional trailer or truck that is used, a cost ofP, as a penalty, is added to the objective function in order to make the solutions of this type undesirable [9].


## 3.2. Simulated Annealing to Solve TTRPSD

Simulated annealing uses some parameters for improving the solutions, such asT
0, T
f, I
i
t
e
r, K, N
n
o
n, α, P, n
p
o
p, n
m
o
v
e, I
i
t
e
r
p
e
r
t
e
m
p, where T indicates the thermodynamic temperature, which is used to simulate the system at different temperatures, whereas T gradually cools down from an initial high temperature (T
0) to a final low temperature (T
f). It means that the procedure will be stopped when the temperature reaches lower than T
f. I
i
t
e
r represents the number of iterations during the procedure. K is the Boltzmann constant which is used in the probability function. N
n
o
n denotes the maximum number of allowable iterations in temperature. The algorithm will be terminated if the number of reductions exceed the N
n
o
n without improving the best cost. α is the coefficient for controlling the cooling scheme. P is the penalty cost associated with the number of extra vehicles used. Hence it is not allowed to use more than the number of available vehicles in 21 benchmark instance problems that are used for TTRPSD, the penalty cost is assigned to be very large (e.g., P
=
10
5 is considered for this model) to prevent from this issue. n
p
o
p is the number of initial solution that should be considered for producing new solutions. n
m
o
v
e is the number of move from current solution to generate new solutions. It means that each current solution can produce n
m
o
v
e new solutions and then the best one should be chosen as new solution. I
i
t
e
r
p
e
r
t
e
m
p is the number of iterations in each temperature T.In M-SA, it is important to define an appropriate method for finding the effective neighbors to improve an existing solution. A random neighborhood structure including swap, reversion, insertion, and change of service vehicle type of VCs is used for generating an appropriate neighborhood from a current solution. In each iteration, the algorithm generates a new solutionY from the current solution X by using swap, reversion, insertion and change of service vehicle type randomly.The swap is carried out by selecting two customers randomly and swapping them to generate a new solution from the current solution. The reversion is performed by selecting two numbers (customers) from the string of numbers representing the current solution then reversing the route from bigger number to smaller one. The insertion is performed by selecting two customers randomly and inserting the first customer immediately after the second one. The change of service vehicle type of VCs is fulfilled by selecting a vehicle customer randomly. If it is serviced by a single truck in the current solution, its service vehicle type will be changed to a complete vehicle and vice versa. It means that the vehicle service type of a selected VC will be changed from 1 to 0 or 0 to 1. For example, VC customer servicing in the main-tour by a complete vehicle will be serviced on a subtour by a single truck and vice versa.The probability of selecting swap, reversion, insertion and change service vehicle type methods will be set to 0.25. For each of them as it is assumed that each of these methods has an equal chance to improve the solution. In addition, failure may occur. Therefore, the failure should be diagnosed in the route and recourse action cost should be calculated for each failure.The M-SA procedure is as follows. At first, the current temperature is set to be the initial temperature and the algorithm generates the initial solutionsX
i randomly. In addition, the best solution X
i
b
e
s
t and the best objective function considering recourse cost, if occurred, are set to be X
i and o
b
j
(
X
i
,
P
), respectively. In each iteration, the algorithm generates new solutions Y
i from its neighborhood and evaluates its objective function value. Let Δ
i represents the differences between objective functional values from X
i and Y
i. Therefore, Δ
i
=
o
b
j
Y
i
,
P
-
o
b
j
(
X
i
,
P
). The probability of replacement of the current solution with the next solution is as follows(18)
1
,
Δ
≤
0
,
exp
⁡
-
Δ
i
K
T
,
Δ
>
0
.If Δ
i
≤
0, it means that the next solution is better than the current solution and should be replaced for the current solution. If Δ
i
>
0, it means that the next solution is not appropriate. However, as it mentioned, since escaping from trapping in local search, these kinds of solutions will be accepted with exp
⁡
-
Δ
i
/
K
T probability by generating a random number r
∈
[
0,1
] and replacing the solution X
i with Y
i if r
<
exp
⁡
-
Δ
i
/
K
T.After accomplishing the number of iterations, the current temperature will be decreased according to the formulaT
=
α
T, where 0
<
α
<
1. As there are more chances to find a better solution, the algorithm uses a local search procedure, which sequentially performs 2-Opt, swap, reversion, insertion and change of service vehicle types in every three temperature reductions [9].If the current temperatureT becomes lower than T
f or the number of reductions exceeds the N
n
o
n without improving the best cost, the algorithm will be terminated.


## 3.3. Tabu Search Algorithm to Solve TTRPSD

Similar to M-SA, TS also needs to start its procedure from an initial solution in the solution space. This initial solution can be infeasible. For this purpose, the initial solution method which is explained in previous section can be used for TS. TS algorithm tries to encourage the procedure to explore that part of the solution space which has not been visited yet. This purpose can be attained by preventing the reverse moves. The reversal of previous moves should become tabu for prohibiting a return to the previous solutions [39]. Cycle avoidance can be obtained over the tabu tenure. At each iteration, the inverse modification is identified as a tabu and occurred in the tabu list. In this paper, four operators (swap, reversion, insertion, and change of service vehicle type of VCs) have been considered for finding new solutions. In addition, tabu tenure ɳ is generated randomly with integer uniform distribution from [
r
,
s
]. Four tabu lists T
A
B
U
i
i
=
1,2
,
3,4 are considered to store the relevant inverse modifications for four operators. If the tabu status T
A
B
U
i
(
j
,
k
) is less or equal 0, it means that arc(
j
,
k
) is not tabu and vice versa.Tabu moves can be overridden if the aspiration criterion is satisfied. If the tabu solution has less objective function value, it can be overridden.The TS algorithm is explained as follows. At first, the relevant symbols are mentioned as follows:(i)
S: set of candidate solutions;
(ii)
t: iteration counter;
(iii)
t
n
o
n: current number of iteration without improvement;
(iv)
t
n: maximum number of allowable iterations;
(v)
t
m
a
x: maximum number of iterations;
(vi)
n
c: current number of generated candidate solutions;
(vii)
n
m
a
x: maximum number of candidate moves;
(viii)
x
∗: best-found solution for complete vehicle;
(ix)
y
∗: best-found solution for single truck;
(x)
T
x
,
y: the objective value;
(xi)
P
x
,
y: the penalized objective value.Tabu Search Heuristic
(1)
Generate an initial solution.
(2)
Initialize the parameters:A
1
=
A
2
=
1, A
3
=
A
4
=
1000, S
=
Ø, t
=
0, n
c
=
0, t
n
o
n
=
0, δ
i
=
0, n
m
a
x, t
m
a
x, t
n, r, s.
(3)
Update the data:x
→
x
∗, y
→
y
∗, T
x
,
y and P
x
,
y.
(4)
Ifn
c
≥
n
m
a
x, go to (3).
(5)
Generate random neighborhoodN
(
x
,
y
) and moves from current solution to new solution.
(6)
Updatet
+
1
→
t, n
c
+
1
→
n
c, S
=
S
∪
x
n
c
∪
{
y
n
c
}.
(7)
Evaluate all solutions inS by procedure 1 and set i
=
1, x
b
e
s
t
=
Ø, y
b
e
s
t
=
Ø and P
x
,
y
=
∞ and insert all nontabu solutions in N
’
(
x
,
y
).
(8)
Ifi
>
n
c, go to (10).
(9)
If(
x
i
,
y
i
)
∈
N
’
(
x
,
y
) and P
x
i
,
y
i
<
P
x
b
e
s
t
,
y
b
e
s
t, x
b
e
s
t
,
y
b
e
s
t
=
(
x
i
,
y
i
) and P
x
b
e
s
t
,
y
b
e
s
t
=
P
x
i
,
y
i; else if the solution is feasible, and T
x
i
,
y
i
<
T
∗
(
x
∗
,
y
∗
), some parameters should be updated as follows. x
b
e
s
t
,
y
b
e
s
t
=
(
x
i
,
y
i
), P
x
b
e
s
t
,
y
b
e
s
t
=
P
x
i
,
y
i and i
+
1
→
i and go to (8).
(10)
IfP
x
b
e
s
t
,
y
b
e
s
t
≥
P
∗
(
x
,
y
), set t
n
o
n
+
1
→
t
n
o
n, else update: x
b
e
s
t
,
y
b
e
s
t
→
(
x
,
y
), x
b
e
s
t
,
y
b
e
s
t
→
P
∗
(
x
,
y
).
(11)
Ifx
b
e
s
t
,
y
b
e
s
t is feasible and T
x
b
e
s
t
,
y
b
e
s
t
<
T
∗
(
x
∗
,
y
∗
), then set x
b
e
s
t
,
y
b
e
s
t
→
(
x
∗
,
y
∗
). In addition, T
x
b
e
s
t
,
y
b
e
s
t
→
T
∗
(
x
∗
,
y
∗
) and update the values of E
k, F
k, G
k, H
k and δ
i.
(12)
SetT
A
B
U
i
j
,
k
-
1
→
T
A
B
U
i
(
j
,
k
). It should be mentioned that the corresponding inverse modifications are tabooed for a tabu tenure ɳ, which is generated randomly from [
r
,
s
].
(13)
Ifm
o
d
t
,
10
=
0, then A
i should be updated. In addition, δ
i
=
0 and t
+
1
→
t.Ift
≤
t
m
a
x and t
n
o
n
≤
t
n, then go to (5); else the algorithm is terminated and the best solution found is (
x
∗
,
y
∗
).


## 3.4. Memetic Algorithm to Solve TTRPSD

To employ the MA, various crossovers, mutations and local search approaches were applied. The order crossover (OX) and partial-mapped crossover (PMX) were considered as the representations of permutation. Goldberg and Holland (1988) were the first ones to propose the PMX [40]. The operator of the PMX is an extended version of the two-point crossover. The two-point crossover is applied for a binary string and by using this operator, the occurrence of some infeasible solutions might be expected. However, some procedures are used by PMX for repairing these infeasible solutions by fixing the solutions resulted from the two point crossover [41, 42]. First, a random selection of two positions in the chromosome is performed and then the subchromosomes, which have been located in between these positions, are replaced. After that, the relations of the two sections of mapping are decided. Finally, the rest of the customers are randomly arranged according to the discovered relations in the remaining positions [40]. For the TTRPSD solution, this operator of PMX has been considered. Figure 1 depicts a PMX sample representation. This figure indicates that two samples, which are considered as the chosen parent chromosomes and according to PMX, two new offspring are generated. At first, the subchromosomes (the shadow customers) are selected. After that, all of the customers in this subchromosome are substituted and at the end, the rest of the customers are arbitrarily allocated consistent with the developed relations.Figure 1
Partial-mapped crossover for TTRPSD.The OX operator uses different methods for repairing the procedure. This operator has been successfully applied for the TSP and VRP [37]. First, a subchromosome from one parent is randomly selected. Then the new offspring is made by dropping the subchromosome into the same position. Finally, the remaining customers are arranged according to their positions in the other parent. A sample representation of an OX is illustrated in Figure 2.Figure 2
Order crossover for TTRPSD.Various mutation types such as swaps, insertions, inversions, displacements and changes in service vehicle type for the VCs have been used by MA in the production of offspring. By random selection of subchromosomes and putting them in a random situation, the implementation of the displacement mutation is performed. Swap is implemented by random selection of two customers and changing them for creation of new solutions out of the current solution. Reversion is performed by selection of two chromosome customers and overturning the route from the greater number to the smaller number. The achievement of insertion is obtained by selecting two random customers and inserting the first customer exactly after the second one. The fulfilment of the change in service vehicle type for the VCs is achieved by choosing a random VC. If in the current solution, it is serviced by a single truck, then its complete vehicle type will be changed to a service vehicle type and vice versa. This means that the selected VC vehicle service type will be changed from 0 to 1 or 1 to 0. For instance, a VC that has been serviced by a complete vehicle in the main tour will be serviced by a single truck on a ST and vice versa. Different types of mutations applied in the MA are illustrated in Figures3–6.Figure 3
Displacement mutation for TTRPSD.Figure 4
Swap mutation for TTRPSD.Figure 5
Reversion mutation for TTRPSD.Figure 6
Insertion mutation for TTRPSD.After applying the crossovers and mutation procedures, the approaches of LS are used for the improvement of the chromosomes in the candidates’ pool. Three different procedures of LS, change ST root, two-point exchange (TPE), and one-point movement (OMP) are applied in the MA. Each one of the approaches can be selected randomly with probability of equal value.A customer is moved in OPM from one particular route to one another. If this movement decreases the cost, then this movement can be accepted. In addition, only one customer is considered at a time. In movement execution, it is forbidden to move a TC customer to the main tour of the CVR or two PVRs. At first, the customers on the main tour and the PTR are examined by the algorithm and then the examination of the ST customers is started.Two customers of two different routes need to be replaced in a TPE. However, it is forbidden to move a TC customer into a PVR or into the main tour in the CVR. Considering the capacity of the truck, when customers are switched between two routes, all of the subroute continues to be feasible. In the aforementioned procedures, the position of the root node is never changed. This is due to an improvement in the solution when a number of root nodes are substituted. In this case, customer resequencing or root reselection is considered.


## 4. Computational Study

In this section, a real case study from an Iranian dairy company is provided. This study has been carried out with the aid of Pegah Co., a large dairy distribution company in Iran, whose products are distributed to more than 50,000 retailers (customers) in Iran and some other countries. Iran Dairy Industries Co. (IDIC) is the largest dairy producer in Iran with “PEGAH” brand. This factory produces some dairy products such as Pasteurized and UHT milk, flavored milk, pasteurized and UHT cream, a variety of cheese (process, slice, pizza, UF), different kinds of yoghurt, probiotic products (such as yoghurt, cheese, ice cream), and drinking yoghurt.To implement the TTRPSD model for this case, 100 customers were selected based on their stochastic demands and the types of customers. To select the customers, the last 20 demands of each customer were realized and the customers of stochastic demands with the Poisson distribution were selected for this research. Then the customers’ locations were determined to compute the travel distance matrix between the customers and the depot. Furthermore, the type of each customer was distinguished and the truck customers (TCs) and the vehicle customers (VCs) were classified into their respective groups, where 30 customers were TC and the remaining customers were VC. The demands were measured based on their weights. The company considered 5 trucks and 3 trailers to serve these customers. The capacity of a truck and a trailer are 2000 and 3000 kilograms, respectively.


### 4.1. Computational Results

The MA, M-SA and TS were coded by MATLAB 7.9.0 using a computer with a 2.4 GHz dual processor and 4 G RAM. This type of problem has not yet been considered and there is no scope to compare this solution with an existing one. To overcome this issue, the special 21 TTRPSD benchmark instances are modified which are derived from Chao [7] for the TTRP with deterministic demands. First, the benchmarks were solved in order to increase the validity of the aforesaid algorithms and to show the consistency of the results. Then the case study problem was solved using MA, M-SA and TS. Further, the case study problem was checked by sensitivity analysis to confirm the results.To implement the benchmarks, the coordinates of vertices are kept the same as in the Chao’s instances [7]. As said earlier, the customer demands are stochastic and are generated with the Poisson distribution having a mean value equal to the customer demand. However, due to increasing possibility of failure on the route, the capacity of a truck and trailer is less than Chao’s benchmark, because the failure has rarely occurred with its capacities. Therefore, they are reduced considering P
ξ
i
≤
Q
k
≅
1 for each customer.Parameters used in the model may affect the quality of computational results. These parameters may be used for MA, M-SA or TS. In order to obtain better solutions, different values have been tested in initial experiments. These are:P
=
10000; α
=
0.99,0.975,0.95; I
i
t
e
r
=
20000,40000,60000,80000,100000; I
i
t
e
r
p
e
r
t
e
m
p
=
2,4
,
6,8
,
10; K
=
1,0.9,0.8,0.7,0.6,0.5; n
p
o
p
=
1,3
,
5,7
,
9; n
m
o
v
e
=
1,3
,
5,7
,
9; T
0
=
1000, T
F
=
0.1 or T
0
=
10, T
F
=
0.001; N
n
o
n
=
100,200; r
,
s
=
5,10; n
m
a
x
=
150; t
m
a
x
=
1000; t
n
=
50,100; P
c
=
0.7,0.75,0.8; P
m
=
0.2,0.3,0.4; P
l
=
0.3,0.4,0.5; n
=
50,100, where P
c, P
m, P
l and n are the probability of crossover, mutation, LS improvement and population size, respectively. Since it is not possible to use more than the predetermined number of vehicles, the penalty cost is considered too high. The parameters have been examined with different values and set to be α
=
0.95, I
i
t
e
r
=
60000, K
=
0.7, T
0
=
10, T
F
=
0.001, n
p
o
p
=
5, n
m
o
v
e
=
9, I
i
t
e
r
p
e
r
t
e
m
p
=
6 and N
n
o
n
=
200, t
n
=
50, P
c
=
0.7, P
m
=
0.3, P
l
=
0.5 and n
=
100. They seem to give the best results and will be used for further computation.Table1 illustrates the benchmark results. Each set has been run 10 times and the best solutions of MA, M-SA and TS from the 10 runs are shown.Table 1
The best solutions of TTRPSD benchmarks.
Problem ID
Number of VC
Number of TC
Truck  number
Truck capacity
Trailer number
Trailer capacity
MA results
M-SA results
TS results
1
38
12
5
60
3
60
635.45
639.38
643.30
2
25
25
5
60
3
60
671.84
659.77
666.56
3
13
37
5
60
3
60
677.44
683.76
686.42
4
57
18
9
60
5
60
929.11
931.63
942.32
5
38
37
9
60
5
60
921.01
919.10
924.55
6
19
56
9
60
5
60
971.90
991.28
1007.85
7
75
25
8
100
4
60
934.34
931.83
926.79
8
50
50
8
100
4
60
879.80
898.86
911.33
9
25
75
8
100
4
60
956.87
949.57
954.27
10
113
37
12
100
6
60
1117.32
1126.66
1133.38
11
75
75
12
100
6
60
1220.02
1198.94
1228.06
12
38
112
12
100
6
60
1221.29
1230.32
1243.02
13
150
49
17
100
9
60
1312.69
1309.84
1305.68
14
100
99
17
100
9
60
1374.32
1377.01
1381.58
15
50
149
17
100
9
60
1473.65
1475.00
1482.67
16
90
30
7
100
4
60
1092.37
1129.42
1113.86
17
60
60
7
100
4
60
1139.52
1114.28
1118.90
18
30
90
7
100
4
60
1118.93
1127.65
1120.05
19
75
25
10
100
5
60
843.46
843.46
852.14
20
50
50
10
100
5
60
853.72
867.77
876.15
21
25
75
10
100
5
60
867.62
965.24
948.41
Avg.
1010.13
1017.65
1022.25Table1 indicates that the average results generated by the MA are better than M-SA and TS. However, the differences between these results are insignificant and the results show that these algorithms are effective for solving related problems.Table2 indicates that the average results obtained by the proposed MA were improved about 1.27 and 0.78 percent comparing with tabu search and multipoint simulated annealing, respectively. Consequently, the performance of the proposed MA is found slightly better than SA and TS. As the differences between these results are insignificant, the results obtained by MA, M-SA and TS can be accepted as the new solutions. Therefore, the results indicate that the algorithms are efficient and effective in solving TTRP with stochastic demands.Table 2
Compared the results obtained from proposed MA with TS and M-SA.
Problem ID
Comparison of best solution between proposed MA and TS results (MA-TS)/TS (%)
Comparison of best solution between proposed MA and M-SA results (MA-SA)/SA (%)
1
-
1.22
-
0.61
2
0.79
1.83
3
-
1.31
-
0.93
4
-
1.40
-
0.27
5
-
0.39
0.21
6
-
3.60
-
1.96
7
0.81
0.27
8
-
3.58
-
2.12
9
0.27
0.76
10
-
1.44
-
0.84
11
-
0.66
1.73
12
-
1.78
-
0.74
13
0.53
0.22
14
-
0.53
-
0.20
15
-
0.61
-
0.09
16
-
1.93
-
3.28
17
1.81
2.21
18
-
0.10
-
0.78
19
-
1.02
0
20
-
2.56
-
1.62
21
-
8.52
-
10.11
Average
-
1.27
-
0.78For every parameter, the algorithms were tuned sequentially, leaving the remaining parameters unchanged. For this problem, in total about 400 runs were undertaken during the sensitivity analysis (including running a parameter setting 10 times on the problem for each algorithm).The results in kilometers which are obtained by MA, M-SA, and TS are presented in Tables3–5, respectively. Each set has been run 10 times and the best (best-sol), the worst (worst-sol) and the average solutions (avg.-sol) from the 10 runs are shown. Also, the time taken for the best solutions is presented in the last column. As this type of problem was not solved earlier, the results cannot be compared with any data or earlier heuristic solutions. However, the problems were checked by sensitivity analysis. Now, the results can be compared.Table 3
The MA results and sensitivity.
Best-sol
Worst-sol
Avg-sol
Sensitivity analysis
Time (second)
205.484
214.498
208.389
204.673
647Table 4
The M-SA results.
Best-sol
Worst-sol
Avg-sol
Sensitivity analysis
Time (second)
207.230
213.827
209.027
206.004
810Table 5
The TS results.
Best-sol
Worst-sol
Avg-sol
Sensitivity analysis
Time (second)
209.015
218.984
211.582
208.716
759Tables2 to 4 show the best solutions obtained by MA, M-SA and TS are 205.484, 207.230, and 209.015, respectively. In addition, the comparison between the obtained results and the sensitivity analysis results confirm a very slight difference. Therefore, the applied MA, M-SA, and TS are efficient with confidence consistency of a reasonable time.To indicate the convergence of the proposed approach, trends are shown in Figure7. This study has presented the relation between the number of iteration and the obtained objective function value. As it can be noted from the figure, the improvement rate of the solution reduces as the number of the iteration increases and after a particular number of iteration, the achieved solution converges. Therefore, the quality of the solution may not be enhanced by a greater number of iteration.Figure 7
Convergence trend for the algorithms’ solution.In brief, the analyses show that all of the methods of TS, M-SA, and MA are greatly efficient and their efficiency is almost equal; however, it seems that MA indicates slightly more efficiency than the others.


## 4.1. Computational Results

The MA, M-SA and TS were coded by MATLAB 7.9.0 using a computer with a 2.4 GHz dual processor and 4 G RAM. This type of problem has not yet been considered and there is no scope to compare this solution with an existing one. To overcome this issue, the special 21 TTRPSD benchmark instances are modified which are derived from Chao [7] for the TTRP with deterministic demands. First, the benchmarks were solved in order to increase the validity of the aforesaid algorithms and to show the consistency of the results. Then the case study problem was solved using MA, M-SA and TS. Further, the case study problem was checked by sensitivity analysis to confirm the results.To implement the benchmarks, the coordinates of vertices are kept the same as in the Chao’s instances [7]. As said earlier, the customer demands are stochastic and are generated with the Poisson distribution having a mean value equal to the customer demand. However, due to increasing possibility of failure on the route, the capacity of a truck and trailer is less than Chao’s benchmark, because the failure has rarely occurred with its capacities. Therefore, they are reduced considering P
ξ
i
≤
Q
k
≅
1 for each customer.Parameters used in the model may affect the quality of computational results. These parameters may be used for MA, M-SA or TS. In order to obtain better solutions, different values have been tested in initial experiments. These are:P
=
10000; α
=
0.99,0.975,0.95; I
i
t
e
r
=
20000,40000,60000,80000,100000; I
i
t
e
r
p
e
r
t
e
m
p
=
2,4
,
6,8
,
10; K
=
1,0.9,0.8,0.7,0.6,0.5; n
p
o
p
=
1,3
,
5,7
,
9; n
m
o
v
e
=
1,3
,
5,7
,
9; T
0
=
1000, T
F
=
0.1 or T
0
=
10, T
F
=
0.001; N
n
o
n
=
100,200; r
,
s
=
5,10; n
m
a
x
=
150; t
m
a
x
=
1000; t
n
=
50,100; P
c
=
0.7,0.75,0.8; P
m
=
0.2,0.3,0.4; P
l
=
0.3,0.4,0.5; n
=
50,100, where P
c, P
m, P
l and n are the probability of crossover, mutation, LS improvement and population size, respectively. Since it is not possible to use more than the predetermined number of vehicles, the penalty cost is considered too high. The parameters have been examined with different values and set to be α
=
0.95, I
i
t
e
r
=
60000, K
=
0.7, T
0
=
10, T
F
=
0.001, n
p
o
p
=
5, n
m
o
v
e
=
9, I
i
t
e
r
p
e
r
t
e
m
p
=
6 and N
n
o
n
=
200, t
n
=
50, P
c
=
0.7, P
m
=
0.3, P
l
=
0.5 and n
=
100. They seem to give the best results and will be used for further computation.Table1 illustrates the benchmark results. Each set has been run 10 times and the best solutions of MA, M-SA and TS from the 10 runs are shown.Table 1
The best solutions of TTRPSD benchmarks.
Problem ID
Number of VC
Number of TC
Truck  number
Truck capacity
Trailer number
Trailer capacity
MA results
M-SA results
TS results
1
38
12
5
60
3
60
635.45
639.38
643.30
2
25
25
5
60
3
60
671.84
659.77
666.56
3
13
37
5
60
3
60
677.44
683.76
686.42
4
57
18
9
60
5
60
929.11
931.63
942.32
5
38
37
9
60
5
60
921.01
919.10
924.55
6
19
56
9
60
5
60
971.90
991.28
1007.85
7
75
25
8
100
4
60
934.34
931.83
926.79
8
50
50
8
100
4
60
879.80
898.86
911.33
9
25
75
8
100
4
60
956.87
949.57
954.27
10
113
37
12
100
6
60
1117.32
1126.66
1133.38
11
75
75
12
100
6
60
1220.02
1198.94
1228.06
12
38
112
12
100
6
60
1221.29
1230.32
1243.02
13
150
49
17
100
9
60
1312.69
1309.84
1305.68
14
100
99
17
100
9
60
1374.32
1377.01
1381.58
15
50
149
17
100
9
60
1473.65
1475.00
1482.67
16
90
30
7
100
4
60
1092.37
1129.42
1113.86
17
60
60
7
100
4
60
1139.52
1114.28
1118.90
18
30
90
7
100
4
60
1118.93
1127.65
1120.05
19
75
25
10
100
5
60
843.46
843.46
852.14
20
50
50
10
100
5
60
853.72
867.77
876.15
21
25
75
10
100
5
60
867.62
965.24
948.41
Avg.
1010.13
1017.65
1022.25Table1 indicates that the average results generated by the MA are better than M-SA and TS. However, the differences between these results are insignificant and the results show that these algorithms are effective for solving related problems.Table2 indicates that the average results obtained by the proposed MA were improved about 1.27 and 0.78 percent comparing with tabu search and multipoint simulated annealing, respectively. Consequently, the performance of the proposed MA is found slightly better than SA and TS. As the differences between these results are insignificant, the results obtained by MA, M-SA and TS can be accepted as the new solutions. Therefore, the results indicate that the algorithms are efficient and effective in solving TTRP with stochastic demands.Table 2
Compared the results obtained from proposed MA with TS and M-SA.
Problem ID
Comparison of best solution between proposed MA and TS results (MA-TS)/TS (%)
Comparison of best solution between proposed MA and M-SA results (MA-SA)/SA (%)
1
-
1.22
-
0.61
2
0.79
1.83
3
-
1.31
-
0.93
4
-
1.40
-
0.27
5
-
0.39
0.21
6
-
3.60
-
1.96
7
0.81
0.27
8
-
3.58
-
2.12
9
0.27
0.76
10
-
1.44
-
0.84
11
-
0.66
1.73
12
-
1.78
-
0.74
13
0.53
0.22
14
-
0.53
-
0.20
15
-
0.61
-
0.09
16
-
1.93
-
3.28
17
1.81
2.21
18
-
0.10
-
0.78
19
-
1.02
0
20
-
2.56
-
1.62
21
-
8.52
-
10.11
Average
-
1.27
-
0.78For every parameter, the algorithms were tuned sequentially, leaving the remaining parameters unchanged. For this problem, in total about 400 runs were undertaken during the sensitivity analysis (including running a parameter setting 10 times on the problem for each algorithm).The results in kilometers which are obtained by MA, M-SA, and TS are presented in Tables3–5, respectively. Each set has been run 10 times and the best (best-sol), the worst (worst-sol) and the average solutions (avg.-sol) from the 10 runs are shown. Also, the time taken for the best solutions is presented in the last column. As this type of problem was not solved earlier, the results cannot be compared with any data or earlier heuristic solutions. However, the problems were checked by sensitivity analysis. Now, the results can be compared.Table 3
The MA results and sensitivity.
Best-sol
Worst-sol
Avg-sol
Sensitivity analysis
Time (second)
205.484
214.498
208.389
204.673
647Table 4
The M-SA results.
Best-sol
Worst-sol
Avg-sol
Sensitivity analysis
Time (second)
207.230
213.827
209.027
206.004
810Table 5
The TS results.
Best-sol
Worst-sol
Avg-sol
Sensitivity analysis
Time (second)
209.015
218.984
211.582
208.716
759Tables2 to 4 show the best solutions obtained by MA, M-SA and TS are 205.484, 207.230, and 209.015, respectively. In addition, the comparison between the obtained results and the sensitivity analysis results confirm a very slight difference. Therefore, the applied MA, M-SA, and TS are efficient with confidence consistency of a reasonable time.To indicate the convergence of the proposed approach, trends are shown in Figure7. This study has presented the relation between the number of iteration and the obtained objective function value. As it can be noted from the figure, the improvement rate of the solution reduces as the number of the iteration increases and after a particular number of iteration, the achieved solution converges. Therefore, the quality of the solution may not be enhanced by a greater number of iteration.Figure 7
Convergence trend for the algorithms’ solution.In brief, the analyses show that all of the methods of TS, M-SA, and MA are greatly efficient and their efficiency is almost equal; however, it seems that MA indicates slightly more efficiency than the others.


## 5. Conclusion and Further Research Directions

This research introduced, formulated, and solved the real-world stochastic TTRP problem. Indeed, in case of most real-life problems, demands are varying and TTRP should be expanded for meeting stochastic demands in the name of TTRP with stochastic demands (TTRPSD). Most researchers, however, solved TTRP with deterministic demands.This study modelled the TTRPSD problem within the framework of a stochastic program with recourse action (SPR). Three metaheuristic algorithms, namely, MA, M-SA, and TS have been applied to solve the problem. Firstly, twenty-one TTRP benchmarked problems have been modified for this model. Then MA, M-SA, and TS algorithms have been applied to solve these benchmarks to show the validity, consistency, and efficiency of the algorithms. The results are obtained in a reasonable time. As the differences between the results are not so much and the achieved solution converges after certain iterations, the aforesaid algorithms are efficient with confidence consistency.Each set of three case studies problem ran for 10 times and the best, worst, and average results were compared. As the differences between the best, worst, and corresponding average solutions are insignificant, the algorithm is found capable of producing TTRPSD solutions consistently and the results are useful.Furthermore, the problems have been tested by sensitive parameter analysis to realize the impact of the parameters. Since the results between the best, worst, and corresponding average solutions and the results obtained by sensitive parameter analysis are found closer; the results approved and indicated that the MA, M-SA, and TS are efficient and effective in solving this case study problem.Future research may attempt to extend the TTRPSD by introducing other practical or real-world conditions, such as multiple time windows, time dependent travel times and multidepots. Also, work on TTRP with stochastic travel times to solve real-life problems in manufacturing and transportation systems can be considered. Consequently, new TTRP benchmark instance problems need to be modified for this purpose.



---

*Source: 494019-2015-10-27.xml*