Evaluating factorial experiments with simulated annealing for territorial optimization

18 Pages·2012·0.94 MB·English

Journal of Applied Operational Research (2012) 4(2), 91–108
© Tadbir Operational Research Group Ltd. All rights reserved. www.tadbir.ca

ISSN 1735-8523 (Print), ISSN 1927-0089 (Online)
Evaluating factorial experiments with
simulated annealing for territorial
optimization
María Beatriz Bernábe Loranca 1 and David Pinto Avendaño 2,*
Benemérita Universidad Autónoma de Puebla, Puebla, México
Abstract. A particular method for solving the problem of territorial optimization requires a classification process based on
clustering in which multiple comparisons must be performed in order to fulfill an objective function that minimizes the distances
among objects with the purpose of achieving geographical compactness. The computational complexity of this problem is
known to be in the NP-complete category, therefore, we have tackle the problem by means of heuristic methods. We have
selected the simulated annealing technique to be applied inside the clustering algorithm, because of its capability of finding
high quality suboptimal solutions. The clustering algorithm proposed considers both the properties of the classical partitioning
and main aspects of dynamic cloud methods. In general, in this paper we present a mathematical model and a geographical
partitioning algorithm for solving aggregation, a particular feature that exists in every problem of the territorial kind. The
algorithm is combinatorial and, therefore, it obtains approximate solutions throughout the execution of the simulated annealing
algorithm. The experiment carried out needed a factorial statistical analysis in order to model the parameters of the employed
heuristic. In this regard, we have constructed and analyzed the results of various factorial experiments Box Behnken and 2k
that have been applied to different instances obtained from the geographical partitioning algorithm which is implicit in the
problem of territorial optimization. The aim of these experiments is to identify the appropriate combination of parameters
for this particular problem. The dataset are of the census type, with a well defined spatial component and a vector of descriptive
variables.
Keywords: clustering; experiment; heuristic; territorial optimization
* Received October 2011. Accepted February 2012
Introduction
The territorial optimization (TO) problem consists in the classification of geographic units (GU) subjected to the
fulfillment of certain criteria, such as the geometric compactness which has been recently presented by (Bernábe
et al, 2009). This work describes a partitioning algorithm that fulfills the property of geometric compactness for
Geostatiscal Basic Areas (AGEBs) with the addition of Simulated Annealing (SA) to escape from local-optimal
solutions. At the same time, the performance of the partitioning algorithm has been improved. This algorithm is
used in regional partitioning problems, in particular to solve some population-related problems. On the one hand,
the implemented clustering algorithm uses as objective function the minimization of distances among objects

* Correspondence: David Pinto Avendaño, Facultad de Ciencias
de la Computación, Benemérita Universidad Autónoma de Puebla,
14 Sur & Av. San Claudio, CU, Edif. 104A, 72570, Puebla, Puebla,
México. E-mail: [email protected]
92 Journal of Applied Operational Research Vol. 4, No. 2

