# Single-Commodity Vehicle Routing Problem with Pickup and Delivery Service

Authors: Goran Martinovic; Ivan Aleksi; Alfonzo  Baumgartner
Journal: Mathematical Problems in Engineering
 (2008)

---

## Abstract

We present a novel variation of the vehicle routing problem (VRP). Single commodity cargo with pickup and delivery service is considered. Customers are labeled as either cargo sink or cargo source, depending on their pickup or delivery demand. This problem is called a single commodity vehicle routing problem with pickup and delivery service (1-VRPPD). 1-VRPPD deals with multiple vehicles and is the same as the
single-commodity traveling salesman problem (1-PDTSP) when the number of vehicles is equal to 1.
Since 1-VRPPD specializes VRP, it is𝒩𝒫 hard in the strong sense. Iterative modified simulated annealing (IMSA) is presented along with greedy random-based initial solution algorithm. IMSA provides a good approximation to the global optimum in a large search space. Experiment is done for the instances with different number of customers and their demands. With respect to average values of IMSA execution times, proposed method is appropriate for practical applications.


---

## Body

## 1. Introduction

In this paper, we consider a real-world problem with a loaded
vehicle (the term capacitated vehicle is also used in
the literature). The vehicle is available for pickup and delivery
service. Single-commodity cargo, that is, same cargo type, is present at the
customer sites. The vehicle is able to carry out a
Pickup and delivery service of such cargo among the customers. By starting
and ending at the origin depot, the vehicle follows
a route without having a predefined sequence of Pickup and delivery services. This
problem is called aone-commodity vehicle
routing problem (1-VRPPD). The main feature of 1-VRPPD is that delivery
customers can be served with a cargo gathered from Pickup customers. This
problem was first defined in [1] as a one-commodity
Pickup and delivery traveling salesman problem (1-PDTSP), which is the
same as 1-VRPPD, when the number of vehicles is equal to 1. Although, 1-VRPPD
deals with multiple vehicles, in this paper we consider a single vehicle. In
order to solve this problem, the main task of our route planner is to calculate
a single suboptimal vehicle route, which starts and ends at the origin depot
and satisfies capacity constraints, minimizing the travel distance at the same
time. For the first time iterative modified simulated annealing (IMSA) is used
to solve this combinatorial optimization problem. In general, the VRP problem
and its variations are 𝒩𝒫 hard. 1-VRPPD is 𝒩𝒫 hard problem
in the strong sense since it specializes VRP. The main VRP feature that differs
from traveling salesman problem (TSP) is vehicle capacity constraint. Experimental
results given by IMSA confirm that complexity of VRP depends on the vehicle
capacity. When the capacity tends to infinity, VRP tends to TSP.The organization of this paper is as follows. In Section2, we
discuss several variations of VPR: the open VRP (OVRP), the distance-constrained capacitated VRP
(DCVRP), the VRP with backhauling (VRPB), the VRP with simultaneous Pickup and
delivery (VRPSPD), and several Hybrid variations of VRP, since they are most
similar to 1-VRPPD. Section 3 deals with problem and model definition. Solution
to the problem and its implementation are provided in Section 4. In Section 5 we describe
our simulation environment, as well as its input and output parameters. Solution
performance is illustrated in Section 6 with experimental results.


## 2. Related work

