2009 International Conference on New Trends in Information and Service Science

Two-Phase Particle Swarm Optimization

for Multi-depot Location-Routing Problem
PENG Yang
School of Computer & Information Engineering
Zhejiang Gongshang University
Hangzhou, P R China
[email protected]

CHEN Zi-Xia
School of Computer & Information Engineering
Zhejiang Gongshang University
Hangzhou, P R China
[email protected]
independently, and Min H,J.[7] synthesizes the past research
and suggests some future research directions for the LRP.
This paper divides the original LRP into two level
problems, i.e., one is the location-allocation problem and the
other is general vehicle routing problem. Each sub-problem
is then solved in a sequential and iterative manner by the
particle swarm optimization method embedded in the
general framework for the problem-solving procedure.
Particle swarm optimization(PSO) method is developed by
Kennedy and Eberhart [10], PSO is an evolutionary
algorithm that simulates the social behavior of bird flocking
to a desired place, which follows a collaborative populationbased search model, each individual of the population,
called a particle, flies around in a multidimensional search
space looking for the optimal solution. Particle, then, may
adjust their position according to their own and their
neighboring-particles experience, PSO starts with initial
solutions and updates them from iteration to iteration. PSO
has many advantages over other heuristic techniques such
that it can be implemented in a few lines of computer code,
it requires only primitive mathematical operators, and it has
great capability of escaping local optima. From our survey,
ours is the first attempt in proposing a PSO algorithm for
the Location-Routing problem.

AbstractThe location-routing problem is one of the

problems in distribution network designing and logistic
management, whose mathematical model is proposed in this
paper firstly. Since finding the optimal solution for this model
is an NP(non-polynomial)problem, we divide the original
problem into two sub-problems, i.e., one is the locationallocation problem and the other is general vehicle routing
problem. Each sub-problem is then solved in a sequential and
iterative manner by the particle swarm optimization method
embedded in the general framework for the problem-solving
procedure. From an illustrational example, the availability and
performance of proposed PSO can be approved.
Keywords- Particle swarm optimization; LocationRouting problem; logistics; Vehicle-routing
Introduction

The conceptual foundation of LRP studies date back to

more than thirty years ago, although those earlier studies are
far from capturing the total complexity of LRP, they already
recognized the close relation between location and
transportation decisions. In practical situation, the location
of facilities not only affects locational costs, but also has a
major impact on routing costs(Webb[2], Salhi and Rand
[3]), the overall problem is, in fact, a multilevel LRP that
should be taken jointly both location and routing problems
into consideration. Location-Routing Problem(LRP) is one
of the problems in integrated logistics optimization. Which
is defined to find the optimal number and locations of the
facilities, and simultaneously to get the vehicle schedules
and distribution routes so as to minimize the total system
cost (Tai-His Wu,C,J [1]). In the last years, however,
several works have addressed LRP of various
characteristics. Two sub-problems, LAP (location allocation problem ) and VRP (vehicle routing problem) of
LRP have been recognized as NP-hard, thus LRP also
belongs to the class of NP-hard problem. Since the exact
algorithm seems infeasible to solve the problem, it is
worthwhile to develop some heuristic methods. For
instance, Nagy and Salhi[4,5] adopted the concept of nested
methods for LRP; Tai-His Wu[1] proposes Simulation
annealing method which decomposes the LRP into a LAP
and a VRP; Maria A S, Juan A [8] also present a two-phase
tabu search architecture for the solution of the LRP. Many
researchers(Chen TW[6],S.c.Liu,S.B Lee[9], etc.)resort to
iterative heuristic approaches, in which at each iteration a
location problem plus a routing problem are solved

978-0-7695-3687-3/09 $25.00 2009 IEEE

DOI 10.1109/NISS.2009.223

I.

DESCRIPTION OF LRP AND ITS

MATHMATICAL MODEL

Location-Routing Problem can be stated as follows.