with the purpose of achieving compactness between AGEBs. On the other hand, the application involved in this
work, approximate solutions for “population” problems, an open problem in Territorial Design.
In previous works, a partitioning algorithm with exponential computational complexity (Bernábe et al, 2008)
named PAM (Partitioning Around Medoids) proposed by (Kaufman and Rousseeuw, 1987) has been implemented.
The results obtained by PAM (optimum value for certain number of groups) have allowed doing a comparison
with respect to the solutions generated by the heuristic proposed in this paper. The data considered for the
classification task correspond to the 423 AGEBs of the Metropolitan Zone of the Valley of Toluca (MZVT), and
the classification variables correspond to 57 socio-economical variables available for those areas (INEGI, 2000).
The Simulated Annealing algorithm together with the partitioning algorithm was used in order to cluster AGEBs.
The centroids are chosen randomly in order to determine the number of groups (clusters). There exists a swap between
objects and centroids (AGEBs) until the cost function converges, in order to consider all the possible cases (allowed
by the heuristic parameters) in the classification process, finding a minimum for the implicit cost function to be
calculated. Finally, in order to determine the quality of the solutions obtained by the proposed algorithm, we have
collected results of previous works which are then introduced in a number of factorial experiments presented in
this paper.
This document is organized as follows: Section 1 Introduction, settles down the purpose of this paper, Section 2
describes related works; Section 3 introduces the mathematical formulation of geographical units and the associated
algorithm; Section 4 presents the processing of the generated instances in a set of factorial experiments; Section 5
describes factorial experiment 2k for Simulated Annealing algorithm and finally, in Section 6 the conclusions are
given.
Related works
The Territorial Optimization (TO) problem falls within the Territorial Design (TD) category, a broad research
line that rise diverse regionalization problems. It is important to point out that there are some features in territorial
design that are common with TO, for instance, definitions, scope, mathematical models, solving methods, and result
validation, among other things. However, some differences can be found in the definitions of regional partitioning,
zoning, territorial aggregation, regional distribution, zone categorization, geographical clustering, territorial
rearrangement or regional design. As an example, the term alignment is often used instead of design, and aggregation
is used instead of classification.
In general terms, territory design can be seen as a problem of grouping small geographical areas (basic areas)
in larger geographical clusters called territories, in such way that an acceptable grouping is the one that meets
certain predetermined criteria according to the problem at hand (Zoltners and Sinha, 1983). Depending on the
context, these criteria may be economically motivated (average potential sales, jobs or number of salesmen) or
may have a demographic background (number of inhabitants, voting population).
The criteria or properties to meet in TD problems depend on each problem, with spatial restrictions as continuity
and geometric compactness being in high demand. In this sense, these problems tend to begin with a description
in terms of TD with an optimization approach subjected to certain criteria that allow mathematical modeling with
the presence of a cost function together with the properties of the problem that are formalized as restrictions.
Additionally, the corresponding implementation has to be aided by heuristic methods in order to obtain an approximate
solution close to the optimum.
There are international efforts aimed to automatically generate geographical groupings. From the seminal
work of (Garfinkel and Nemhause, 1970), to the more recent of (Ricca and Simeone, 1997), (Bozkaya et al.,
1999), (Bozkaya et al., 2003), (D’Amico et al, 2002), (Romero, 2002) among others, have been directed mainly
towards solving political districting. In the case of homogeneous territorial design problem, the work of (Duque,
2004) stands out; and, as far as we know, no one has approached territorial aggregation with a combinatorial
optimization method as an aid in the generation of groups considering the AGEBs as territorial units. In Mexico,
(Romero, 2002) is the main initiator of this research line for parceling of the national territory. Romero takes city
blocks as the geographical units to group (which makes easier to establish the compactness and contiguity by resorting
to computational geometry); however, the known adjacency methods to obtain geometrical compactness do not
help when such units are separated by non-uniform distances, and this is precisely the spatial nature of the
AGEBs in Mexico.
MBB Loranca and DP Avendaño 93