According to [2], over the last 30 years, less than 10% of the VRP
research was devoted to problems where vehicles have two stages, that is, where
vehicles are able to Pickup and deliver. The vehicle routing problem with
Pickup and delivery is not explored nearly as much as the other variants of
the vehicle routing problem.Recently, [3] described a practical example of heuristics applied
on VRP resulting in the transport network operational cost savings of 10% for
Melbourne mail distribution at Australia Post. A hybrid genetic algorithm was used to solve
VRP in [4] providing better solutions than those obtained using other
genetic algorithms.Thevehicle routing
problem with backhauling (VRPB) considers a vehicle servicing all delivery (Linehaul) customers with cargo loaded at
the depot, followed by Pickup (Backhaul)
customer services. Practical applications of this VRP variation are found in
grocery industry. A memetic algorithm was used in [5] to solve VRPB. In [6], the authors present a new tabu search algorithm that
starts from pseudo-lower bounds and is able to match almost all the best
published solutions and find many new solutions.Thevehicle routing
problem with simultaneous Pickup and delivery (VRPSPD) was introduced in
1989 by Min [7]. It considers both Pickup and delivery service at
each customer, while each collected cargo must be returned to the origin depot.
This problem is present in milk bottles transporting while empty ones must be
returned to the origin depot. A tabu search algorithm, with and without maximum
distance constraints, was recently developed in [8] to solve VRP-SPD.The main difference between VRPB, VRPSPD, and 1-VRPPD is
that in 1-VRPPD cargo picked up from Pickup customers can be delivered to
delivery customers. In 1-VRPPD, the predefined sequence of servicing customers
is not a constraint.Thedistance-constrained
capacitated vehicle routing problem (DCVRP) with a flexible assignment of
start and end depots is proposed in [9] and solved by means of a heuristic algorithm on the network
model.Theopen vehicle
routing problem (OVRP) considers solving a problem similar to VRP. The difference
is that the OVRP route ends when the last customer is served, that is, after
servicing customers the vehicle does not return to the depot. An overview of
OVRP algorithms is provided in [10]. A heuristic algorithm applied in [11] is comparable in terms of solution quality to the
best performing published heuristics. Practical applications of OVRP are in
home delivery of packages and newspapers with third-party contractors who use
their own vehicles and do not return to the depot.Hybrid variations of VRP with Pickup and delivery service have
also been present in literature. A hybrid variation of VRP is considered in [12] and it is used for blood transportation in health
care, where the vehicle route consists of the following sequence: delivery only, simultaneous Pickup and delivery, Pickup only. The proposed genetic algorithm (GA) called CLOVES
solves the problem in four steps.The 1-PDTSP, defined in [1, 13], is solved with a branch-and-cut algorithm. A special
case of the 1-PDTSP, when vehicle capacity is either 1 or ∞ and customers'
demands are either +1 or –1, was defined in
[14] and solved for a path and a tree graph topology.Practical applications of 1-VRPPD are found in the transport
of material of the same type where delivery customers correspond to cargo sinks
and Pickup customers correspond to cargo sources. That is the case in transportation
of gas, earth, sand, eggs, money, and so forth. In such a way, some customers produce
goods while others demand those goods.In this paper, we present a scenario with an employee utilizing
our VRP application. In order to find a feasible suboptimal vehicle route, the
user imports customer's coordinates and demands, and runs the IMSA algorithm.
Via graphical user interface, described in Section4, the user can select one
part of the resulting route in order to find a shorter route.In this paper, we explore the benefits of solving 1-VRPPD
with simulated annealing and its variations.


## 3. Problem and model definition

LetG be a
complete graph, V={v0,…,vn} a vertex set, E={{vi,vj}:vi,vj∈V} an edge set, A={(vi,vj):vi,vj∈V} an arc set. G has dual, directed and
undirected, representations, as Figure 1 illustrates. Vertices vi,i∈{1,…,n} correspond to the customers and vertex v0 corresponds to the depot vertex. Arc a∈A with start in vi and end in vj is denoted by (vi,vj). Non-negative cost cij
is associated with each
edge e={vi,vj}∈E that corresponds to the travel distance
between two vertices vi
and vj. Also, a binary
number xe∈{0,1} is assigned to e, where 1 represent visited edge and 0 represents unvisited edge.
Thus, cost matrix c is symmetric,
having the same cost in both directions, that is, cij=cji∀vi,vj∈V.
The use of loop edges is not allowed, that is, cii=+∞∀vi∈V.Figure 1
Dual, directed and
undirected, graph representation.Customer demandqi>0 indicates cargo that needs to be delivered, that is, the customer
requires delivery service. Customer demand qi≤0 indicates cargo that needs to be collected, that is, the customer
requires a Pickup service. When qi=0, Pickup service is considered in order to visit the costumer.For a given vertex setS⊆V,
let q(S)=∑i∈Sqi denote the total demand of the set. Let δ(S) and E(S) denote the set of
edges e∈E that have only one or both endpoints in S respectively.Letδ+(S):={(vi,vj)∈A:vi∈S,vj∉S}, ∀S⊂V,
denote the set of vertices j directly
reachable from vertex i. Let δ−(S):={(vi,vj)∈A:vi∉S,vj∈S}, ∀S⊂V,
denote the set of vertices j from
which vertex i is directly reachable.
Let δ(S):={{vi,vj}∈E:vi∈S,vj∉S}, ∀S⊂V,
denote the set of vertices j directly
reachable from vertex i. δ+(S) and δ−(S) are given with respect to A, while δ(S) is given with respect to E.
To ensure cargo conservation let the depot cargo
(3.1)q0=−∑i=1i=nqi,
and let
(3.2)K=∑∀vi∈V:qi>0qi=−∑∀vi∈V:qi≤0qi.
A vehicle with a fixed positive vehicle capacity Q is available at the depot. To ensure
feasibility, it is assumed that vehicle cargo qi∈[0,Q] is a positive whole number and never exceeds
its maximal value. An example of this model is illustrated in Figure 2.Figure 2
Vehicle routes of the one-commodity vehicle routing problem with Pickup and
delivery (1-VRPPD).Bearing in mind the aforementioned, a mathematical
formulation of 1-VRPPD can be defined as(3.3)Minimize∑e∈Ecexe(3.4)subjectto∑e∈δ(vi)xe=2∀vi∈V,(3.5)∑e∈δ(S)xe≥2∀S⊆V∖{v0},(3.6)∑e∈δ(v0)xe≤2M,(3.7)∑a∈δ+(vi)ga−∑a∈δ−(vi)ga=qi∀vi∈V,(3.8)0≤ga≤Q2xe∀a∈A,∀e∈E,(3.9)xe∈{0,1}∀e∈E.
Equation (3.4) ensures that Hamiltonian
cycle exists. With (3.5) a vehicle is forced to visit and leave
all customers along the route. (3.6) ensures that at least one of M vehicles visits and leaves depot
vertex. (3.7) describes the flow of Pickup and delivery commodities,
respectively. It ensures that cargo demands are appropriately cared for, that is,
the cargo is gathered or delivered in its entirety. (3.8) ensures that
vehicle capacity value is never less than zero and that it never exceeds its
maximal value. It defines a feasible 1-VRPPD route, where ga
is a continuous variable that has mathematical
meaning of a vehicle load going through an arc a. (3.9) represents xe
as a binary type variable assigned to e whose value is 1 if e is visited and
0 otherwise.In this section, the 1-VRPPD model is presented as a
combination of a generic VRP model described in [15] and added vehicle capacity constraints. References [1, 5, 13] describe similar VRP models with Pickup and delivery
services.