A feasible set of potential facility sites and location and
expected demands of each customer are given, which will
meet its demands. The shipments of customer demand are
carried out by vehicles, which are dispatched from the
facilities and operated on routes that include multiple
customers. There is a fixed cost that is associated with
opening a facility at each potential site. A distribution cost
associated with any routing of vehicles includes the cost of
delivery operations. The cost of delivery operations is linear
in the total distance the location of the facilities and the
vehicle routes from the facilities to customers, the objective
is to minimize the sum of the location an distribution costs
so that the facilities capacities are not exceeded. The
hypotheses are as follows:
(1) The transportation is just in time.
(2) The facility is both starting point and destination of
circular vehicle routing; each facility serves more than
two customers.
240

(3)
(4)
(5)
(6)

Nature of demand/supply is deterministic.

There are multiple facilities.
Size of vehicle fleets is a single vehicle.
Vehicle capacities are determined. The total amount
of goods is limited in every route by each vehicles
capability.
(7) facility capacities are undetermined; not all facilities
have been chosen in every decision.
(8) Each customer is served by one and only one vehicle.
(9) Each facility is considered as a separate entity, not
linked to the other facilities.
(10) The objective is to minimize total cost.
The sets, parameters and variables used in the
mathematical model are defined as follows:
G={i|i=1,,M}: the set of M feasible sites of candidate
facility;
H={j|j=1,N}: the set of N customers to be served;
S = {G} {H } :the set of all feasible sites and
customers, it is also referred to all nodes;
V={k|k=1,,K}: the set of K vehicles available for
routing from the facilities;
Cij: distance between node i to node j, i S , j S
Gi: fixed costs of establishing and operating a depot at
site r(r=1,,m);
Fr: fixed costs of using vehicle k(k=1,,K);
qj: average demands of customer j(j=1,2,,N);
Ui: maximum throughput at depot i;
Qk: capacity of vehicle k(k=1,,K);
Decision variables:
Xijk=1, traveling by vehicle k from node i to customer j ; 0
otherwise.
Yi=1, depot i is established or opened; 0 otherwise.
Mathematical Model for LRP is defined as follows.

II.

kV

(1)

iG jH

subject to

X pjk = 0, k V , p S

ijk

= 1,

j H

(2)

kV iS

iH jH

ipk

iS

X ijk Qk ,

k V

(3)
(4)

jS

ijk

1,

k V

(5)

pqk

+Y p +Y q 2 p, q G

(6)

ijk

Yi 0, i G

TWO-PHASE METHOD FOR LOCATION-ROUTING

PROBLEM

Like evolutionary algorithm, PSO conducts search using

a population(called swarm) of individuals(called particle)
that are updated from iteration to iteration. Each particle
represents a candidate position(i.e., solution) to the problem
at hand, resembling the chromosome of GA. A particle is
treated as a point in an M-dimension space, and the status of

kV

(8)

B. Fundamental principle of PSO

rG jH

k V , i G

A. Two-phase approach
In this paper a heuristic method which decomposes the
MWLRP into the following two sub- problems is proposed:
Phase 1. Location-allocation problem (LAP). Phase 2.
Vehicle Routing problem (VRP).These sub-problems can be
solved by either optimization techniques or heuristic
methods iteratively and at the same time, the dependence
between each sub-problem can be considered. In each
iteration, current routes are unified to include more
customers and allocated to selected.
warehouses while
taking into account the capacity restrictions of the
warehouses. Since both LAP and VRP are well-known
combinatorial optimization problems, heuristic methods that
give quick and good solutions have been presented. In the
1rst iteration, the solutions of the LAP are some selected DC
sites (only necessary ones, not all) and a plan for allocating
customers to each chosen DC site. Notice that the calculation
of the distance for the LAP in this iteration (the "rst) is based
on the `moment-suma equation used by the traditional
location problems. These solutions are then used as input to
the VRP to generate a starting feasible set of routes. At this
moment it is very possible that the number of DCs
established can still be reduced. To achieve this objective,
each current route consisting of several customers is viewed
as a single node with the demand represented by the sum of
demands of all customers in that route. These `big nodesa are
then consolidated for reducing the number of DCs
established and, thus, the total cost. The consolidation
process starts from the second iteration and is performed by
the LAP module, followed by the VRP module. This
procedure is repeated until the convergence criterion is
met.Maintaining the Integrity of the Specifications.

iG jH kV

+ Fk xijk )