The problem of TO stands out as a branch of Territorial Design (TD). Problems such as point of sale assignment,
location of services either population-related or census-related and electoral districting are approached from the
TD standpoint. The classification of geographical units (term also known as Territorial Aggregation (TA) is a
complex process:
a) There exists a large quantity of data of spatial nature to classify, which implies that the classification process
is slow in terms of computing units b) since the classification must consider mainly the characteristics and restrictions
of the problem, those demand to be interpreted and modeled mathematically as similitude and dissimilitude
measurements subjected to an optimization model (which in turn determines the fundamental part of the geographic
cluster algorithm). However, such model reveals the computational complexity of the problem, confirmed during
the implementation of the algorithm, making necessary the inclusion of heuristic methods.
The mathematical optimization models formulated in this way are always supported by the inherent properties
that define the TD problem: compactness, homogeneity, connectedness and contiguity.
Territorial design can be seen as a problem of grouping small geographical areas (basic areas) into larger
geographical clusters, in such a way that an acceptable grouping is the one that meets certain predetermined criteria
for the problem at hand. Depending on the context, these criteria may be of economic motivation (average potential
sales, jobs, number of salesmen) or it may have a demographic backing (number of inhabitants, voting population).
Spatial restrictions such as continuity and compactness are often demanded (Kalcsics et al, 2005).
The problem of geographic zone design or planning (also known as generalized districting or a particular case
of it) occurs when n units of area are aggregated into k zones in such a way that a value function is optimized,
depending on the restrictions on the zones’ topology, e. g. internal connectivity (Bação et al, 2004).
It is important to emphasize that the territorial planning problems require a territory aggregation model used to
build “small” groups of zones to ease both the analysis and solution of the problem (Bação et al, 2005).
In this section a brief summary has been presented that highlights the implications of solving problems of territory
aggregation, the most demanded task is a solution to the zone classification process based on restrictions, where
the grouping demands a large number of operation, a NP-Hard category according to Bação et al. (2004).
As might be expected, this demands a combinatorial optimization model supported by a meta-heuristic to reach
a grouping that is a good solution. This problem is recognized as a special case of the knapsack problem or the
clustering problem (Bação et al, 2004), and it is in this direction that the present work fits, resorting to simulated
annealing as a solving method.
Recently Duque has produced an excellent taxonomy of the geographical aggregation problems, where several
references to the application of simulated annealing in this context can be found (Duque et al, 2007). As we have
already mentioned, the problem of Territorial Optimization requires finding a correct clustering of data based on
their geometric compactness, which is understood as the process of clustering objects by using the minimization
of distances among objects as the optimization criterion. Thus, we have studied the properties of partitioning
algorithms in order to solve geographical compactness in OP. In (Kaufman, 1987) it is described an exhaustive
proposal of partitioning around medoids; this algorithm has been implemented in previous research works with
the purpose of executing initial tests for finding an exact value which may be used as a comparison reference
with respect to the optimum value in further implementations (Bernabe et al, 2006). However, the PAM (Partitioning
Around Medoids) algorithm is exhaustive and, therefore, its complexity time may be prohibitive for being used in
TO. Alternatively, in Piza et al. (1999) it has been presented a taxonomy of the different classification techniques
that have been employed in the aforementioned algorithm of partitioning in TO.
Mathematical formulation
The following considerations are introduced before presenting the model, so that the reader may be aware of these
basic and useful restrictions of common use in the task of modeling of aggregation for territorial problems.
Let the initial set of area units be X=x ,x ,...,x  where x is the ith area unit. The number of zones is
1 2 n i
designated by k and Z is the set of all area units that belong to the zone Z then
j i
Z for i = 1,...,k
i
Z Z  for i  j
i j
94 Journal of Applied Operational Research Vol. 4, No. 2

k Z X
i1 i
These considerations constitute the set of constraints that can be applied equally in clustering and in zone design
(Bação et al, 2005). On the other hand, it is clear that these restrictions correspond to the properties of the classical
partitioning, and, we must remember that we are interested for a unique object partition with k disjoint classes in
classical partitioning. The traditional theory of partitioning methods include the following ones: k-means of
Forgy, dynamic clouds of Diday, algorithm of transfers of Régnier, among others (Piza et al, 1999).
Let us consider a set of finite objects (n) to be classified x ,x ,...,x  with k < n the number of target classes.
1 2 n
In the classification by partitioning task, a partition PC ,C ,...,C  of Ω in k classes C ,C ,...,C is
1 2 k 1 2 k
characterized by the following two conditions:
k
1) C 2)C C  with i j
i j
i
i1
In the different problems of TD, the compactness and contiguity is described according to the problem to
solve, some topological definitions are taken as basis and then adapted to the problem. The spatial nature of the
geographical data to group determines the way that these points are defined and on the other hand how they are
modeled.
This work will not present the mathematical details of geometric compactness, however, it is informally understood
as compactness the fact that objects (geographic units) are close one another, for this reason the implicit objective
function of the geographical clustering algorithm minimizes distances between objects. In the problem at hand,
the AGEBs, are composed of a set of 55 city blocks in average and separated by different distances, so
that the adjacency is not adequate to formulate this problem.

Problem formulation
Let T  R2be a territory made up of AGEBs (denoted by A), i.e.,
j
 