## 4. Proposed solution and implementation

Iterative modified simulated annealing (IMSA) algorithm is a
heuristic algorithm with iterative approach strategy used for solving
combinatorial optimization problems. IMSA and exact permutation algorithm (EPA)
are used to solve the 1-VRPPD𝒩𝒫
hard problem.In this section, we describe the implementation of IMSA
algorithm. The number of nodes, that is, customers, is denoted byn. The number of vehicles M is 1.


### 4.1. Exact Permutation Algorithm (epa)

For the purpose of improving the optimal route length for
small graph instances in which the number of nodes isn≤15,
we use Exact Permutation Algorithm (EPA). Instances with graphs, which have more
than 15 nodes, are not practical and solving EPA for an arbitrary n is impossible. Its computational
complexity is O(n!), therefore, it cannot be enumerated in polynomial time. EPA
represents the Exeter's
permutation algorithm described in [16]. EPA is used for improving existing solution via
interactive graphical user interface (GUI). Implementation of EPA is depicted at
the end of Section 5.


### 4.2. Simulated Annealing (sa)

Simulated annealing (SA) was first described in [17]. SA is a heuristic algorithm used for solving global
optimization problems [13, 15]. In the beginning of the algorithm, at high
temperature T, SA acts as a random
algorithm. Subsequently, the current feasible solution X changes to another feasible neighbor solution Y in a very random way. In our
implementation, neighbor feasible solution Y is created with nSAChanges random
feasible substitutions made on current feasible solution X. In this work, we consider nSAChanges=2n.
As system anneals proportionally with the number of iteration ϑ,
the algorithm becomes more deterministic. Criterion, which increases algorithm determinism,
as ϑ growths and T falls, is the main part of the SA algorithm and is defined with
temperature-probability function f(Δ) (4.1).
(4.1)f(Δ)={1,ifΔ≥0,eΔ/Tϑ,ifΔ<0,
where  Δ=RouteLength(X)-RouteLength(Y) is the amplitude defined with the difference between route
lengths of X and Y. Notation Tϑ∈[Tstart,Tend] represents the current temperature at iteration ϑ.
Temperature Tϑ falls linearly from Tstart to Tend as ϑ growths. Temperature-probability function f(Δ) results with continuous probability [0,1] for
rejecting new solutions. When Δ≥0,
the current feasible solution X is
replaced with better feasible solution Y.
Otherwise, random number r∈[0,1] is generated and compared with f(Δ).
If r≥f(Δ),
new feasible solution Y is rejected and
the current solution X remains
unchanged. If r<f(Δ),
then the current feasible solution X with
shorter route length is replaced with feasible solution Y with longer route length. Occasional
acceptance of a solution that leads to a longer route prevents the algorithm
from becoming stuck in a local minimum.


### 4.3. Modified Simulated Annealing (msa)

Modified simulated annealing (MSA) algorithm inherits
features of SA algorithm with a difference in creating neighbor solutionsY. MSA deals with multiple feasible neighbor
solutions Y made from the current feasible
solution X. Neighbor solutions Y are created in MSA nChangesmax inner iterations,
at the same temperature Tϑ,
by making nChangesmax random
substitutions on X. At the current temperature Tϑ,
the first inner iteration provides Y made
with 1 random substitution over X. The
second inner iteration provides another Y made with 2 substitutions over X. The
last, nChangesmax iteration
provides another Y created from X by making nChangesmax substitutions over X. In accordance with SA algorithm, at the same temperature more
neighbor solutions Y are created with
MSA algorithm. The current solution X is changed accordingly with the SA algorithm, for all Y solutions.


### 4.4. Iterated Modified Simulated Annealing (imsa)

(IMSA) is a heuristic
algorithm. IMSA algorithm inherits features of MSA algorithm. In order to find better
solution, IMSA algorithm has outer loop in which MSA repeatsnCountsmax times. Chances for
finding a better optimum are higher when the execution time needed for finding
suboptimal solution is longer or the number of iterations is greater. Algorithm 1
illustrates IMSA algorithm.Algorithm 1:(IMSA) algorithm.
1.lnitialize(nCountsmax, nChangesmax, Tstart, Tend, B,L,α)2.x=lnitialSolution( );3.while(Feasible(x) == FALSE )4.x=InitialSoiution( );5. minDist =Dist(x);6.for(nCounts = nCountsmax; nCounts > 0; nCounts ⋅⋅)7.for(nChanges = nChangesmax; nChanges > 0; nChanges ⋅⋅)8.for( k=0,T=Tstart; (k<B) && (T>Tend)); k++, T=T/α)9.for (z=0; z<L; z++)10.y=RandomNeighborSequence(x, nChanges);11.if(Feasible(y)== TRUE)12.Δ= minDist − Dist(y);13.if(p=exp(Δ/T)) > 1.0)14.p=1.0;15.if (RandomNumber([0,1]) <= p)16.if (Δ>0)17.x=y;18.minDist = Dist(x);Functions used in IMSA
algorithm, shown in Algorithm1, are:
Initialize: sets parameters to their default user defined
values.InitialSolution: returns an initial feasible route.Feasible: returns TRUE if the route is feasible,
otherwise returns FALSE.Dist: returns route length.RandomNeighborSequence: returns neighbor sequence Y that differs from the current sequence X on nChanges made.RandomNumber: returns a real number for a given interval.InitialSolution is
based on Greedy-random sequence (GRS) algorithm. In order to create an initial
solution, GRS starts from the start node, also denoted as current node. From
current node, all unvisited neighbor nodes, denoted as next node, are
considered with respect to their vehicle capacity Q constraint and the route length d. The feasible neighbors,
which satisfy vehicle capacity constraint, are sorted with respect to travel
length from the current node to the next node. The neighbour with the shortest
path would be selected with pure greedy algorithm. However, GRS algorithm
calculates initial solution from the graph. Starting with the current node, GRS
selects feasible neighbor nodes with respect to vehicle cargo. Again, the
nearest feasible neighbor would be selected with the pure greedy algorithm. But
GRS selects the nearest feasible neighbor solution with 80% probability and the
second feasible solution with 20% probability. This way, the greedy algorithm
becomes random.The values of parameters are the same for all variations of
SA algorithm. Starting temperatureTstart has value 1000, from which the system anneals until the temperature reaches
final temperature Tend=10.
Constant L with default value 100
represents the number of neighbor solutions Y considers at the same temperature T.
A bounded number of iterations B is set to be 1000 in order to prevent a long execution time if the input temperature interval [Tstart,Tend] is too wide. Temperature reducing factor α equals 1.1 in order to decrease the
temperature value by 10% in the next iteration, and is usually set to be a bit
greater than 1.In order to get the basic SA algorithm, the following
changes need to be implemented in the algorithm from Algorithm1. Line 6 does not
exist in the SA algorithm, and it should be removed. Accordingly, parameter Changesmax in line 10 should
be replaced with nSAChanges. Parameter Countsmax,
in line 7, should be set to value 1 or the line 7 should be removed,
respectively. The MSA algorithm can be extracted from IMSA illustrated in Algorithm 1, by removing the line 6 of the proposed algorithm. Thus, implementation of
the SA algorithm is simpler. It uses a fixed number of random substitutions nSAChanges, while MSA and IMSA have an
iteratively changeable number of changes from 1 to nChangesmax. MSA differs from IMSA in the number of Countsmax used to repeat the MSA algorithm. Finally,
a small solution space is explored with
SA, more is explored with MSA, and even more with IMSA algorithm. Parameters
that are different for different methods are illustrated in Table 1.Table 1
Parameters of implemented algorithms.
MethodSAMSAIMSAnSAChanges2n——nChangesmax—2n2nnCountsmax11100


