Problem
OSMAN ALP and ERHAN ERKUT
School of Business, University of Alberta, Edmonton, Alberta, T6G 2R6, Canada
ZVI DREZNER
Department of Management Science/Information Systems, California State University, Fullerton,
CA 92634-6848, USA
Abstract. We propose a new genetic algorithm for a well-known facility location problem. The algorithm
is relatively simple and it generates good solutions quickly. Evolution is facilitated by a greedy heuristic.
Computational tests with a total of 80 problems from four different sources with 100 to 1,000 nodes indicate
that the best solution generated by the algorithm is within 0.1% of the optimum for 85% of the problems.
The coding effort and the computational effort required are minimal, making the algorithm a good choice
for practical applications requiring quick solutions, or for upper-bound generation to speed up optimal
algorithms.
Keywords: facility location, p-median, genetic algorithm, heuristic
1. Introduction
In this paper we propose a new genetic algorithm for the p-median problem, which is
arguably the most popular model in the facility location literature. The goal of the model
is to select the locations of p facilities to serve n demand points so as to minimize the
total travel between the facilities and the demand points. This is a combinatorial opti-
mization problem shown to be NP-hard by Cornuejols, Fisher and Nemhauser (1977).
Many researchers devised heuristic methods to solve large instances of this problem to
near optimality with reasonable computational effort.
Genetic algorithms (GAs) are heuristic search methods that are designed to mimic
the evolution process. New solutions are produced from old solutions in ways that are
reminiscent of the interaction of genes. GAs have been applied with success to problems
with very complex objective functions. While a number of applications to combinatorial
optimization problems have been reported in the literature, there are few applications of
genetic algorithms to facility location problems. The GA we describe in this paper has
a number of desirable properties: it is simple, it generates excellent solutions, and it is
fast.
In section 2, we introduce the p-median problem formally and briefly describe the
relevant literature. In section 3, we discuss the properties of GAs in general terms and
review the cross-section of the p-median and GA literatures. We describe our GA in
22 ALP, DREZNER AND ERKUT
section 4, and provide a small numerical example in section 5. Section 6 contains the
results of our computational experience with the GA.
2. The p-median problem
The p-median model is a location/allocation model, which locates p facilities among n
demand points and allocates the demand points to the facilities. The objective is to min-
imize the total demand-weighted distance between the demand points and the facilities.
The following formulation of the p-median problem is due to ReVelle and Swain (1970).
min
n
summationdisplay
i=1
n
summationdisplay
j=1
w
i
d
ij
x
ij
s.t.
n
summationdisplay
j=1
x
ij
= 1 ∀i,
x
ij
lessorequalslant y
j
∀i,j,
n
summationdisplay
j=1
y
j
= p,
x
ij
= 0or1 ∀i,j,
y
j
= 0or1 ∀j,
where
n=total number of demand points,
x
ij
=
braceleftbigg
1 if point i is assigned to facility located at point j,
0 otherwise,
y
j
=
braceleftbigg
1 if a facility is located at point j,
0 otherwise,
w
i
=demand at point i,
d
ij
=travel distance between points i and j,
p=number of facilities to be located.
This is an uncapacitated facility location model where every demand point is served
by one facility and trips to demand points are not combined. It is a useful strategic
planning tool for many single-level distribution systems (for application examples see
(Fitzsimmons and Austin, 1983; Erkut, Myroon and Strangway, 2000)).
Teitz and Bart (1968) presented one of the oldest and most popular heuristic algo-
rithms for this problem. It is essentially an exchange heuristic that starts with a random
solution and improves it iteratively by swapping facilities in and out of the solution.
If one uses multiple starts, this method generates very good solutions when applied
to smaller problems. Mathematical programming-based heuristics were suggested by
EFFICIENT GENETIC ALGORITHM 23
Narula, Ogbu and Samuelsson (1977), which combines Lagrangian relaxation with sub-
gradient optimization, and Galvão (1980), which utilizes the dual of the problem. More
recently, Densham and Rushton (1992) proposed a two-phase search heuristic, and Piz-
zolato (1994) suggested a decomposition heuristic that works on a forest of p trees.
Modern heuristics have been applied to the p-median problem as well. A Tabu
search algorithm was designed by Rolland, Schilling and Current (1997). Simulated
annealing algorithms were designed by Murray and Church (1996) and Chiyoshi and
Galvão (2000). Rosing, ReVelle and Schilling (1999) proposed a composite heuristic
which found optimal solutions with regularity. We discuss GAs for the p-median prob-
lem in the next section.
Despite the combinatorial nature of the problem, optimal algorithms can solve
some instances of the p-median problem with reasonable effort. A typical approach
is to combine Lagrangian relaxation with branch-and-bound (see, for example, (Daskin,
1995)). Galvão and Raggi (1989) complemented this approach with a primal–dual al-
gorithm. Similarly, Koerkel (1989) combined a primal–dual algorithm with branch-and-
bound. In contrast, Avella and Sassano (2001) presented two classes of inequalities and
use them in a cutting plane algorithm.
3. Genetic algorithms
GAs are computational solution procedures that have been inspired by biological pro-
gression. They are intensive search heuristics that allow solutions to evolve iteratively
into good ones. They were first proposed as problem-solving methods in the 1960s, but
they have become popular in the operations research literature more recently. Thorough
treatments of GAs can be found in (Reeves, 1993; Dowsland, 1996).
The chromosomes in a GA correspond to solutions in an optimization problem.
The one-to-one relationship between a chromosome and the corresponding solution is
determined by an appropriate encoding. There is a fitness function that evaluates the
quality of a chromosome. Pairs of chromosomes are selected by the fitness function
and crossed to produce new chromosomes – this is the primary mechanism by which
solutions evolve. Mutations are used to promote genetic diversity (hence to facilitate
a thorough search of the feasible region). GAs work well for complex optimization
problems since they preserve the common sections of the chromosomes that have high
fitness values. They consistently disregard poor solutions and evaluate more and more
of the better solutions.
Although the general principles of GAs are common to all GA applications, there
is no generic GA, and the user has to custom-design the algorithm for the problem at
hand. This is not a trivial task since a GA requires many design decisions. The encoding
is critical since a poor choice may result in a poor algorithm regardless of its other
features. Other important decisions include the size of the population in one generation
(i.e., the number of solutions), the selection of parents (selection of two old solutions to
produce new solutions), the crossover operator (the method by which new solutions are
24 ALP, DREZNER AND ERKUT
produced by old solutions), the replacement of one generation by the next, the method
and frequency of the mutations, and the number of generations.
In principle, GAs can be applied to any optimization problem. Given that GAs
have been applied to many combinatorial optimization problems with success (Reeves,
1993), we expect them to work well on location/allocation problems. Yet little effort
has been directed towards designing GAs for location problems in general, and for the
p-median problem in particular. In the first paper applying the GA framework to the
p-median problem, Hosage and Goodchild (1986) encoded the solutions of the problem
as a string of n binary digits (genes). This encoding does not guarantee the selection of
exactly p facilities in each solution, and the authors use a penalty function to impose this
constraint. Primarily due to this encoding choice the algorithm performs rather poorly
even on very small problems.
Dibble and Densham (1993) describe another GA for a multi-criteria facility loca-
tion problem, solving a problem with n = 150 and p = 9 with a population size of 1000
and 150 generations. They report that their algorithm finds solutions that are almost as
good as the solutions generated by an exchange algorithm, but with more computational
effort. Moreno-Perez, Moreno-Vega and Mladenovic (1994) design a parallelized GA
for the p-median problem, where multiple population groups exist and individuals are
exchanged between these groups. This feature has effects similar to mutation since it
prevents premature homogenization in the population groups. They do not report com-
putational results.
In the most recent application of GA to the p-median problem, Bozkaya, Zhang
and Erkut (2002) describe a fairly complex GA and report the results of an extensive
computational experiment with algorithm parameters. This algorithm is able to produce
solutions that are better than the solutions of an exchange algorithm. However, con-
vergence is very slow. In contrast, the algorithm described in the next section is much
simpler and produces good solutions very quickly.
4. A new genetic algorithm
4.1. Encoding
We use a simple encoding where the genes of a chromosome correspond to the indices
of the selected facilities. For example, (5,7,2,12) is a chromosome that corresponds
to a feasible solution for a 4-median problem where demand points 2, 5, 7, and 12
are selected as facility locations. This encoding ensures that the last constraint in the
formulation is always satisfied.
4.2. Fitness function
The fitness function is the objective function of the p-median problem. The fitness of
a chromosome is identical to the objective function value of the solution it corresponds
to, and it can be calculated easily using the problem data. The calculations assume
EFFICIENT GENETIC ALGORITHM 25
that every demand point would be allocated to one facility, namely to the closest open
facility. This ensures that the first two constraints in the formulation are always satisfied.
Hence, the selections of the fitness function and the encoding satisfy all constraints and
no additional effort is needed in the implementation of the algorithm to enforce the
constraint set.
4.3. Population size
The GA works on a population with fixed size. Large populations slow down the GA
while small populations may not have sufficient genetic diversity to allow for a thorough
search over the feasible region. Hence, the selection of the population size is important.
We target population sizes with the following two properties:
1. Every gene must be present in the initial population. Since the GA iterates by creating
new combinations of the genes in the original gene pool, an incomplete gene pool will
result in a partial search over the feasible region, and the GA would have to rely on
mutations to bring missing genes into the pool. The minimum number of members
to represent each gene in the initial population is ceilingleftn/pceilingright, the smallest integer greater-
than-or-equal to n/p.
2. The population size should be proportional to the number of solutions. In general, the
larger the feasible region of a problem, the more difficult it is to find the best solution.
Our goal is to devise a formula to generate population sizes with these properties
for different problem parameters. Let S = C
According to this formula, the population size is an integer multiple of d, and this satis-
fies the first property. In fact, the result of the max operator is at least two, guaranteeing
that every gene appears at least twice in the initial population. (If each gene appears
only once in the initial population, it might disappear – once and for all – in an early
iteration which may make it difficult to find optimal solutions.) The ln(S) term provides
the increase in the population size in proportion to S, the number of solutions. Yet the
increase is very slow (due to the ln operator), which keeps the population size manage-
able even for very large problems such as n = 1000 and p = 100. Figures 1–3 show
how this population formula size behaves as a function of n and p.
There are many other formulas that have the two properties listed above. The one
we suggest is merely one that works. It generates very good solutions over a large
spectrum of problem parameters in our empirical tests. It may be possible to find a
formula that generates better solutions or works faster – we do not claim that our formula
is “optimal” in some sense.
26 ALP, DREZNER AND ERKUT
Figure 1. Population size as a function of n for three different values of p.
Figure 2. Population size as a function of p for three different values of n.
4.4. Initializing the population
As discussed above, every gene should be present in the initial pool. It is also desirable
to have each gene approximately with the same frequency in the initial pool so the gene
pool is not biased. Suppose that the population size is equal to kd, for some constant k.
For the first set of n/p members, we assign the genes 1,2,...,p to the first member,
the genes p + 1,p+ 2,...,2p to the second member, and so on. For the second set
of n/p members, we distribute the genes similarly, but we use an increment of two in
EFFICIENT GENETIC ALGORITHM 27
Figure 3. Population size as a function of n for problems with n = 5p.
the sequences. For example, we assign the genes 1,3,5,...,2p−1tothefirstmember
in the second group. Similarly, in the kth group we distribute the genes to the members
sequentially with an increment of k. For example, for (n,p,k)= (12,4,2), we would
have the following initial population: (1,2,3,4), (5,6,7,8), (9,10,11,12), (1,3,5,7),
(9,11,2,4), (6,8,10,12).Ifn/p is an integer then each gene is represented in the initial
population with an equal frequency. If n/p is not an integer then after distributing all of
the genes from 1 to n to each group, we allocate random genes to fill empty slots.
4.5. Selecting the parents
The parents are selected randomly from the population. We experimented with biased
selection mechanisms where fitter parents are more likely to be selected for reproduc-
tion, but experienced no significant impact on the performance of the algorithm. While
both versions of the algorithm produced excellent results, we obtained marginally better
results with randomly selected parents.
4.6. Generating new members
In a typical GA two parents are selected for reproduction, and their chromosomes are
merged in a prescribed way to produce two children. Usually the chromosomes of the
parents are split into two, creating four partial chromosomes, and then these four pieces
are combined to create two new chromosomes. For example, the parents (1,2,3,4,5)
and (6,7,8,9,10) would create the children (1,2,3,9,10) and (6,7,8,4,5) if a
crossover after the third gene is used. We do not use such a traditional crossover op-
erator. Instead, we take the union of the genes of the parents, obtaining an infeasible
28 ALP, DREZNER AND ERKUT
solution with m genes where m>p. Then, we apply a greedy deletion heuristic to
decrease the number of genes in this solution one at a time until we reach p.However,
we never drop genes that are present in both parents. To reduce the number of genes by
one we discard the gene whose discarding produces the best fitness function value (i.e.,
increases the fitness value by the least amount). We call the infeasible solution obtained
after the union operation the “draft member” and the feasible solution generated by the
heuristic as the “candidate member”.
The generation operator can be summarized as follows.
Generation Operator.
Input: Two different members.
Step 1. Take the union of the input members’ genes to obtain a draft member.
Step 2. Let the total number of genes in this draft member be m. Call the genes that are
present in both parents fixed genes and the rest free genes.
Step 3. Compute the fitness value of this draft member.
Step 4. Find the free gene that produces the minimum increase in the current fitness
value when deleted from the draft member, and delete it. Repeat this step until
m = p. Let this final solution be a candidate member.
Output: A candidate member.
Our children generation method is distinctly different from a crossover operator.
We use a greedy heuristic instead of blindly creating children by cut-and-paste of chro-
mosome parts. This increases time demands, but it also improves the quality of the
children. While this can be viewed as compromising a major strength of the GA, our
computational experience shows that this operator works well. Berman, Drezner and
Wesolowsky (2002) report success with a similar merge-drop operator in solving an-
other location problem.
4.7. Mutation
The mutation operator is an important component of a GA. This operator helps the algo-
rithm to escape from local optima. Although we experimented with different mutation
techniques (for example, in 2% of new member generation iterations we added a facility
not found in the union of the two parents before starting the generation operator, and
did not allow the deletion of this facility), they did not improve the performance of the
algorithm. Hence we decided not to use the mutation operator.
4.8. Replacement
We admit a candidate member into the population, if it is distinct (i.e., not identical to
an existing member), and if its fitness value is better than the worst fitness value in the
EFFICIENT GENETIC ALGORITHM 29
population, by discarding the worst member. Such a policy improves the average fitness
value of the population gradually while maintaining genetic diversity. We store the worst
and best members of the population after every population update.
The steps for the replacement operator are as follows.
Replacement Operator.
Input: One candidate member.
Step 1. If fitness value of the input candidate member is higher than the maximum fitness
value in the population, then discard this candidate member and terminate this
operator.
Step 2. If the candidate member is identical to an existing member of the current popu-
lation, then discard this candidate member and terminate this operator.
Step 3. Replace the worst member of the population with the input candidate member.
Step 4. Update the worst member of the population.
Step 5. Update the best member of the population.
Output: Population after introducing the candidate member.
4.9. Termination
After experimenting with different iteration counts that depend on the problem parame-
ters, we opted for a convergence-based stopping criterion. The algorithm terminates after
observing ceilingleftn
√
pceilingright successive iterations where the best solution found has not changed.
One iteration consists of one use of the generation and replacement operators. For exam-
ple, for a problem with (n,p) = (100,20), the algorithm terminates if 448 successive
children fail in improving the best solution. (For all of our test problems, we had n>2p.
For a problem with n lessorequalslant 2p, we would suggest stopping after n
√
(n−p) non-improving
iterations.) An alternative termination rule (stop after a given number of successive it-
erations where the entire population has not changed) generated better solutions, but
required more computational effort.
4.10. Algorithm
The overall algorithm can be stated as follows.
Algorithm.
1. Generate an initial population of size P(n,p)as described in section 4.4.
2. Initialize a variable for keeping track of successive iterations where the best solution
found has not changed, MaxIter.
3. Repeat while MaxIter lessorequalslantceilingleftn
√
pceilingright:
30 ALP, DREZNER AND ERKUT
3.1 Randomly select two members from the current population.
3.2 Run the Generation Operator: input these two members and obtain a candidate
member.
3.3 Run the Replacement Operator: input candidate member if possible.
3.4 If the best solution found so far has not changed, then increment MaxIter.
4. Select the best member as the final solution of the algorithm.
5. Numerical example
We use a very small problem with n = 12 and p = 3 to provide a numerical example.
The coordinates of the 12 demand points are given in table 1, and the Euclidean metric
is used to measure distances. For P = 8 the initial population along with their fitness
values is shown in table 2.
In the first iteration, members 1 and 4 are combined. The resulting draft member
is 1-2-3-10-11-12 with a fitness value of 136. Three genes have to be dropped from this
member to produce a feasible solution. Dropping 12 results in the smallest increase in
the fitness value (to 156). Then we drop 2 and 1, and generate the candidate member
3-10-11 with fitness value 241. This is not a current member of the population and
its fitness value is lower than that of the worst. Hence, we delete 2-4-6 and introduce
3-10-11 as member 7 of the population in its place.
In iteration 2, members 4 and 7 are selected for merger. The resulting draft member
is 3-10-11-12 with a fitness value of 210. Genes 10 and 11 are fixed since they appear
in both parents. The greedy drop heuristic selects gene 12 for deletion, and the resulting
Table 1
The coordinates of the 12 demand points (unweighted).
Point12345678910112
x 2 2292 267983325953631
y 55 91 91 99 70 99 52 6 44 88 36 71
Table 2
The initial population.
No. Member Fitness
1 1-2-3 352
2 4-5-6 316
3 7-8-9 348
4 10-11-12 257 (best)
5 1-3-5 358
6 7-9-11 365
7 2-4-6 391 (worst)
8 8-10-12 271
EFFICIENT GENETIC ALGORITHM 31
Table 3
Summary of the first 10 iterations of GA.
Iter. Merge Draft Candidate Fitness Replace Comment
1 1 and 4 1-2-3-10-11-12 3-10-11 241 7 new best
2 4 and 7 3-10-11-12 3-10-11 241 – identical to 7
3 8 and 2 4-5-6-8-10-12 5-8-10 266 6
4 1 and 2 1-2-3-4-5-6 1-3-6 282 5
5 2 and 7 3-4-5-6-10-11 3-10-11 241 – identical to 7
6 3 and 5 1-3-6-7-8-9 3-7-9 245 1
7 1 and 8 3-7-8-9-10-12 3-9-10 236 3 new best
8 4 and 5 1-3-6-10-11-12 3-10-11 241 – identical to 7
9 2 and 6 4-5-6-8-10 5-8-10 266 – identical to 6
10 1 and 2 3-4-5-6-7-9 2-6-9 262 2
member is 3-10-12. This is identical to the current member 7, and the population remains
unchanged. In iteration 3, members 2 and 8 are merged to create 4-5-6-8-10-12. From
this draft, the greedy heuristic generates candidate member 5-8-10 with fitness value 266,
which replaces member 6 in table 2. The first 10 iterations are summarized in table 3.
The optimal solution is 3-9-10, which is found in iteration 7.
6. Computational study
6.1. Accuracy and speed of the GA
To test the performance of our algorithm (abbreviated as ADE here), we solved four sets
of test problems:
1. OR Library: There are 40 p-median problems with known optimal solutions in the
OR Library (Beasley, 1990). The problem sizes range from n = 100 to 900 and
p = 5 to 200.
2. Alberta: We generated a 316-node network using all population centers in Alberta.
We computed distances using shortest paths on the actual road network and used
populations for weights. We solved 11 different problems on this network with values
of p ranging from 5 to 100, and we found the optimal solutions using Daskin’s (1995)
SITATION software.
3. Galvão: We solved 16 random problems from (Galvão and ReVelle, 1996) with n =
100 and n = 150 and various values of p.
4. Koerkel: We solved 13 problems from (Koerkel, 1989) for p values ranging from 2
to 333 on a network with 1,000 nodes.
(The Alberta, Galvão, and Koerkel problem sets are available through the Internet:
www.bus.ualberta.ca/eerkut/testproblems)
32 ALP, DREZNER AND ERKUT
In sum, we solved 80 problems of different sizes (ranging from 100 to 1,000 nodes)
from different sources, including random as well as realistic problems. The computa-
tional results for 10 replications of ADE, produced by a C++ code on a Pentium III
733 MHz computer with 128 MB memory, are summarized in table 4. In 28 of the
40 problems, ADE finds an optimal solution. The average deviation between the best
(worst) solution value found in 10 replications and the optimal objective function value
is 0.036% (0.291%). Another measure of ADE’s performance is the average deviation
between the average of the best solution values and the optimal values in 10 replica-
tions: 0.11%. In other words, the average of the best solutions over the 10 replications
is 100.11% of the optimal. These near-optimal solutions are generated with little com-
putational effort. ADE is very fast with an average per-replication time of 18.4 sec-
onds for the 40 test problems. The per-replication times range from 0.1 seconds for
(n,p) = (100,5) to 2.2 minutes for (n,p) = (900,90).
Table 5 summarizes the results for the other three sets of test problems. All Al-
berta problems and approximately half of the Galvão and Koerkel problems are solved
optimally. The average gap between the best solution found and the optimum is 0.05%
for the Galvão problems and 0.10% for the Koerkel problems. The per-replication time
spent on the Alberta and Galvão problems are 3.2 seconds and 0.5 seconds, respectively.
The 1,000-node Koerkel problems take longer: an average of 3 minutes per replication.
These times are similar to those spent on OR Library problems of comparable sizes.
Based on our experience with the 80 test problems, it is fair to say that ADE is
quite fast and find excellent solutions. The best solution found is optimal for more than
half of the problems. The average gap between the best solution and the optimum for
all the 80 problems is 0.045%, and the worst gap is 0.4%. Most of the replications take
under a minute and the longest one takes 3 minutes.
6.2. Comparison with other heuristics
In this section we first compare ADE using the OR Library problems with three other
heuristics: the genetic algorithm (BZE) by Bozkaya, Zhang and Erkut (2002), the
Gamma heuristic (Gamma) suggested by Rosing, ReVelle and Schilling (1999), the sim-
ulated annealing (SA) algorithm by Chiyoshi and Galvão (2000). In addition we compare
ADE with the three heuristics in the SITATION code (myopic, exchange, neighbour-
hood) using the Alberta and Galvão problem sets. (The 1,000-node Koerkel problems
are too large for the SITATION code.)
Although BZE and ADE are both GAs, there are a number of differences between
them. BZE uses the same encoding as ADE. Its initial population is generated at random
while ADE uses a systematic way to seed the initial population. BZE favors parents with
better fitness values whereas ADE selects parents at random. BZE uses three different
crossover operators while ADE uses a merge-drop heuristic to generate a child. BZE
uses mutations and invasions (a more intense form. of mutation) to maintain genetic
diversity, but ADE has no such features. BZE replaces the entire population with the next
generation while keeping the best solution in the population. ADE introduces children
EFFICIENT GENETIC ALGORITHM 33