T  A ,A ,...,A ,...,A such as
1 2 j n
T  n A and A A i  j.
j1 j i j
Let C = {c , c ,…, c ,…,c } be the centroid set for each AGEB of T, where c = (x , y). We want to obtain a set
1 2 j n j i j
P, with p > 1, of AGEBs, P G ,G ,...,G ,...,G  such that P be a partition of T with p < n and each
1 2 k p
G is a grouping of AGEBs. Under the geographical compactness criterion, the objective function consists of
k
minimize the Euclidean distance of one of its p centroids c with respect to each other centroids c of the same
j i
group G, and the solution space Ω is the set of all the partitions of T with cardinality p.
i
n (1)
min  d(c ,c )k 1,2,...p
j i
P j1ciGk
Equation 1 corresponds to the objective function of the geographical clustering problem, which considers the
properties of the classical partitioning. Note that in the case of an exhaustive search in Ω, we may calculate the
number of alternative solutions that must be evaluated, by the following equation:
n! (2)
 n 
p p!n p!
where n is the number of AGEBs, and p is the number of expected groups. Equation 2 indicates that the
complexity algorithm is in O(np) for p << n, thus, polynomial in n for a given value of p, but the number of
combinations is huge. For instance, if n = 50 and p = 10, then
50 1027227819 1010.
10
In other words, if we evaluate 1,000,000 combinations per second, then the algorithm will take at least 3 hours
for enumerating all the possible combination. But if n = 100 and p = 15, the algorithm would take more than
a millennium for evaluating all the possibilities. Thus, for variable values of p, the problem still being NP-
complete because the required time for solving is by exhaustive enumeration is:
MBB Loranca and DP Avendaño 95

n (3)
n   2n 1O(2n)
p
p1
The problem is exponential in n as may be deduced from Equation 3, which justifies the use of meta-heuristics
for finding solutions close to the optimum in a reasonable time of execution. The solutions to the territorial
partitioning problem are subset of C = {c | j=1, 2, …,n} of the type {c , c , ..., c } with p < n, where p is the number of
j 1 2 p
partitioning groups, n is the number of territorial AGEBs, and c are the centroids of each AGEB. The implementation of
j
the algorithm for constructing the groups considers the following restrictions:
X 1 i,e X 1 (4)
ij 0 ij
if the AGEB i belong to group j, and its 0 otherwise
n (5)
X 1j1,2,..n
ij
i1
n (6)
X  pi 1,2,..n
ij
i1
Equation 4 is decision variable for determining wheter or not an AGEB is included in a given group. Equation
5 guarantees that an AGEB may belong exclusively to one group. Equation 6 determines the exact number of
groups. Solutions require the conformation of p groups. Each group is made up of AGEBs with the closest centroids
to each one of the p centroids.
The initial solution is generated randomly and represented as a vector. Thus, let be x = (c , c , ..., c ) an initial
0 1 2 p
solution, in which each component is a centroid of an AGEB that we will call gravity center because each group
will be conformed around this centers.
SA heuristic to obtain local-optimal solutions
The SA algorithm is a neighbor-based search method characterized by an acceptance criterion for neighboring
solutions that adapts itself during running time. It uses a variable called temperature (T), which value determines
to what degree can neighboring solutions worse than the current one be accepted. The temperature variable is
initialized with a high value, called initial temperature To and is reduced with each iteration through a cooling
temperature mechanism alpha (α) until a final temperature Tf is reached. During each iteration a concrete number
of neighbors L(t) is generated, which may be fixed for all the execution time or may change for each iteration.
Each time a neighbor is generated, an acceptance criterion is applied to determine if it will substitute the current
solution.
To escape local optima solutions worse than previously seen can be accepted. The likelihood of movements to
worst solutions decreases as the search proceeds and, at the end, only good solutions are accepted. This is controlled
by a probability density function following the work of Metropolis in the field of statistical thermodynamics (Kirkpatrick
et al, 1983). Basically, he modelled a cooling process by simulating the energy changes in a system of particles with
decreasing temperature, until it converges to a stable state. The laws of thermodynamics say that at a temperature
T, the probability of an energy increase δE can be approximated by P [δE] = exp (- E / kt), where k is the Boltzmann
constant. In the model of Metropolis, a random perturbation is introduced in the system and the resulting energy changes
are calculated: if the energy drops, the change is accepted, otherwise the change is accepted with probability P [δE].
The objective function previously described (1) will be implemented with the following SA algorithm:

INPUT To, α, L(t), Tf
T To Initial value for the control parameter
Sact  Generate initial solution
WHILE T >= Tf DO Stopping condition
BEGIN
96 Journal of Applied Operational Research Vol. 4, No. 2

FOR cont  1 TO L(t) DO Cooling speed (T)
BEGIN
Scand  Select solution N(Sact) Creation of a new solution
 cost(Scand) - cost(Sact) Computation of cost difference
IF U(0,1) < e(/T) OR <0 Application of acceptance criterion
THEN Sact Scand
END
T(T) Cooling mechanism
END
{Write as solution the best of the visited Sact}

Different classification algorithms agree with the central idea of objective function (1) to minimize the distance
between the objects to their centroid. Along the same idea, a clustering algorithm has been built exclusively for
grouping geographical units with the addition of SA. The following pseudo code has been commented to highlight the
operation of both the loops and the calculation of the cost function.
Simulated annealing algorithm embedded in the geographical partitioning for territorial
optimization
The partitioning method of dynamic clouds has highly influenced the resolution of aggregation in TO. This method
looks for a partition of Ω in K classes well aggregated, separated and disjoint. The number K of classes is given
in advance disregarding the domain of the data. Dynamic clouds algorithm, introduced by (Diday, 1980) as a
generalization of the Forgy’s k-means clustering method, is based on the assumption that each class must have a
representative (named kernel) which are then used for conforming portioning iteratively until a certain criterion is
optimized. The algorithm is described as follows:
1. An initial partition of Ω is given
2. The kernels are computed by means of a representation function
3. A partition is constructed by assigning each element to the closest kernel by means of an assignation function
4. Steps 2 and 3 are repeated until there exist some degree of convergence in the classes obtained.
The initial kernels are usually randomly selected. In the general case, m elements are chosen K times from the
individual by using some additive criterion such as the one that follows:
k
W(P)=   D(x,N )
i k
k1xC
i k
where, N is the kernelof C (made up of m objects), and D is a measure of dissimilarity (for instance, an
k k
aggregation) between the objects x and the kernels N (which are also a set of objects). The kernel N is defined
i k k
as the subset of C with m elements minimizing
k

 D(x,N ).
iC i k
k

It may be proved that these formulae converge, improving W and consequently, the previous partitioning. It
may be seen that the method of dynamic clouds and the method of k-means are identical when the kernels are
gravity centers. In our particular problem, kernels and centroids are analogous.
The following algorithm written in pseudo-code has been constructed following the fundamentals of dynamic
clouds, adding simulated annealing as the optimization method:

Let n be the number of objects to classify Ug stands for the object i being assigned to the centroid j
ij
i = 1,..,n and j = 1,…,k
Let M = {M ,M ,,M } be a solution with K centroids
1 2 k
To initial temperature
Tf final temperature
L(t) Number of iterations that will be performed with the same temperature
MBB Loranca and DP Avendaño 97

Begin
TTo
‘Obtain the initial solution
Randomly generate initial centroids M = {M M , … , M }
1, 2 k
‘Any AGEB could be a randomly chosen to be centroid

cur_cost  Cost (M)*
‘This assignment already represents an initial solution; it is a proposed solution generated
in the previous step. The next steps will generate another solution to determine how good
it is in relation to the current one and to decide whether the current solution is changed
or not.

While T  Tf
‘while the system is not cold
For count=1 to L(t) do
‘number of loops to perform with the same temperature (SA parameter)
C Generate a random solution
‘a solution to compare with * is generated
cand_cost Cost (C)
‘the cost for the candidate solution just generated is calculated
  cand_cost – cur_cost