## 4.1. Exact Permutation Algorithm (epa)

For the purpose of improving the optimal route length for
small graph instances in which the number of nodes isn≤15,
we use Exact Permutation Algorithm (EPA). Instances with graphs, which have more
than 15 nodes, are not practical and solving EPA for an arbitrary n is impossible. Its computational
complexity is O(n!), therefore, it cannot be enumerated in polynomial time. EPA
represents the Exeter's
permutation algorithm described in [16]. EPA is used for improving existing solution via
interactive graphical user interface (GUI). Implementation of EPA is depicted at
the end of Section 5.


## 4.2. Simulated Annealing (sa)

Simulated annealing (SA) was first described in [17]. SA is a heuristic algorithm used for solving global
optimization problems [13, 15]. In the beginning of the algorithm, at high
temperature T, SA acts as a random
algorithm. Subsequently, the current feasible solution X changes to another feasible neighbor solution Y in a very random way. In our
implementation, neighbor feasible solution Y is created with nSAChanges random
feasible substitutions made on current feasible solution X. In this work, we consider nSAChanges=2n.
As system anneals proportionally with the number of iteration ϑ,
the algorithm becomes more deterministic. Criterion, which increases algorithm determinism,
as ϑ growths and T falls, is the main part of the SA algorithm and is defined with
temperature-probability function f(Δ) (4.1).
(4.1)f(Δ)={1,ifΔ≥0,eΔ/Tϑ,ifΔ<0,
where  Δ=RouteLength(X)-RouteLength(Y) is the amplitude defined with the difference between route
lengths of X and Y. Notation Tϑ∈[Tstart,Tend] represents the current temperature at iteration ϑ.
Temperature Tϑ falls linearly from Tstart to Tend as ϑ growths. Temperature-probability function f(Δ) results with continuous probability [0,1] for
rejecting new solutions. When Δ≥0,
the current feasible solution X is
replaced with better feasible solution Y.
Otherwise, random number r∈[0,1] is generated and compared with f(Δ).
If r≥f(Δ),
new feasible solution Y is rejected and
the current solution X remains
unchanged. If r<f(Δ),
then the current feasible solution X with
shorter route length is replaced with feasible solution Y with longer route length. Occasional
acceptance of a solution that leads to a longer route prevents the algorithm
from becoming stuck in a local minimum.


## 4.3. Modified Simulated Annealing (msa)

Modified simulated annealing (MSA) algorithm inherits
features of SA algorithm with a difference in creating neighbor solutionsY. MSA deals with multiple feasible neighbor
solutions Y made from the current feasible
solution X. Neighbor solutions Y are created in MSA nChangesmax inner iterations,
at the same temperature Tϑ,
by making nChangesmax random
substitutions on X. At the current temperature Tϑ,
the first inner iteration provides Y made
with 1 random substitution over X. The
second inner iteration provides another Y made with 2 substitutions over X. The
last, nChangesmax iteration
provides another Y created from X by making nChangesmax substitutions over X. In accordance with SA algorithm, at the same temperature more
neighbor solutions Y are created with
MSA algorithm. The current solution X is changed accordingly with the SA algorithm, for all Y solutions.


## 4.4. Iterated Modified Simulated Annealing (imsa)

(IMSA) is a heuristic
algorithm. IMSA algorithm inherits features of MSA algorithm. In order to find better
solution, IMSA algorithm has outer loop in which MSA repeatsnCountsmax times. Chances for
finding a better optimum are higher when the execution time needed for finding
suboptimal solution is longer or the number of iterations is greater. Algorithm 1
illustrates IMSA algorithm.Algorithm 1:(IMSA) algorithm.
1.lnitialize(nCountsmax, nChangesmax, Tstart, Tend, B,L,α)2.x=lnitialSolution( );3.while(Feasible(x) == FALSE )4.x=InitialSoiution( );5. minDist =Dist(x);6.for(nCounts = nCountsmax; nCounts > 0; nCounts ⋅⋅)7.for(nChanges = nChangesmax; nChanges > 0; nChanges ⋅⋅)8.for( k=0,T=Tstart; (k<B) && (T>Tend)); k++, T=T/α)9.for (z=0; z<L; z++)10.y=RandomNeighborSequence(x, nChanges);11.if(Feasible(y)== TRUE)12.Δ= minDist − Dist(y);13.if(p=exp(Δ/T)) > 1.0)14.p=1.0;15.if (RandomNumber([0,1]) <= p)16.if (Δ>0)17.x=y;18.minDist = Dist(x);Functions used in IMSA
algorithm, shown in Algorithm1, are:
Initialize: sets parameters to their default user defined
values.InitialSolution: returns an initial feasible route.Feasible: returns TRUE if the route is feasible,
otherwise returns FALSE.Dist: returns route length.RandomNeighborSequence: returns neighbor sequence Y that differs from the current sequence X on nChanges made.RandomNumber: returns a real number for a given interval.InitialSolution is
based on Greedy-random sequence (GRS) algorithm. In order to create an initial
solution, GRS starts from the start node, also denoted as current node. From
current node, all unvisited neighbor nodes, denoted as next node, are
considered with respect to their vehicle capacity Q constraint and the route length d. The feasible neighbors,
which satisfy vehicle capacity constraint, are sorted with respect to travel
length from the current node to the next node. The neighbour with the shortest
path would be selected with pure greedy algorithm. However, GRS algorithm
calculates initial solution from the graph. Starting with the current node, GRS
selects feasible neighbor nodes with respect to vehicle cargo. Again, the
nearest feasible neighbor would be selected with the pure greedy algorithm. But
GRS selects the nearest feasible neighbor solution with 80% probability and the
second feasible solution with 20% probability. This way, the greedy algorithm
becomes random.The values of parameters are the same for all variations of
SA algorithm. Starting temperatureTstart has value 1000, from which the system anneals until the temperature reaches
final temperature Tend=10.
Constant L with default value 100
represents the number of neighbor solutions Y considers at the same temperature T.
A bounded number of iterations B is set to be 1000 in order to prevent a long execution time if the input temperature interval [Tstart,Tend] is too wide. Temperature reducing factor α equals 1.1 in order to decrease the
temperature value by 10% in the next iteration, and is usually set to be a bit
greater than 1.In order to get the basic SA algorithm, the following
changes need to be implemented in the algorithm from Algorithm1. Line 6 does not
exist in the SA algorithm, and it should be removed. Accordingly, parameter Changesmax in line 10 should
be replaced with nSAChanges. Parameter Countsmax,
in line 7, should be set to value 1 or the line 7 should be removed,
respectively. The MSA algorithm can be extracted from IMSA illustrated in Algorithm 1, by removing the line 6 of the proposed algorithm. Thus, implementation of
the SA algorithm is simpler. It uses a fixed number of random substitutions nSAChanges, while MSA and IMSA have an
iteratively changeable number of changes from 1 to nChangesmax. MSA differs from IMSA in the number of Countsmax used to repeat the MSA algorithm. Finally,
a small solution space is explored with
SA, more is explored with MSA, and even more with IMSA algorithm. Parameters
that are different for different methods are illustrated in Table 1.Table 1
Parameters of implemented algorithms.
MethodSAMSAIMSAnSAChanges2n——nChangesmax—2n2nnCountsmax11100