Yi 0,

In this model, the objective function minimizes the total

cost of routing, establishing and operating the facilities.
Constraint (2) ensures that each customer is served by one
and only one facility. Constraint (3) ensures that the total
demands served by one facilities do not surplus its capacity,
while (4) is the route continuity constraints, which implies
that the vehicle should leave every point entered by the same
vehicle. Constraint (5) guarantees that each vehicle is routed
from one facility. Constraint (6) guarantees that there is no
link between an two facilitys, Constraints (7) and (8) require
that a vehicle only be from a facility if that facility is opened.

f ( x) = min( Gi y i + C ij X i jk
iG

ijk

jH

(7)

kV jH

241

a particle is characterized by its position and velocity.

Initialized with a swarm of random particle, PSO is achieved
through particle flying along the trajectory that will be
adjusted based on the best experience or position of the one
particle(called local best) and the best experience or position
ever found by all particles(called global best). The Mdimension position for the ith iteration can be denoted as
Xi(t)={xi1(t),xi2(t),xiM(t)},similarly, the velocity(i.e.,distance
change),also a M-dimension vector, for the ith interationcan
be described as Vi(t)={vi1(t),vi2(t),,viM(t)}, the particleupdating mechanism for particle flying(i.e., search process)
can
be
formulated
as
follows:

rij = ( xi x j ) 2 + ( yi y j ) 2 , i G, j H
1) if

i, and then the customers served by facility i ,

Fi=(H1,H2,..,Hr)
2) if

j =1

k =1

d j > qk or

> U i , then transfer

j =1

some customers to other sub-nearest facility, until the

restrictions of capacity are satisfied.
(1) Add move of location
For a given number of facilities, the location phase can
perform add moves, which opens one more facility randomly
from other closed facilities. This operation will increase the
fix cost of open facilities, but it also influence the
distribution cost, so the total cost possibly decrease.
If perform add move of location, firstly open a facility
randomly, and then re-allocation the customers to new set of
facilities as the operations of the initialization, which realize
a updating of the location configuration.
(2) Delete move of location
For a given number of facilities, the location phase can
also perform add moves, which closes one facility randomly
from now open facilities, the operation influence the total
cost obviously. Which firstly closes a facility randomly, and
then re-allocation the customers to new set of facilities as
add move operation, which also modifies the location
configuration.

Vi (t ) = w(t )Vi (t 1) +c 1 r1 ( X iL X i (t 1))