‘cost difference to obtain the likelihood of accepting the candidate solution
If U(0,1) < e/T or < 0 do
‘if the likelihood of accepting is still high
MC
‘if the candidate solution is accepted
cur_cost  cand_cost
End If
End For
T (T)
‘the system is cooling down
End While

2. Cost Function (Sol)
‘determine how good is the solution SOL, i. e. how much minimizes the objective

i  1
‘initialize first object
cost  0
While i  n
‘for each object in Ugdo
if Ug is not a centroid then
i
dmin  dist(Sol ,Ug )
1 i
‘represent the distance from object to Sol (first centroid where Sol represents the set
i 1
of all centroids. Calculate the distance from each object to its nearest centroid (distance
of an object that is not a centroid to Sol , that is the centroid 1)
i 1
j  2
‘changes to the second centroid
While j k
If dist (Sol ,Ug ) < dmin
j i
‘calculate the distance from object i to Sol (another centroid)
j
dmin  dist (Sol ,Ug )
j i
End If
j  j + 1
‘changes to the next centroid
End While
cost  cost + dmin
End If
i i + 1
End While
Cost(Sol)  cost
98 Journal of Applied Operational Research Vol. 4, No. 2

Results of simulated annealing for territorial optimization
In this section some instances of SA in the clustering algorithm are presented. In this process the initial temperature
To and final temperature Tf are used as control variables, α is the control factor for the likelihood of accepting
solutions worse than a local optimum, L(t) is the control factor for the number of iterations the improvement of a
local solution and Nc the number of groups to consider.
Random experiment
In previous works a variant of SA, called systematic SA, was reported (Bernábe et al, 2008). Notice that when a
new random solution is generated, it is likely to create a solution that has already been tried, the algorithm was
modified so that all neighboring solutions are unique: a systematization of the neighboring generated solutions.
We do not identified significant differences in the quality of the solutions between RS and systematic RS, but the
experience has been important in order to identify initial values for a given statistical experiment used for evaluating
the quality of solutions obtained by the algorithms of clustering that we have implemented (Bernábe et al, 2008).
In this point, two sets of experiments where performed, in the first one the values given to the parameters were
selected randomly. In the instances the parameters were considered according to the specifications of SA: Initial
temperature To: 5000. To must have a high value for the exploration of the solution space to be acceptable, alpha
(α): 0.98. In each loop, the current temperature is multiplied by alpha. With this cooling mechanism, at the beginning
(when the temperature is very high) the system cools down quickly (when there is a greater probability of accepting
a solution that is worse than the current one), when the temperature is low, the value of the current temperature
decreases less each time. If the value of alpha is small, the system cools down quickly, as it gets nearer to zero,
the number of explorations will decrease, even with a high initial temperature. For example, with alpha α = 0.1
only about four iterations will be performed. As the value of alpha gets nearer to one, the number of iterations
will be greater and so more solution space will be explored. Final temperature Tf: was chosen at 0.001. In theory
the algorithm must finish when the current temperature is equal to zero, but since the cooling of the system
is performed by multiplying the current temperature by alpha, its value will be progressively smaller without
reaching zero, so a small value is adequate as finishing criterion.
Statistical experiment
The second set of instances was performed as a statistical experiment that allowed to modeling the effects of the
parameters upon the cost function. Experimental design is a method helpful for obtaining better solutions as a result of
a combination of diverse factors. This method uses the factors themselves, the levels of factors, treatments,
experimental unit, response variable, and test size. Once established the components of Experimental Design,
they must be brought together for their use on diminishing or controlling the experimental error. In the replication
it is possible to estimate precision of the response and randomness allows a homogenization of the result or the
group error obtained. After randomization of the elements, two methods for determining factors dependency must
be applied. In the first method of regression analysis, an independent value is taken, whereas, in the second method
(variance analysis - ANOVA), the independent variables are compared in order to find the difference among them.
The experimental design models are adjusted to the different combinations of factors because the first model
handles the fixed effects in which a fixed factor must be compared with other variables; the second model compares
different variable factors; the third model combines the variable factor with fixed ones adjusting the mixed model
to these factors. In experimental design, we also specify an initial model and acceptation for the models which
certainly helps to have congruency on the experiments, as well as helping for starting the ANOVA method, which
is completed by obtaining the identity of sum of squares, the expected value of sum of squares of the treatment,
the mean square of the treatment, and the mean square of error. Finally, the acceptation or rejection of the result
is determined by a certain criterion.
MBB Loranca and DP Avendaño 99

Box-Behnken
The first experiment developed has used the methodology Box-Behnken (BB). Due to its features, this type of
design is easy to carry out by defining adequate levels for the design parameters, additionally, it is a rotary design, i. e.
with the same variance for all the experimentation points located at the same distance from the design center, and
on the other hand it is possible to perform either sequential experiments to study the individual effects of each
control parameter or simultaneously study their combined effects (Montgomery, 1991).
Another advantage of using this design is that it is possible to model the results with a second order function
and therefore to analyze the behavior of the cost function using the response surfaces methodology.
According to the parameters of SA, five control parameters were chosen, this gave rise to an experiment with
46 experimental runs where four central points where used. The results obtained previously by the heuristic
method were taken into account to choose the levels for the parameters used; this has allowed an experimentation
region to be defined. The levels used in the experiments are presented in Table 1. Related results have also been
reported in (Bernabe et al, 2009). With the aim of representing the behavior of the different experiments for validating
the solutions of geographical/territorial partitioning, we have analyzed the different levels of experimentation, the
random proposal of runs executed with BB, el regression model that adjust the parameters and some conclusions.
The first factorial experiment considered the following values:
Table 1. Levels for the BB experiment.
Param. High Central Low Level
To 5500 5250 5000
Tf 0.1 0.055 0.01
Alpha 0.99 0.985 0.98
L(t) 5 4 3
G 24 18 12

With these levels and the BB design 46 experimental runs were carried out. The notation used in the table is:
test, To (Initial Temperature), Tf (Final Temperature), α (alpha), L(t), G (groups), CF (Cost Function), OP (Optimum).
The number of groups was defined according to previous works where it was determined that 24 is the best
number of groups for the statistical clustering method. Since the experiment requires another two levels, it was
decided to randomly choose two even numbers smaller than 24 (18 and 12). As might be expected, this choice
must be more accurate, and is the subject of ongoing work. The next table 2 shows the details of the experiment
just described.
Table 2. BB experiment.
test To Tf α L(t) G CF test To Tf α L(t) G CF OP
1 5000 0.01 0.985 4 18 13.5279 24 5250 0.1 0.99 4 18 13.7725 11.098
2 5500 0.01 0.985 4 18 13.5878 25 5000 0.055 0.985 3 18 13.6595

3 5000 0.1 0.985 4 18 14.0342 26 5500 0.055 0.985 3 18 13.5348

4 5500 0.1 0.985 4 18 14.1222 27 5000 0.055 0.985 5 18 14.0258

5 5250 0.055 0.98 3 18 13.9166 28 5500 0.055 0.985 5 18 13.0667

6 5250 0.055 0.99 3 18 14.1286 29 5250 0.055 0.98 4 12 16.8496 14.12
7 5250 0.055 0.98 5 18 13.2346 30 5250 0.055 0.99 4 12 17.1076

8 5250 0.055 0.99 5 18 13.8935 31 5250 0.055 0.98 4 24 12.2146 9.279
9 5250 0.01 0.985 4 12 16.2161 32 5250 0.055 0.99 4 24 11.7276

10 5250 0.1 0.985 4 12 16.553 33 5000 0.055 0.985 4 12 16.6959

11 5250 0.01 0.985 4 24 11.5392 34 5500 0.055 0.985 4 12 16.7826

12 5250 0.1 0.985 4 24 12.0292 35 5000 0.055 0.985 4 24 11.8841

13 5000 0.055 0.98 4 18 16.3016 36 5500 0.055 0.985 4 24 11.2403

14 5500 0.055 0.98 4 18 14.1095 37 5250 0.01 0.985 3 18 13.5575

15 5000 0.055 0.99 4 18 13.9159 38 5250 0.1 0.985 3 18 13.2107
100 Journal of Applied Operational Research Vol. 4, No. 2

test To Tf α L(t) G CF test To Tf α L(t) G CF OP
16 5500 0.055 0.99 4 18 13.9545 39 5250 0.01 0.985 5 18 13.6996

17 5250 0.055 0.985 3 12 15.6353 40 5250 0.1 0.985 5 18 14.7604

18 5250 0.055 0.985 5 12 16.0845 41 5250 0.055 0.985 4 18 13.9274

19 5250 0.055 0.985 3 24 12.3314 42 5250 0.055 0.985 4 18 13.8217

20 5250 0.055 0.985 5 24 11.6377 43 5250 0.055 0.985 4 18 13.5833

21 5250 0.01 0.98 4 18 13.5198 44 5250 0.055 0.985 4 18 13.9886

22 5250 0.1 0.98 4 18 14.304 45 5250 0.055 0.985 4 18 13.6392

23 5250 0.01 0.99 4 18 13.3445 46 5250 0.055 0.985 4 18 12.9008

The terms of hardware for testing of PAM were performed on a computer Pentium ® 4 CPU 2400 GHz, 480
MB RAM. In trial 36 there are 24 groups with the parameters of To=5500, Tf=0.055, α=0.985, L(t)=4, a cost of
11.2403 was generated for the objective function, the closest to the optimum obtained by PAM that was 9.279. In
contrast with the time required by PAM to generate the exact solution, which was 7 hours, our algorithm with
3,049 iterations, 2,183 accepted solutions, reduced the computational cost to only 2 seconds.
In the same manner, for trial 17 with 12 groups and the parameters To=5250, Tf = 0.055, α=0.985, L(t)=3 a
cost of 15.6353 for the objective function was generated. The optimum obtained by PAM was 14.120 in 6 hours
time. In contrast with Simulated Annealing with 2278 iterations and 1568 accepted solutions that required a total
of 1 second. The trial 49 for 18 groups, initial temperature 5250, final temperature .055, alpha .985, L(t): 4,
number of iterations 3037, number of accepted solutions 2112. The cost of the objective function with SA was
12.09. With PAM was 11.098. The model that fits the experiment follows.

Estimated Regression Coefficients for FC

Term Coef SE Coef T P
Constant 16963.0 6707.0 2.529 0.018
To -0.5 0.2 -2.360 0.026
TF 358.1 1106.2 0.324 0.749
alpha -31921.4 13443.1 -2.375 0.026
l(t) -16.2 49.8 -0.326 0.747
grupos 6.2 8.3 0.742 0.465
To*To 0.0 0.0 1.301 0.205
TF*TF -0.4 84.0 -0.005 0.996
alfa*alpha 15020.7 6801.6 2.208 0.037
l(t)*l(t) -0.1 0.2 -0.692 0.495
grupos*grupos 0.0 0.0 2.323 0.029
To*TF 0.0 0.0 0.028 0.978
To*alpha 0.4 0.2 2.220 0.036
To*l(t) -0.0 0.0 -0.831 0.414
To*grupos -0.0 0.0 -0.727 0.474
TF*alpha -395.8 1116.3 -0.355 0.726
TF*l(t) 7.8 5.6 1.401 0.173
TF*grupos 0.1 0.9 0.152 0.880
alpha*l(t) 22.3 50.2 0.445 0.660
alpha*grupos -6.2 8.4 -0.742 0.465
l(t)*grupos -0.0 0.0 -1.138 0.266

S = 0.5023 R-Sq = 93.8% R-Sq(adj) = 88.8%

Based on the previous model, we have obtained and reported different results as solution surfaces. We have selected
the instance for 24 groups (clusters) because we found it to have results close to the optimum by varying all the
parameter, i.e., avoiding fixing a particular variable. The real cost of the objective function is 9.27 for 24 clusters,
the minimum value is reached is y = 10.2597 by using the parameters To = 5477.6723, Tf = 0.102, α = 0.980 and
L(t) = 4.9775.