## 5. Simulation environment

For the purpose of experimental analysis, we developed a simulation
environment usingMSVisual C++ 6.0.
Figure 3 depicts functional components and the usage of the simulation
environment. During the initialization
phase in step 1, the graph with or without cargo, vehicle capacity Q, and common simulation parameters are
defined. In order to search for an optimal solution, the methods discussed in the
previous section are available in the search method phase, simulation step 2. When a suboptimal
solution is found, output parameters, route length d and execution time t, are evaluated in simulation step 3.Figure 3
Simulation environment.A solution
improvement method can be applied by using an interactive GUI. An interactive
GUI enables a user to improve the route performances in simulation step 4 (Figure3). By selecting a subset of vertices M⊆V on a route, the search improvement method provides some local improvements of
selected subroutes. For example, let M={a,b,c,d,e,f,g} be a set of user-selected vertices. Such set M represents an input parameter for the solution
improvement method that combines selected vertices in the resulting route. In
that way, neighbor solutions Y differ
from solution X only by vertices from
set M, as illustrated in Figure 4.Figure 4
Route improvement method based on an interactive GUI.In this paper, EPA that is defined in Section4, is used to
enumerate the routes with improved performance, if such routes exist. After the
suboptimal solution
phase in simulation step 3 (Figure 3), a suboptimal route can be improved by
selecting vertex set M and starting the
next simulation phase. Vertices v∈M from the suboptimal route are then permuted in
the solution improvement method phase in step 5 (Figure 3). If the new route is feasible and shorter than
the previous one, then it is set to be the route with improved performance.