+ c 2 r2 ( X G X i (t 1)
X i(t ) = Vi (t ) + X i (t 1)

rij Ri then set the customer j served by facility

(11)

(12)
Where i=1,2,,P, and P means the total number of the
particles in a swarm, which is called population size;
t=1,2,,T, and T means the iteration limit;
L
X iL = {xiL1 , xiL2 ,L, xiM
} represents the local best of the
ith particle encountered after t-1 iterations, while
X G = {x1G , x 2G ,L, x mG } represents the global best among
all the swarm of particles achieved so far. c1 and c2 are
positive constants(namely, learning factors), and r1 and r2 are
random numbers between 0 and 1; w(t) is the inertia weight
used to control the impact of the previous velocities on the
current velocity, influencing the trade-off between the global
and local experiences.

D. Routing phase
In each iteration, followed by the location phase, i.e. ,
given a current facilities configuration, a routing phase is
started. This paper presented a novel particle swarm
optimization algorithm for the VRP.
(1) The particle representation
In this section, we describe the formulation of a PSO
algorithm for the VRP, one of the key issues in designing a
successful PSO algorithm is the representation step, i.e.
finding a suitable mapping between problem solution and
PSO particle. In this paper, we reference the method in
literatures[11,12] and setup a 2N-dimension search space, N
is the number of customer to be served, 2N-dimension
particle X can be regarded as two N-dimension vectors: the
elements of Xv stand for the customers are served by
different vehicle, and the element of Xr stand for the serves
priority in the route of relevant vehicle. For example, in a
given location configuration, there are 4 vehicle in depot i
and 10 customers to be served, if a solution can be
represented as the particle in table 1, then we should turn the
Xv into integer, and if Xvi=Xvj(i,j=1,2,10),sort the order
according to Xr, which is showed in table 2.

C .location phase
In the location phase, our purpose is to determine a good
configuration of facilities to used in the distribution. For
each of the location configurations visited during the
location phase, the routing phase is started in order to update
the routing according to the new configuration. Since only a
certain section of the routing is affected by the change in the
depot configuration, it is possible to restrict the search to
routing, thus, the routing phase is a localized search, as
opposed to a global exploration of all routing moves, this
eliminates a lot of unnecessary computation, and allows the
two phase algorithm to find good solutions within
reasonable computation time. Now, we define the element
of location phase.
(0) initialization of location:
Suppose the positions of m candidate facilities are
PF1(x1 ,y1), PF2(x2,y2),,PFm(Xm,ym), the distance of any
PFi(xi ,yi) and PFj(xj ,yj) is dij:

d ij = ( xi x j ) 2 + ( yi y j ) 2
Ri is the radius of cover area of facility i,

Ri = 1 / 2 min(d ij | j G {i}) i G

The set of customers H=(H1,H2,,Hn), the position of

Hj=(xj,yj),

242

A PARTICLE REPRESENTATION FOR SOLUTION

TABLE I.

Step2: Repeat until the number of generation equal to

maximum generation number or satisfy the stop situation.
1) For each particle, update X and V according to
equation (2) and equation (3), the updated particle position
must be subject to the limit[Xmin,Xmax]=[1,M], each element
of the 2N-dimension particle that is beyond [Xmin,Xmax] can
be adjusted as: if Xi(t)>Xmax, then Xi(t)=Xmax ; else if
Xi(t)<Xmin, then Xi(t)=Xmin; and each element of the 2Ndimension particle velocity that is beyond [-Vmax,Vmin]=[N,N] should be adjusted as: if Vi(t)>Vmax then Vi(t)=Vmax,
else if Vj(t)<-Vmax, then Vj(t)=-Vmax.
2) Turn the 2N-dimension particle X into appropriated
representation as section III(B).

customer 1 2 3 4 5 6 7 8 9 10
Xv
1. 2 2. 1 1. 1 2. 5 3. 2 2. 2 3. 7 2. 7 3. 1 1. 2
Xr
1. 2 3. 4 5. 6 8. 1 2. 3 4. 6 2. 1 6. 2 4. 8 8. 1
SHORT THE ORDER OF THE SOLUTIN IN

TABLE II.

PARTICLE

customer 1
Xv
1
Xr
1

2
2
1

3
1
2

4
2
4

5
3
2

6
2
2

7
3
1

8
2
3

9 10
3 1
3 3

3) Evaluate the fitness of all particle in the new

population;
4) update the Gbest, if one particles fitness is better than
current Gbest; for each particle, if the fitness is better than its
Pbest, then update it;
Step 3: Stop the PSO and get the approximately optimal
solution
The PSO execution should be terminated if the current
iteration meets any one of the stop signals. The stop signals
considered here include: (1) maximum number of iterations
with steady global best(i.e., maximum number of iterations
since last updating of the global best) and (2) maximum
total number of iterations. Then the optimal solution to LRP
has been got, which corresponds to the global best particlerepresentation.

The particle represents that related routes for vehicle

are: route 1(vehicle 1): 0-1-3-10-0; route 2(vehicle 2): 0-26-8-4-0; route 3(vehicle 3): 0-7-5-9-0 , where 0 denote the
depot, which is the begin and end of the route.
This representation makes each customer can be served by
only one routing/vehicle, and number and amount of
selected vehicle also can be embodied. Though number of
dimension is increased, PSO shows the good characteristic
in multi-dimension optimization, the computation
complication does not increase observably.
(1) the modified objective function
LRP is a restricted combined optimization problem,
formulation 3 scripts the restriction of route capacity, that is,
the sum of demand to be served of one vehicle can not
surplus its capacity. This paper adopts a modified objective
function that includes a penalty term associated with
infeasible, R is a large number used as penalty coefficient,
which multiples the overload, then is added to the objective
function. The modified objective function as follow:

E. Total procedure of proposed method

Total procedure of proposed method can be scribbled as
follows:
Step 1:finding the initial location configuration
Perform the initialization operation of location phase;
Step 2: finding the routing of now location configuration
Perform the routing phase using our proposed PSO
algorithm, which is stated as section 3.4.
Step 3 Calculate the total costs by original objective
function which is stated in formula (1), records the current
best costs and the solution.
Step 4: updating the location configuration
For the given location configuration, perform randomly
the add moves of location or delete moves, the numbers of
add moves and delete moves of location can be generated
randomly from [0,P], P is the maximum number of
permission to add or delete the facility.
Step 5: return to step 2 until the number of generation
equal to maximum generation number or satisfy the stop
situation, which can be considered that maximum number of
iterations since last updating of the best solution.

F ( x) = min( Gi y i + C ij X i jk +
iG

F x
k

kV

iG jH

iG jH kV

ijk

i =1

i =1

+ max( C ij X ijk q,0))

(2) Procedure of the PSO algorithm for VRP

The procedure of the PSO can be stated as follow:
step1: initialize the particle swarm,
1) For each 2N-dimension particle in the population, the
position vector Xv={X1,X2,,XN}, each element is randomly
generated from 1~M, and Xr={XN+1,XN+2,,X2N} randomly
from 1~N;
2) Each velocity vector is also 2N-dimension,
Vv={V1,V2,,VN}, each element is randomly generated from
-M~M, and Vr={VN+1,VN+2,,V2N} randomly from -N~N;
3) Evaluate the fitness of each particle in the population,
the fitness function is described as formula (13)
4) Initialize the Gbest with the best fitness (lowest cost)
among the population.
5) Initialize the Pbest with a copy the population;

III.

4.EXPERIMENTAL RESULT

Numerical experiments were conducted to examine the

computational effectiveness and efficiency of our proposed
method. The heuristic methods are coded using the Visual
243

C++ programming language and the tests are carried out on a

PC Pentium 1.8GHZ.
The demand for each customer in any test is randomly
selected from a uniform distribution U[450, 600] for each
month. The demand during lead time for each customer is
randomly selected from a uniform distribution U [0, 10].
The location (x, y) of each customer and candidate depot is
randomly selected from a uniform distribution U[0, 100].
The other relevant cost values are listed in table and
some parameters used in PSO are listed in table .

other heuristic method described earlier, which is presented

in literature 13. It is found that our method is better than
Tabu-Search based method in total costs. And also, when
the number of candidate depots increases, the total system
costs for all the two solutions decrease. When the number of
customers increases, the total system costs for all the three
heuristic solutions increase.. However, in terms of average
CPU time. these two methods are efficient, for all average
CPU times in Table 3 are less than 30 seconds.
IV.

TABLE III. OTHER RELEVANT COST VALUES USED IN THE