## 6. Experimental results

We experimented with the simulation environment on a PC with
Intel T7200 processor under Windows XP. We chose several instances from [13] with positive demand at the depot. In this way, we
considered real world problems with the vehicle starting loaded or empty at the
depot. We compare the results with those from [13].Figures5–7 illustrate results
for the graph instance with 25 vertices and vehicle capacities 10, 15, and 20. Circles
represent vertices, and arrows represent the route. Numbers in vertices
represent the amount of product qi>0 that needs to be delivered, and qi≤0 indicates that the product needs to be collected. Shaded circle with
demand 0 is the depot vertex and corresponds to the vehicle start and end
positions. Figure 5 illustrates an optimal route for the smallest vehicle
capacity value Q=10. Figures 6 and 7 represent
the resulting routes for Q=15 and Q=20, respectively. The longest route
length corresponds to the smallest vehicle capacity. Larger vehicle capacity
results in a shorter route.Figure 5
IMSA route whenQ=10.Figure 6
IMSA route whenQ=15.Figure 7
IMSA route whenQ=20.Possible practical application of such a graph
representation is when there are several places with same sort of material storage
and other places that demand such material. For instance, the 1-VRPPD application may be used in newspaper redelivery
service. A newspaper company puts a newspaper in circulation, that is, they produce
certain amount of newspapers and deliver them to the newspaper stands. In the
morning, shops receive newspapers and sell them during the day with different
quantities. In the end of the day, newspapers are not appropriate for tomorrow
sales. In order to reduce company losses, shops that posses lesser newspaper
quantities require newspapers from shops that posses more newspapers in stock. A
small vehicle with capacityQ, has a
task to redeliver newspapers among newspaper stands.The heuristic algorithm presented
in Algorithm1 is appropriate for solving combinatorial optimization problems, which
may be 𝒩𝒫-hard. Besides 1-VRPPD, other problems can be solved with proposed
algorithm. Such problems should have feasible initial solution, which is later
improved when solution space is more explored. Inputs to the algorithm are: the
algorithm parameters, a graph with vertexes which represents customers and
their demands, and finally a depot vertex which represents start and goal
position. As a result, the proposed algorithm always returns a Hamiltonian
cycle with a suboptimal distance between vertexes.Another practical application is in
urgent medicament transport, when several hospitals have more medicaments of
one type than the others do. Practical transportation of H5N1 medicament is feasible
in an unpredictable epidemical contamination. Another practical application is
in concrete industry where customers are served with a truck which can load
concrete at several secondary depots, that is, at Pickup customers, and unload it
according to customer demands. Another practical application is a sensor
network, where sensors need to be transported in order to improve network
efficiency. This is the case in military and meteorology applications. Due to
some natural forces, a sensor network may be interrupted. Thus, our algorithm
may reduce network sensor transportation expenses.When vehicle capacityQ tends to infinity, 1-VRPPD requires less computational power than in cases
where vehicle capacity Q is small. According
to SA execution times shown in Table 2, feasible solutions are found quicker when
vehicle capacity Q is larger. Thus for
the same graph, when SA is used, execution times are highest when capacity
constraints are stringent, when Q is
smaller (Figure 5). As mentioned before, MSA has larger neighborhood space to explore
than SA. With the number of customers n>50,
the average route length and execution time is smaller when vehicle capacity is
larger. Consequently, IMSA average route length d is the smallest for largest Q.
According to Table 2, the time needed for
enumeration of the optimal feasible route is proportional to the size of
explored neighborhood space.Table 2
Experimental results.
nKQSAMSAIMSAd*drsd (%)t  (ms)d*drsd (%)t (ms)d*drsd (%)t (ms)153110266932748.193258228785.4410257127394.28833153115245629249.652232926487.6011232925556.30938153120244329037.993232926356.3412232924905.421081153125244328578.252232925816.6712227324895.431158153130242128409.303227325886.4013227324954.561175254510425450019.924392645658.2023373442586.021888254515392146117.044356842016.4723354139956.512088254520335245459.695327039208.0924333238827.4623132545253273444911.754341840087.0527339739076.4925622545303892459210.324320339238.5628327038637.142637406610595470427.313543362826.0731563960724.782782406615546065518.625512858505.9433512056505.382983406620525963837.584497756856.5936495354665.453251406625527364297.555466054965.8638458253966.013753406630520664059.465461854806.8341473852965.813934558410766286916.065694479446.5744669575015.943897558415640677338.426619272166.1146609767765.344173558420637975647.516578767776.3150568966155.884365558425656374337.006569164636.4755578663535.065392558430598072287.986579965576.0958550262625.9159517011310847797926.697842396307.6169789988826.4852407011315804291167.237742185286.1861686782267.6154797011320788791506.837725683715.6264713680265.3257787011325760087626.127718981825.4869706477835.4564737011330757888287.589700480486.3875678778656.076988Average values523462048.165487056186.6038478453945.843484Table2 presents our experimental results, where graph
instance features are denoted as follows. n is the number of vertices, K is
quantity of products that needs to be collected and picked-up (3.2), Q is vehicle capacity, d* is best solution
found, rsd is the route length's
relative standard deviation and d and t are average route length and average
execution time. Vehicle capacity Q has an experimental variable in range [10, 30]  in steps of 5. Proposed algorithm
is random based and average values obtained from 100 experiments.Route calculations with the IMSA last longer than with the MSA.
Thus, using IMSA is more appropriate in the case when the user working with our
application has customer demands several days before a vehicle starts with
Pickup and delivery services. Then the user can setup even a daylong search
strategy in order to get shorter routes. The numbernCountsmax (Algorithm 1), different in MSA and IMSA, needs to
be set appropriately in order to satisfy user needs. If the route length is
more important than the execution time, nCountsmax must be set to a relatively large value. Otherwise, if the execution time is
preferable, nCountsmax must be set to a relatively small value.In this paper, the IMSA is demonstrated as the algorithm resulting
in short routes, and the MSA as the algorithm with the short execution time.In Table3, dependencies between the number of customers and
average length and execution times for EPA, SA, MSA, and IMSA are proposed. For
instances, with 15 customers, EPA execution time is about 3 hours long, making
EPA impractical for larger instances.Table 3
Average route lengths and execution times
based on 100 experiments.
nd*dt (ms)EPASAMSAIMSAEPASAMSAIMSAEPASAMSAIMSA15227324862368235522732960266625541⋅107312103725—373834773455—464041233981—425229840—543049635006—656257595576—436334155—659860835954—773069916701—651475670—791774597151—913085528156—7685992Table3 illustrates the characteristics of methods analyzed in
this paper, average values of route lengths and execution times gathered in Table 2. In this way, IMSA is presented as the algorithm that has the shortest route length
and the longest execution time. The SA algorithm is the algorithm with the shortest
execution time and the longest route length.Dependency between average route lengthd and vehicle capacity is proposed in Figure 8. Figure 8 represents
dependencies between average execution time t and vehicle capacity Q. Figures 9(a) and
9(b) illustrate the values of average lengths and execution times with respect
of the number of customers.Average route lengthd (pixels) and average execution time t (ms) output parameters for various
graph instances with variable vehicle capacity Q values.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Best found route lengthd* (pixels), average route length d (pixels) and average execution time t (ms) output parameters for various
graph instances with variable n number of customers.
(a)(b)(c)The chart in Figures8(a), 8(c), 8(e), 8(g), 8(i) illustrate how a
change in vehicle capacity affects average route length d. The legend in Figure 8 illustrates which method is used and how
many customers are involved in a single route. The smallest instance with 15
customers is marked as SA_15, MSA_15 and IMSA_15, respectively. Other instances
with a different number of customers are marked accordingly. In the charts Figures 8(a), 8(c), 8(e), 8(g), 8(i), we see that the route length becomes shorter as the vehicle
capacity growths. This is the case for all methods applied. Figure 9(a)
illustrates the average route length values shown in Table 3. It shows that the
average route length depends linearly on a number of customers that needs to be
served. Route lengths are highest in case of SA, than there comes MSA, and
finally the smallest is for IMSA. This means that SA provides suboptimal routes
that are on average the longest suboptimal routes. IMSA routes are on average
the shortest ones. Resulting MSA routes are on average little longer than the
IMSAs, while they are much shorter than the SAs suboptimal routes.The charts in Figures8(b), 8(d), 8(f), 8(h), 8(j) illustrate how
vehicle capacity affects execution time t.
Since IMSAs execution lasts much longer than SA and MSA, the execution time for
IMSA is not displayed in Figures 8 and 9. As illustrated in Figure 9(b), the execution
time of SA is little shorter than of MSAs. The shortest execution times are on
average for the SA search method. MSAs execution times are little longer than
SA's, and much shorter than IMSA's, (see Table 3).Results proposed in Table4 compare our route length results
with heuristic algorithm used in [13]. H represents heuristic algorithm proposed in [13], d* is the shortest
route length, d is an average route
length, and t is the average execution
time. Used graph instances with Pickup and delivery demands are available at [18]. We chose instances with positive vehicle capacity at
the depot. These instances are presented in Table 4. Among tested instances, MSA
and IMSA algorithms have longer route lengths, while execution times are
shorter with MSA algorithm for several instances. As the number of customer
increases, compared results have proportional growth of average execution time.Table 4
IMSA algorithm output parameters compared with heuristic
algorithm [13].
NamenKQHMSAIMSAdt (ms)d*dt (ms)d*dt (ms)n30q20B3062205109905883671425577863682158n30q20E3066204916606545764123647370472101n40q20B40852053341106031692833591468172993n40q20I40902052621407299816033627479272922n50q20A50992059081807128821145704679854158n50q20C501122069624608616998442857594773857n60q20A601262066964308059932145776388894136n60q20B601462067304108663964144795993494012n30q40B3062404529505421635930513060382751n30q40E3066404822505530629930545161112795n40q40B4085405315906028669740580267193987n40q40I4090404967906026747641586571253791n50q40A50994058161406867806958684976915671n50q40C501124062841907400871148735985894493n60q40A601264061562406839839461672980725699n60q40B601464065242307643895954732886075199Average values…57081856874797341664376763795As illustrated in Table4, IMSA heuristics is the algorithm
that can solve 1-VRPPD and is comparable with heuristic H presented in [13]. MSA has 450% shorter execution time than H, while H provides
35% shorter route length. IMSA repeats MSA for 100 times. Thus, IMSA execution
time is approximately 100 times longer than MSA.


## 7. Conclusion

In this paper, we have presented the (IMSA) algorithm. IMSA is an iterative approach heuristic algorithm
used for solving the combinatorial optimization problem of 1-VRPPD. Both MSA
and IMSA inherit features of the simulated annealing (SA) algorithm with
different calculation of neighbor solutions.After the suboptimal route is calculated, several local
improvements along the route are possible. By using an interactive graphical
user interface, a user has the ability to select a sub graph on which solution
improvement method is applied. The exact permutation algorithm (EPA) was used
as a solution improvement method. In order to have the results in reasonable
time interval, selected sub graph should have less than 15 nodes, due to the
time complexity of EPA. Such improvements are useful because of large search
space of 1-VRPPD, since 1-VRPPD generalize VRP, which is already𝒩𝒫 hard. As the experimental
results substantiates, when the execution time is long enough or if the number
of iterations is large enough, MSA and
IMSA provide a good
approximation to the global optimum in large search space.



---

*Source: 697981-2009-03-25.xml*