Vehicle type
Vehi cl e capaci ty
Fixed cost of vehicle using

Type 1

Type 2

Type 3

100

60

40

Distribution Cost ratio

4 / Ton. Km

Depot establishing Cost

1200(high), 700(low)

TABLE IV.

PARAMETERS USED IN PSO FOR ROUTING PHASE

W(t)

C1

C2

Population
Size

1.5

100

Penalty
coefficient R
108

Maximum
iterations
50

REFERENCES

TABLE V. RESULTS FOR THE ILLUSTRATIONAL PROBLEM

COMPARED WITH OTHER METHOD

number number of
of depots customers
10
100
10
150
10
180
10
200
10
220
20
100
20
150
20
180
20
200
20
250
25
200
25
220
25
250
25
300
25
350

CONCLUSION

In this paper, a two-phase method based on particle

swarm optimization for Location-Routing problem is
proposed. In the second phase, the paper adopts novel
particle representation to correspond to the integrated
optimization problem in logistics, and modifies the
objective function with additive penalty function part in
order to get the feasible solution. From an illustrational
example, the availability and performance of proposed
method can be approved.
Further two related research directions are as follows:
(1)develop a model and methods for the multi- product
multi-depot location-routing problem, and (2) develop the
hybrid PSO method with other possible intelligent
optimization algorithm to solute more complex locationrouting problem.

[1]Tai-His Wu, Chinyao Low, Jiunn-Wei Bai. Heuristic

solution to multi-facility location-routing problems.
Computers &operations Research.2002;29:1393-1415.
[2]Webb MHJ. Cost functions in the location of facility for
multiple delivery journal Research Quarterly. 1968; 19:31120.
[3]Salhi S,Rand GK. The effect of ignoring routes when
locating facilitys. European Journal of Operational
Research.1989;19:311-20.
[4]Nagy G,Salhi S. Nested heuristic methods for the
location-routing problem. Journal of operational Research
Society.1996;47:1166-74.
[5]Nagy G, Salhi S. A nested location-routing heuristic
using route length estimation. Studies in Locational
Analysis. 1996;10:109-27.
[6]Chen TW. Heuristic procedures for practical-sized
uncapacitated location-capcitated routing problem. Decision
Sciences. 1993;24:995-1021.
[7]Min H,Jayaraman V,Srivastava R. Combined locationrouting problem:a synthesis and future research directions.
European Journal of Operational Research. 1998;108:1-15.
[8]Maria Albareda-Sambola, Juan A, Diaz, Elena
Fernandez. A compact model and tight bounds for a
combined location-routing problem. Computer &
Operations Research. 2005;32:407-428.

Tabu Search
Based PSO
cost cpu
cost
cpu
79542
8 79238
9
121564
10 110256
10
156325
11 152357
11
162321
13 160241
14
173243
14 172314
14
67356
7 67135
7
102325
9 110234
8
121325
11 121034
10
134241
12 128967
12
142315
14 147782
15
124531
14 122215
14
131025
15 130021
16
150124
15 145781
15
162437
17 161325
17
180219
19 179324
18

Table shows the average solutions and average CPU

times for this proposed two-phase PSO based method and

244

[9]S.C.Liu, S.B.Lee. A two-phase heuristic method for the

muti-facility location routing problem taking inventory
control decisions into consideration. Int J Adv Manuf
Technol. 2003;22:941-950.
[10]Kennedy J, Eberhart RC. Particle swarm optimization.
Proceedings of the IEEE International Conference on Neural
Networks.(1995)1942-1948.
[11]Eberhart RC.Shi Y. Particle swarm optimization:
developments, applications and resources. Proc .Congress
on evolutionary computation, 2001, Seoul, Korea.P.81-6.
[12]Ayed Salman,Imtiaz Ahmad,Sabah Al-Madani. Partical
swarm optimizaion for task assignment problem.
Microprocessors and Microsystems. 2002,26:363-371

245