CHEMICAL REACTION OPTIMIZATION METAHEURISTIC FOR LOCATING SERVICE STATIONS THROUGH THE CAPACITATED P-MEDIAN PROBLEM

Chemical Reaction Optimization (CRO) is a metaheuristic for solving optimization problems, which mimics the interactions between molecules in a chemical reaction with the purpose of achieving a stable, low-energy state. In the present work, we utilize the CRO metaheuristic to solve, in an efficient manner, the capacitated p-median problem, in order to locate service stations. Results from solving small to medium-sized problems available in the literature, with up to 724 notes and 200 medians, are compared to their optimal or best-known values. Results show that CRO results are comparable, in terms of accuracy and execution time, to many existing successfulmetaheuristics, as well as exact and hybrid methods, having exceeding those in some cases.


INTRODUCTION
Facility layout and planning is an important topic that has a wide variety of applications in real life.Both private and public sectors are frequently faced with problems involving facility layout decisions.Facility location is concerned in finding the best locations for facilities based on supply-demand requirements.This problem has many applications in real life including locating retail stores, ambulance centers, schools, gas stations, electric vehicles charging stations, hospitals, fire stations, ATM machines, and wireless base stations.Design parameters of the facility location problem may include how many facilities should be sited, where should each facility be located, how large each facility should be, and how should demand be allocated.
Modeling of the facility location problem has been widely investigated in the literature.One of the best-known facility location models is the capacitated p-median problem (CPMP), which is a location problem where a set of objects (e.g., customers) has to be partitioned into a fixed number of disjoint clusters.Each object has an associated weight (or demand) and must be assigned to exactly one cluster.For a given cluster, the median is that object of the cluster from which the sum of the dissimilarities to all other objects in the cluster is minimized.The dissimilarity of a cluster is the sum of the dissimilarities between each object who belongs to the cluster and the median associated with the cluster.The dissimilarity is measured as a cost (e.g., distance) between any two customers.Each cluster has a given capacity, which must not be exceeded by the total weight of the customers in the cluster.The objective is to find a set of medians, which minimizes the total dissimilarity within each cluster.
The CPMP can be found in the literature under various different names, such as the capacitated warehouse location problem, the sum-of-stars clustering problem, the capacitated clustering problem, and others.Another way to look at the problem is to consider it a variation of the classic p-median problem, which was first formulated by Hakimi (1964) and consists of locating p service stations to serve n customers, or nodes, in such a way that the average weighted distance between these customers and their closest stations is minimized.This model, which is also known as uncapacitated p-median problem, has been widely used to locate service stations and was proven, by Kariv & Hakimi (1979), to be NP-Hard.The CPMP extends the original p-median problem by adding a fixed demand to each customer.In addition, a capacity restriction is added to each service station, so that the total demand from all customers, served by a given station, must not exceed its capacity.The CPMP is also known to be NP-Hard (Gary & Johnson 1979).Mulvey & Beck (1984) were the first to extend the uncapacitated p-median problem, by adding a capacity constraint to each service station.In their seminal work, the authors propose a primal heuristic as well as a hybrid method, based on heuristic optimization and subgradients, to achieve good solutions for the problem.Since then, several other researchers have proposed approaches employing exact, heuristic or hybrid methods to provide good quality solutions to the problem, within acceptable computational times.Osman & Christofides (1994) developed a hybrid solution, based on Simulated Annealing and Tabu Search metaheuristics, to provide near optimal to optimal solutions to a group of instances from 50 to 100 nodes and 5 to 10 medians.Maniezzo et al. (1998) presented an evolutionary method combined with an effective local search technique to solve a variety of CPMP problems, including the ones proposed by Osman & Christofides (1994).Baldacci et al. (2002) proposed an exact algorithm for solving the CPMP based on a set partitioning formulation.Lorena & Senne (2003) proposed a local search heuristic for the capacitated p-median problem to be used in solutions made feasible by a Lagrangean/surrogate optimization process, which explores improvements on upper bounds limits of primal-dual heuristics, based on location-allocation procedures that swap medians and vertices inside clusters, reallocate vertices, and iterate until no improvements occur.The authors used instances from the literature as well as real data provided by a geographical information system.A version of a Genetic Algorithm was developed by Correa et al. (2004).
Lorena & Senne (2004) developed a column generation approach integrated with a Lagrangean/ surrogate relaxation to calculate lower bounds.The approach could identify new productive columns, reducing computation time.Computational results were presented on instances created from a geographic database from the city of São José dos Campos, Brazil.Ahmadi & Osman (2005) proposed a metaheuristic called Greedy Random Adaptive Memory Programming Search (GRAMPS) for the capacitated clustering problem.A branch-and-price algorithm for the CPMP was proposed by Ceselli & Righini (2005).Scheuerer & Wendolsky (2006) proposed a scatter search heuristic for the capacitated clustering problem.It was evaluated on instances from the literature, obtaining several new best solutions.Diaz & Fernandez (2006) proposed a hybrid scatter search and path relinking heuristic for the same problem.The authors ran a series of computational experiments evaluating the proposed methods on instances from the literature, including instances corresponding to 737 cities in Spain.Chaves et al. (2007) presented a hybrid heuristic called Clustering search (CS), which consisted in detecting promising search areas based on clustering.Boccia et al. (2008) proposed a cut-and-branch approach, which proved to be effective in solving hard instances, using IBM CPLEX, or reducing their integrality gap.Fleszar & Hindi (2008) developed a hybrid heuristic that utilizes Variable Neighborhood Search to find suitable medians, thus reducing the CPMP to a generalized assignment problem, which was then solved using IBM CPLEX.Stefanello et al. (2015) presented a matheuristic approach, which consisted of reducing mathematical models by heuristic elimination of variables that are unlikely to belong to a good or optimal solution.Additionally, a partial optimization algorithm based on their reduction technique was proposed.Resulting models were solved by IBM CPLEX, with good accuracy and performance.
Chemical Reaction Optimization (CRO) is a recently created metaheuristic for optimization, inspired by the nature of chemical reactions, proposed by Lam & Li (2012).A chemical reaction is a natural process of transforming unstable substances into stable ones.Under a microscopic view, a chemical reaction starts with some unstable molecules with excessive energy.These molecules interact with each other through a sequence of elementary reactions.At the end, they are converted to molecules with minimum energy to support their existence.This property is embedded in CRO to solve optimization problems.In this paper, we adapt the CRO metaheuristic to solve the CPMP problem.Additionally a simple heuristic, based on the first size-reduction heuristic proposed by Stefanello et al. (2015) as part of their "Iterated Reduction Matheuristic Algorithm" (IRMA), along with a modified version of the λ-interchange mechanism, presented by Osman & Christofides (1994) is used during intensification (local search) phases.This adapted version of CRO is used to solve a wide variety of instances available in the literature, with sizes ranging from 50 to 724 customers and 5 to 200 medians.Results show that CRO can be effectively used with the capacitated p-median problem, achieving good results in terms of accuracy and execution time.This paper is organized as follows: Section 2 provides a formal description of the CPMP for locating service stations.Section 3 describes the CRO metaheuristic, focusing on the customizations that we implement to solve the CPMP, including a simple constructive heuristics and a λ-interchange mechanism used during the local search phase.Section 4 contains the computational results of our studies for all tested instances.Conclusions and future developments are presented in Section 5.

THE CPMP PROBLEM
The CPMP problem, depicted in this paper, aims to locate service stations, from a set of candidate locations and a set of customers.The CPMP integer linear programming model is shown below: Decision variables: x j = 1 if candidate station j was selected; or 0, otherwise; (1) Model: min subject to: x j ∈ {0, 1}, ∀ j ∈ J (8) where: J = set of candidate service station locations (median candidates); I = set of customers (nodes); d i j = distance, or cost, from customer i to service station j ; p = number of medians or stations to be opened; Di = demand associated to customer i; cj = capacity of service station j .
Function (3) represents the optimization function.The objective is to minimize the total cost, or distance, from all customers to their assigned service stations.Constraint (4) requires the number of selected stations to be equal to p, whereas, constraint (5) requires the demands from all customers to be met.Constraint (6) ensures that a customer is associated only to a service station that is selected.Constraint (7) makes sure that the total demand from all customers assigned to a station does not exceed its capacity.Finally, constraints ( 8) and ( 9) define the domain of the decision variables x and y.Without loss of generality, we assume that J = I , in this paper.

CRO METHAHEURISTICS FOR THE CAPACITATED P-MEDIAN PROBLEM
The CRO metaheuristic is a technique developed by Lam & Li (2012), which loosely relates chemical reactions with optimization and is based on the first two laws of thermodynamics.The first law, the conservation of energy, states that energy cannot be created or destroyed, but only transformed or transferred from one entity to another.The second law states that entropy, which is the measure of the degree of disorder of a system, tends to increase.
A chemical reaction system consists of the chemicals substances and their environment.The energy of the environment is symbolically represented by a central energy reservoir, i.e., a buffer.A chemical substance is comprised of molecules, which possess potential and kinetic energy.
A chemical reaction occurs when the system is unstable, due to excessive energy.All chemical reaction systems tend to reach a balanced state, in which potential energy drops to a minimum.CRO simulates this phenomenon, by gradually converting potential energy into kinetic energy and transferring energy from the molecules to the environment through consecutive steps or sub-reactions, over several transition states, which result in compounds that are more stable and contain minimal energy.It is an iterative process that seeks the ideal point.
A collision provokes a chemical change in a molecule.There are two types of collisions in CRO: unimolecular and intermolecular.The first ones describes the situation in which a molecule collides with the wall of a container, whereas the latter represents the cases in which a molecule collides with other molecules.Such chemical change is called an elementary reaction.An ineffective elemental reaction is one that results in a subtle change in molecular structure.CRO utilizes four types of elementary reactions: on-wall ineffective collision, decomposition, intermolecular ineffective collision and synthesis.Decompositions and syntheses cause much more vigorous changes in molecular structures.Elemental molecular reactions are summarized in Table 1.In the present work, CRO is implemeted using C# object-oriented programming language, due to the easiness of modeling molecules as instances of a class which contains all the attributes needed for its operation.The molecule and elementary reactions, as well as the main algorithm for the CRO are implemented as methods of classes.The following subsections describe the main components of the CRO focusing on the modifications that are done to solve the CPMP.Further information about the CRO metaheuristics can be obtained in Lam & Li (2010, 2012).

The Molecule
The basic unit of the CRO algorithm is the molecule, which contains several attributes that are essential to its proper operation.For this implementation of CRO, which we call CRO for the CPMP, we define the following attributes: • Molecule ID (MolID): uniquely identifies a molecule in the population of molecules.
• Molecular structure (ω): stores a feasible solution for the problem, which is comprised of the objective function value, as well as the decision variables x(1) and y (2).The x decision variable set stores the current station selection for a given feasible solution.It is implemented as a list of integers, thus storing only the station numbers that are part of the solution.For instance, if the number of candidate stations is 100 and the number of medians is five, x will store five values corresponding to the currently selected stations, such as the list {5, 17, 29, 45, 79}.This has proven to be more effective than storing x as an array of bits of size |J |.Similarly, the decision variable set y, which stores the customer-to-station assignments, is implemented as an array of integers of size |I |, instead of a matrix of bits of size |J | × |I |.
• Potential Energy (PE): it is defined as the objective function value in the molecular structure (ω).If f denotes the objective function, then P Eω = f (ω).
• Kinetic Energy (KE): it is a non-negative number that quantifies the tolerance of a system to accept a solution that is worse than an already existing one.
• Number of Collisions (NumHit): total number of hits (collisions) that a molecule has taken.• Minimum Hit Number (MinHit): it is the collision number when MinStruct was achieved.
• Proximity List Set (PL): it is a set of size |J |.Each element P L j of the set contains a list of stations that are near station j .These lists are populated using a strategy originally presented by Stefanello et al. (2015), for their "Iterated Reduction Matheuristic Algorithm" (IRMA) R1 heuristic.The strategy is modified to provide a list of nearby candidate stations that can replace a selected median, during the intensification phases of CRO (unimolecular and intermolecular collisions), with the purpose of reducing the total number of iterations required by the λ-interchange mechanism.The process of building the Proximity List Set is capacity-based and shown in details in section 3.4.
• κ: it is a parameter used as an expand capacity factor to control the number of nearby candidate stations stored in each of the lists of the Proximity List Set (PL).If f denotes the objective a function that generates a list of nearby stations for each candidate station, then The pseudocode for the "Molecule" class is shown in Figure 1.It contains only properties and a constructor.The molecule's constructor is called from a constructive algorithm, responsible for generating a number of random feasible solutions, and receives a molecule ID number (MolID), a feasible solution and an initial amount of kinetic energy (initialKE).In the constructor's code, a new Proximity List Set with initial capacity factor κ0 is created.

Initialization and constructive phase
The initialization phase consists in setting appropriate values for CRO's operational parameters, as defined in Lam & Li (2012).These parameters are PopSize, KELossRate, MoleColl, buffer, InitialKE, α and β.A brief description of them is given below: • PopSize: it is an integer number that denotes the number of molecules of the initial population.Due to the size of most of the instances evaluated in this paper, small initial populations, from 2 to 10 molecules are used, in order to reduce initialization times.
• KELossRate: used by CRO during on-wall ineffective collisions to determine the minimum amount of kinetic energy that a molecule will retain from its initial energy, after a collision.We set KELossRate to 0.8 for all test instances.
• MoleColl: a real number that denotes the probability of an intermolecular collision to occur.A MoleColl value of 0.1 is used on all test instances.
• buffer: CRO's central energy buffer.Its initial value is set to zero on all test instances.
• InitialKE: denotes the amount of kinetic energy that is given to a new molecule.
• α: sets a limit for the number of times a molecule can undergo a local search without locating a better local minimum, before it becomes entitled for a decomposition.
• β: molecules with too low KE lose their flexibility of escaping from local minima.β denotes the minimum energy that a molecule needs to qualify for a synthesis reaction.
The constructive algorithm creates a new population of molecules containing feasible solutions, i.e., solutions that do not violate constraints (4) to ( 9) of the CPMP model.It is a five-step process based on a primal heuristic proposed by Mulvey & Beck (1984), a neighborhood search algorithm proposed by Maranzana (1964) and a Fast Interchange algorithm, proposed by Whitaker (1983).The process is depicted below: a) A preliminary population 100 times larger than the desired population (PopSize) is created.p stations, or medians, are randomly chosen from the candidate set.Then, a primal heuristic proposed by Mulvey & Beck (1984) is applied to all solutions, as follows: customers are assigned to the selected medians in a decreasing order of their regret value.Regret is defined as the absolute value of the difference, in terms of distance, between the first and second nearest median.
b) Once step a is complete, the molecules population is reduced to PopSize molecules containing the best objective function values among the preliminary population.The remaining molecules are discarded.This procedure, which does not require a higher computational effort, has proven to be effective in increasing quality of solutions.c) Based on the assumption that a set of medians selected for a good capacitated p-median solution should be close to a set of medians selected for a good uncapacitated p-median solution, with respect to the distance between these sets, we optimize the solutions created on step b, ignoring the capacity constraint.In this step, we search for the best local uncapacitated minimum, using a Fast Interchange algorithm, proposed by Whitaker (1983) and further modified by Hansen & Mladenovic (1997).This algorithm is briefly explained in section 3.3.
d) Once a local minimum is found, we apply the same heuristic described on item a, to the uncapacitated solution, in order to obtain a capacitated one.
e) Next, the constructive algorithm examines each one of the nodes assigned to a median, in an attempt to find a candidate node that minimizes the demand-weighted total cost (or distance) d, among all nodes within a particular set.If a better node than the already selected node (median) is found, that node becomes the new median and all other nodes of that set are reassigned to it.At this point, the capacity constraint ( 7) may be violated if the newly found median does not have enough capacity to serve all nodes in the cluster.This procedure was proposed by Maranzana (1964) as part of their neighborhood search algorithm.
f) Once step e is complete, nodes are again assigned to their selected medians in a decreasing order of their regret value, as described on item a. Steps d to e are repeated until a limit of 20 iterations without improvement is reached.
Additionally to obtaining feasible solutions for all molecules, a new set of Proximity Lists, one list for each molecule, is created by the constructive algorithm.The process is described in section 3.4.

The Fast Interchange algorithm
We utilize the Fast Interchange algorithm, proposed by Whitaker (1983) and further modified by Hansen & Mladenovic (1997), as the intensification algorithm for all new solutions created by the constructive algorithm and solutions that have been partly modified by decomposition and synthesis processes.In Whitaker's work, three efficient ingredients are incorporated into the standard interchange heuristics: • Move evaluation, where the best station to be removed is found when the station to be added is known.
• Update of the first and second closest stations to each customer.
• Strategy restricted to first improvement, where each station is considered to be added only once.
We use a version modified by Hansen & Mladenovic (1997) and we obtain the best improvement, instead of the first improvement.
Initially, the algorithm assigns the nearest customers to each station, thus relaxing the capacity constraint ( 7) from the CPMP model.Then, the Fast Interchange operator is used to improve the uncapacitated solution, in order to obtain the best improvement.Finally, the primal heuristic from Mulvey & Beck (1984) is applied to the incumbent solution, as shown in item a of section 3.2.

The Proximity List Set
The Proximity List Set (PL) is a data structure used to limit the number of iterations performed by a λ-interchange mechanism, during the intensification phase of CRO.The λ-interchange mechanism is described in section 3.5.3.Due to the size of some of the instances that are evaluated in our tests, it becomes necessary to employ a strategy to reduce the number of interchange and shift operations that are typical of such algorithm.That is especially important on instances where p is high.In a traditional λ-interchange, every customer or group of customers is systematically relinked to all selected stations, other than the one they are currently assigned to, one station at a time, looking for possible improvements in the objective function value.In our version, we only consider for relinking, the selected stations that are in the vicinity of the station that a customer, or a group of customers, is currently assigned to.To accomplish that, for every node we build a static list containing nodes that are near that particular node.We call it a Proximity List.
A Proximity List Set is a set of Proximity Lists of size J .Each element P L j contains a list of stations that are near station j .The method used to populate the lists of nearby stations is an adaptation of an heuristic proposed by Stefanello et al. (2015), which was originally designed to eliminate decision variables that are unlikely to belong to good or optimal solutions.Consider a decision variable x j .We define a subset P L j ⊆ J of the nearest nodes of j as: where a node t is nearer to j than t if d jt < d jt .Thus, for a candidate median j ∈ J , we include the variable x t in the nearby station list if t ∈ P L j .The parameter κ is an expand capacity factor used to control the size of the proximity lists.Figure 2 shows a feasible solution for one of the instances (lin318 040) we test in this paper.The nodes that are part of the Proximity List of station 122 for κ = 5 are marked with a cross.The dashed line demarcates the boundaries of it, i.e., the farthest nodes from station 122.For the feasible solution shown in the figure, three stations would be considered, by the λ-interchange mechanism, for relinking the customers of station 122.A higher value for κ would include other nearby stations and vice-versa.
A partial pseudocode for a class (ProximityListSet) which stores a proximity list set is shown in Figure 3.The class constructor takes only one parameter, κ, which indirectly controls the size of each of the proximity lists in the set.For every candidate station j , we build an list of edges containing the distance from j to every node i, their distance to candidate station j (d i j ), as well as their demand (D i ).Then, the resulting list of edges is sorted by d i j , in ascending order.Finally, a proximity list is populated picking the first t elements of the list that satisfy (10).

Elementary Reactions
There are four types of elementary reactions that can occur at every CRO iteration.They are used to manipulate solutions (explore the solution space) and redistribute energy between the molecules and the energy buffer.Operators are used to modify solutions or generate new solutions from current solutions.During all operations, energy conservation is always maintained.
For on-wall and intermolecular ineffective collisions we implement a neighborhood search operator based on the λ-interchange mechanism, proposed by Osman & Christofides (1994).For synthesis reactions, a crossover operator commonly used in Genetic Algorithms is used.Finally, in decompositions, we use a Half-Total Change operator.These operators are explained in details in the following sections.If a solution produced by any of these operators causes a capacity constraint violation (7), the solution is rejected.Additionally, for a solution to be accepted, it must pass CRO's own acceptance energy-based mechanisms, which are succinctly described in this paper.More information about CRO's inner workings can be found in Lam & Li (2010, 2012).

On-wall ineffective collision
It happens when a molecule collides with the wall of a container and bounces back, remaining a single molecule.In this type of collision, the existing solution ω is perturbed to become ω , i.e., ω → ω .This is done by generating a ω that is in a neighborhood of ω, through a λ-interchange operator, which was proposed by Osman & Christofides (1994).
Let N (•) be a λ-interchange neighborhood search operator.Therefore, we have ω = N (ω) and P Eω = f (ω ).In this type of reaction, typically, a potential energy loss will occur, i.e., P Eω will be less than P Eω, indicating that a better solution has been obtained.If that does not occur and P Eω is greater than P Eω, then the worse solution can still be accepted, provided that P Eω + K Eω ≥ P Eω .However, every time a reaction occurs, a certain amount of kinetic energy (K E) is transferred to the energy buffer, decreasing the likelihood that worse solutions will be accepted on further iterations.The amount of kinetic energy of the molecule obtained from the ineffective reaction is controlled indirectly by the parameter KElossRate, which is a value between 0 and 1, inclusive, and affects the minimum amount of kinetic energy that is withdrawn from the original solution (ω).

Intermolecular ineffective collision
It occurs when two molecules collide with each other and then bounce away.The number of molecules remains unchanged, i.e., ω 1 + ω 2 → ω 1 + ω 2 .This reaction is very similar to the on-wall ineffective collision, thus the same λ-interchange operator from the on-wall ineffective collision is utilized.Let N (.) be a λ-interchange operator.Therefore, ω 1 and ω 2 are obtained through ω 1 = N (ω 1 ) and ω 2 = N (ω 2 ).Energy management is similar to on-wall ineffective collisions, but does not involve the energy buffer.

The λ-interchange neighborhood search mechanism
For all on-wall and intermolecular ineffective collisions we implement a neighborhood search operator based on the λ-interchange mechanism, proposed by Osman & Christofides (1994) with possibly new different medians ζ i and ζ j , respectively.The new solution becomes S = {C 1 , . . ., C i , . . ., C j , . . ., C p }.The neighborhood N (S) of S is the set of all S solutions generated by the λ-interchange mechanism for a given integer λ and it is denoted by N λ (S).
Let P L i be a proximity list containing nodes near ζ i (the median of cluster C i ) and m i be the number of nodes in P L i which are also medians.The λ-interchange mechanism will only examine the pairs of clusters (C i , C j ) where ζ j ∈ P L i .Therefore, the number of different pairs of clusters (C i , C j ) to be examined is m i , for a given λ.
For any given pair of clusters (C i , C j ), the λ-interchange mechanism utilizes two processes to generate neighborhoods.Let μ be the number of customers from C i or C j to be handled by any of these two processes: a shift process tries to move μ customers from cluster C i to C j , or vice-versa.For μ = 1, a shift process is represented by the (0, 1) and (1, 0) operators.An interchange process, as the name implies, attempts to swap every customer from the first cluster with every other customer in the second cluster and, for μ = 1, it is represented by the operator (1, 1).We implement a different order of search than the one employed by Osman & Christofides (1994).First, we attempt to swap customers before trying to shift them.This has proven to be more effective in our tests.Therefore, the order of operators becomes (1,1), (1,0) and (0,1).
Our implementation of λ-interchange mechanism is described as follows: Upon start, a Proximity List Set is recomputed for an initial κ(κ 0 ).Then, starting with a feasible solution S, the λ-interchange logic is executed for a specific number of iterations, in an attempt to improve the current solution.In the end, the best solution is returned.The iterative process of the λ-interchange is shown below.
Let μ, 1 ≤ μ ≤ λ, be the number of customers from cluster C i to be moved by a shift or interchange process and μ , 1 ≤ μ ≤ λ, be the number of customers to be shifted or interchanged from cluster C j .During a single iteration, and starting with μ = 1 and μ = 1, a search randomly examines all possible cluster pairs (C i , C j ), without repetition, looking for groups of μ customers from C i and μ customers from C j that could be interchanged or shifted.This process continues until all cluster pairs are evaluated for all possible values of μ and μ , without repetition.If no improvement in the objective function (3) is achieved, then it is likely that a local minimum has been reached.In order to escape from that possible local minimum, the Proximity List Set (PL) is recomputed for a new κ, by adding a κ to it.As the proximity lists' sizes grow, new neighborhoods may be reached.To prevent the size of each Proximity List Set to grow too big, thus defeating their main purpose of reducing the overall number of shift and interchange operations, we only let κ grow until the average size of all proximity lists reach 20% of the number of candidate stations, J .When that limit is reached and if there is still no improvement in the objective function, one median ζ i from the incumbent solution S is randomly interchanged with another node ζ j |ζ j / ∈ S, ζ j ∈ P L i .Both measures combined should allow new neighborhoods to be reached on subsequent iterations.
The process finishes when a specific number of iterations without improvement is reached.The number of iterations is proportional to the molecule's NumHit, which is the total number of hits (collisions) that a molecule has suffered, so far.The actual number of iterations, for a given execution of the λ-interchange mechanism, is determined by multiplying NumHit by a minimum number of iterations (Minλ-InterchangeIterations), which is typically a small integer, such as 1 or 2. This allows the first ineffective reactions to occur very fast, ruling out molecules that do not carry promising solutions.These molecules will quickly become targets for decompositions and syntheses, as they will not show improvements fast enough and will lose kinetic energy rapidly.On the other hand, molecules that carry better solutions tend to survive longer, as with every new collision they will be given a higher iteration limit.
According to Osman & Christofides (1994), a solution S is λ-optimal (λ-opt) if, and only if, for any pair of clusters C i , C j ∈ S, there is no improvement that can be made by any λ-interchange move.The authors also state that, in order to produce an efficient λ-opt descent algorithm it is advisable to produce a 1-opt solution.In our version of the λ-interchange mechanism, we follow this recommendation, by always obtaining a 1-opt solution before improving it to a 2opt solution.Due to high computational costs of such mechanism, we limited λ to a maximum value of two.All on-wall collisions obtain a 1-opt solution, whereas intermolecular ineffective collision obtain a 2-opt solution, if necessary.The approximate rate of 1-opt to 2-opt solutions is 10, and it's controlled by setting CRO's parameter Molecoll to 0.1.The pseudocode for the λ-interchange mechanism is shown in Figure 5.

Decomposition
A decomposition occurs when a molecule (ω) collides with a wall and then breaks down in two pieces, producing ω 1 and ω 2 , that is: ω → ω 1 + ω 2 .The purpose of decomposition is to allow the system to explore other regions of the solution space, after having made considerable local search through ineffective collisions.Since more solutions are created, the total sum of P E and K E of the original molecule may not be sufficient.In other words, we can have P Eω + K Eω < P Eω 1 + P Eω 2 .Since energy conservation is not satisfied under these conditions, this decomposition must be aborted.To increase the chance of having a complete decomposition, a small portion of the energy buffer is withdrawn to support the change.
We use the Half-total Change (HTC) operator to generate new solutions through decomposition.As the name implies, a new solution is produced from an existing one, keeping half of the existing values (medians) while assigning new values to the remaining half.Suppose we try to produce two new solutions To obtain ω 1 , first we copy the whole solution (ω) to ω 1 .Then, we randomly select {N /2} elements from ω 1 , where { .} returns the largest integer equal to or less than the argument.For each of these elements, a new median is chosen at random, as long as it does not violate any of the problem constraints.Since elements are chosen randomly, ω 1 and ω 2 are very different from each other and also from ω .In the present implementation, the half-total change operator ensures that the elements of ω are present on either ω 1 or ω 2 , but not on both.Therefore, all medians from ω will exist either on ω 1 or ω 2 .
The HTC operator takes as input a molecule's Minimum Structure (MinStruct), which stores the best solution achieved by the molecule, and returns two output solutions, ω 1 and ω 2 .Initially, the process described above is repeated 100 times for each of the output solutions, ω 1 and ω 2 .Then, we choose the best ω 1 and ω 2 .Before the newly created solutions are returned to the main algorithm, both solutions are improved using the same process depicted on items 3.2b to 3.2e.

Synthesis
It is the opposite of decomposition.A synthesis occurs when two molecules collide and are fused together, that is: ω 1 + ω 2 → ω .In this reaction, a much larger change is allowed for ω with respect to ω 1 and ω 2 , along with a considerable increase in the kinetic energy of the resulting molecule.Therefore, it has a greater ability to explore its own solutions space, due to its higher kinetic energy.
We use a Distance Preserving Crossover (DPX) operator to carry out the synthesis of two molecules in a single molecule.The DPX operator was used by Merz & Freisleben (1997) to solve the quadratic assignment problem (QAP) through a Genetic Algorithm.It has proven to be well adapted to the CPMP problem.Let π 1 and π 2 be valid solutions for a given problem.The distance T , between solutions, is defined as: As defined in (11), T represents the number of selected stations present in π 1 , which are not in π 2 .The DPX operator aims to produce a descendant that has the same distance from each of its parents, being that distance equal to the distance between the parents themselves.Let A and B be two parents, with both A and B containing feasible solutions.First, all stations (medians) which are present on both parents are copied to descendant C. The remaining positions are then randomly populated with stations not yet assigned to any of its parents, ensuring that no station found in only one parent is copied to the descendant.This way, we obtain a descendant C, for which the condition T (C, A) = T (C, B) = T (A, B) holds.Such crossover is highly disruptive, forcing subsequent local searches to explore different regions of the solution space, where better solutions could be found.If a feasible solution cannot be obtained, the synthesis fails.
The DPX operator takes as input the Minimum Structure (MinStruct) of two molecules, ω 1 and ω 2 , containing the best solutions achieved by these molecules, and returns a single molecule, ω .Initially, the process described above is repeated 100 times.Then, we choose the best ω , i.e., the one with the lowest P E. Before the newly created solution is returned to the main algorithm, it goes through an improvement process, which is the same as the one depicted on items 3.2b to 3.2e.

Energy conservation
In CRO, energy cannot be created or destroyed.The total energy of the system is determined by the objective function values, i.e., the P E of the initial molecules' population, whose size is determined by PopSize, the initial K E assigned to them and the initial energy value of the buffer.
In all experiments, we set the buffer's initial value to zero.

Initialization
Upon start, we set scalar variables corresponding to CRO standard operational parameters, which are PopSize, KELossRate, MoleColl, buffer, InitialKE, α and β.In addition to CRO's standard parameters, we set a few other parameters specific to our implementation, as follows: • λ: it controls the highest λ-opt solution to be obtained by the λ-interchange mechanism.For example, if λ = 2, only 1-opt and 2-opt solutions will be computed by the mechanism.
• Minλ-InterchangeIterations: it is the starting number of iterations without improvement to be executed by a λ-interchange mechanism.The actual number of iterations for a given molecule is Minλ-InterchangeIterations * (NumHit + 1).
• MinMol: it is the minimum number of molecules in the population.If the population reaches this threshold, no syntheses will occur.
• MaxMol: it is the maximum number of molecules in the population.If the population reaches this threshold, no decompositions will occur.
• MaxIterations: it is the maximum number of iterations for the main CRO algorithm.
• MaxIterationsWithoutImprovement: it is the maximum number of iterations without improvement in the objective function, for the main CRO algorithm.
• κ 0 : it is the initial value for the Proximity List Set capacity factor (κ).It is used in the λ-interchange mechanism to indirectly influence the size of the lists in the Proximity List Set, as explained in sections 3.4 and 3.5.3.
• κ: capacity factor increment, used by the Proximity List Set in the λ-interchange mechanism, as explained in section 3.5.3.
Next, a constructive algorithm, explained in section 3.2, is invoked to create an initial population of PopSize molecules.

Iterations and finalization
During a single iteration, a molecule can collide with the wall of a container or with another molecule.This is decided by generating a random number b between [0, 1].If b > MoleColl, a unimolecular collision occurs.Otherwise, an intermolecular collision takes place.In a unimolecular collision, we randomly select a molecule from the population and decide if an on-wall inefficient collision or a decomposition occurs, according to a decomposition criterion, which is defined as: NumHit − MinHit > α (12).For an intermolecular collision, two molecules are randomly selected and then we determine if an intermolecular collision or a synthesis happens by checking the synthesis criteria for the chosen molecules, which is defined as K E = β (13).
The inequalities ( 12) and ( 13) control the degree of diversification through a and β parameters.Adequate α and β values allows a good balance between diversification and intensification.
After an elementary reaction occurs, we check if the energy conservation condition is obeyed.If this has not happened, the change is discarded.Next, we check if the solution produced by the collision has a lower objective function value than the best solution we have obtained in the population, thus far.If so, we replace the best solution with the incumbent solution.
If none of the stopping criteria is reached, we begin a new iteration.Otherwise, we exit the main loop and return the best solution found.The number of iterations is controlled by MaxIterations and MaxIterationsWithoutImprovement, as described in section 3.6.The main algorithm for the CRO is shown in Figure 6.
The pseudocodes for both ineffective reactions, synthesis and decomposition have not changed significantly from the tutorial, which the present implementation is based upon, and can be found on Lam & Li (2012).

Benchmark Datasets
For the evaluation of the proposed solution, a number of computational experiments are performed using benchmark data sets taken from the literature.We use five different data sets, with number of nodes ranging from 100 to 724 and number of medians ranging from 5 to 300.
The first dataset was proposed by Osman & Christofides (1994) and has been widely used for benchmarking performance and accuracy of CPMP solutions.The first group of 10 instances from this dataset, named cpmp01 to cmpm10, has 50 nodes and 5 medians, whereas the last group, also with 10 instances, named cpmp11 to cmpm20, has 100 nodes and 10 medians.Results from our implementation of CRO for the CPMP are compared in accuracy with the ones obtained by Osman & Christofides (1994), utilizing a hybrid metaheuristic (SATS) containing elements of Simulated Annealing and Tabu Search.Additionally, the accuracy and computational time of our solution is compared with results obtained from solving the same problems using IBM CPLEX vs.12.6 on the same hardware and, lastly, with the matheuristic IRMA, proposed by Stefanello et al. (2015), on similar hardware and MIP solver (Intel i5 and CPLEX 12.3).The second dataset, proposed by Lorena & Senne (2004), is comprised of six instances, named sjc1 to sjc4b, built from data gathered from a geographic database of the city of São José dos Campos, Brazil.The number of nodes vary from 100 to 400 and the number of medians from 10 to 40.We compare the accuracy of our CRO for the CPMP with results available in the recent literature from Scheuerer & Wendolsky (2006) who developed a Scatter Search heuristic (SS) and with a Variable Neighborhood Search combined with CLEX from Fleszar & Hindi (2008).
We also compare with a Clustering search heuristic from Chaves et al. ( 2007) and a Fenchel cutting planes allied with CPLEX approach (Fen-CPLEX) from Boccia et al. (2008).Additionally, the accuracy and computational time of our solution is compared with results obtained from solving the same problems using IBM CPLEX vs. 12.3 and the matheuristic IRMA, proposed by Stefanello et al. (2015), on similar hardware.
The last three datasets we utilize for accuracy and computational time comparisons, was proposed by Stefanello et al. (2015) and is originally comprised of 6 data sets of 5 instances each, adapted from TSP-LIB, with number of nodes varying from 318 to 4461 and number of medians varying from 5 to 1000.Since the present work is focused on instances of up to 724 nodes, we solve only the first three data sets.Larger instances will be tackled by a parallelized version of our CRO for the CPMP in a future article.However, most of the instances from the first three datasets cannot be solved to optimality by MIP solvers in a timely fashion, which justifies the use of heuristics.
All algorithms are coded in C# programming language and run on an Intel i7 2.3 GHz PC with 16 GB of RAM.After proper tune up, we solve all instances 20 times, recording execution times and objective function values.We also record the best objective and standard deviation.
In addition to solving the instances using our CRO for the CPMP, the first dataset is solved by CPLEX vs.12.6, running on the same hardware.

Parameter tuning procedures
Our parameter tuning procedures are developed based on the premise that a good heuristic should provide acceptable results, that is, optimal or near optimal objective function values in rather short computational times.Therefore, during our tests we privilege low execution times over achieving optimality.Moreover, we try to achieve a gap of less than 1% on all tested instances.In the present work, we define gap as the percentage difference between the objective function values obtained by CRO, or any other metaheuristic or MIP solver, and the best-know value for a given instance.
The CRO is a highly parameterized metaheuristic.There are seven parameters specified in Lam & Li (2012) that may affect accuracy and execution times, or both.These parameters are Pop-Size, KELossRate, MoleColl, buffer, InitialKE, α and β.Other parameters that we introduce that may affect performance are intended to provide stopping criteria (MaxIterations and MaxItera-tionsWithoutImprovement) and control intensification during the execution of the λ-Interchange mechanism (λ and Minλ-InterchangeIterations).Molecule population minimum and maximum sizes are controlled by MinMol and MaxMol parameters, respectively.Finally, κ 0 and κ control the size of the Proximity List Set (P L) lists.Therefore, there is a total of 15 parameters that can affect CRO's performance.Since the combination of parameters exist in a fifteen-dimension space, it is impractical to test all possible combinations.Instead, the parameters are tuned in an ad-hoc manner.Following is a brief discussion on the methodology we use for tuning and some empirical findings that may help in tuning similar CRO implementations.
We have empirically determined that the best results are achieved when the number of syntheses and decompositions is below 10% of the total number of collisions.This is in accordance with Lam & Li (2012), who state that "decomposition and synthesis bring diversification to the algorithm.Diversification cannot take place too often, or CRO will be become a completely random algorithm".We have also empirically found that, for a given KELossRate, MoleColl and InitialKE, the number of decompositions is highly dependent on α.If α is too low, e.g. less than two, molecules carrying solutions, may not have enough time (expressed in number of collisions) to achieve local minima.On the other hand, if a is too high, decompositions may never occur.Similarly, the number of syntheses is highly dependent on β.If it is too close to InitialKE, syntheses will happen more often than desired, thus hurting intensification.If β is too low, the number of synthesis collisions may not suffice.
Another important consideration when choosing appropriate values for the various CRO parameters is population control.If the number of syntheses and decompositions is unbalanced, population will rapidly reach MinMol or MaxMol.
We choose an InitialKE that is about the 10 to 20 times higher than the best-known objective for the problems we tested.When InitialKE is too low, most decompositions may be rejected, due to a lower tolerance to accept poorer solutions.Similarly, syntheses may start occurring too soon.
Finding a suitable value for MoleColl is also very important.We choose 0.1 for all of our experiments, which means the probability of intermolecular collisions to happen is 0.1.Since the number of on-wall collisions is much higher, we always run a 1-opt (λ = 1 or 1-interchange) optimization for those.That has proven to be adequate for several instances to achieve optimality.For intermolecular collisions, we run a 1-opt optimization or a 2-opt optimization if we cannot achieve satisfactory results.Running a 2-interchange logic requires considerably more time, as the number of comparisons increase substantially.In addition, a 2-opt optimization is executed for each of the two molecules involved in the collision.Using a higher λ is impracticable due to the dramatic increase in execution times.
To tune the parameters which provide stopping criteria (MaxIterations and MaxIterationsWith-outImprovement) we solve a particular instance 20 times using our CRO for the CPMP, starting from a rather small value for MaxIterations and MaxIterationsWithoutImprovement, such as 50 iterations.If an average gap, for the objective function, of less than 1% is obtained and at least one run achieves the best-know value for the objective function, we use that number of iterations as MaxIterationsWithoutImprovement and set MaxIterations to be twice as many.Otherwise, we keep increasing MaxIterations and MaxIterationsWithoutImprovement by a small amount, such as 50 iterations, until a gap of less than 1% with one optimal (or best-known) run is obtained or there is no significant improvement in the objective function.If those criteria cannot be fulfilled, we chose the run with the lowest gap.Additionally, we first try to solve an instance running a 1-interchange logic for all on-wall and intermolecular collisions.If the desired gap could not be obtained, then we use a 2-interchange logic for intermolecular collisions only.This is to avoid the performance penalty mentioned earlier.Once we find suitable values for MaxIterations and MaxIterationsWithoutImprovement, we solve the instance another 10-30 times, depending on the data set, and use the results in our comparison with other metaheuristics.
Many other tuning methodologies may be used to tune CRO parameters.The one we present is merely one that works for the instances we test.It generates good solutions over a large spectrum of problem parameters in our empirical tests.It may be possible to find a methodology that generates better solutions or works faster.Therefore, we do not claim that our methodology is "optimal" in any sense.
As an example of the tuning methodology we adopt, Figure 7 shows the objective function gap as a function of MaxIterations, for the second benchmark dataset, proposed by Lorena & Senne (2004), which is comprised of six instances named sjc1, sjc2, sjc3, sjc3a, sjc4a and sjc4b.We use the same value for MaxIterations and MaxIterationsWithoutImprovement on all executions.MaxIterations varies from 50 to 500 iterations.The following CRO parameters are fixed: PopSize = 10, KeLossRate = 0.8, MoleColl = 1, InitialKe = 1, 000, 000, A = 10, B = 5, 000, MinMol = 2, MaxMol = 20 and Minλ-InterchangeIterations = 1.We run a 1-opt optimization (λ-interchange = 1) on all instances, except for instance sjc4a, which we run a 2-opt optimization for all intermolecular collisions.For the same dataset, Table 2 shows the instance name, number of nodes, number of medians ( p), MaxIterations, MaxIterationsWithoutImprovement and optimal value for the objective function as well as the lowest objective achieved by our CRO for the CPMP, the average percentage gap and the standard deviation for 20 runs of each combination of instance and MaxIterations (or MaxIterationsWithoutImprovement).Lines that meet the aforementioned criteria for tuning MaxIterations and MaxIterationsWithoutImprovement (average gap of 1% or less with at least one run achieving the best-know value for the objective function, if possible) are marked in bold.

Results for Dataset 1
We solve the first dataset, comprised of 20 instances proposed by Osman & Christofides (1994), utilizing IBM's MIP solver CPLEX vs.12.6 and our CRO for CPMP.Each instance is solved to   optimality by CPLEX.Table 3 shows the optimal objective function value (opt) and respective execution time to solve it using CPLEX.Then, we use CRO to solve each instance 30 times.The main algorithm stops when MaxIterations or MaxIterationsWithoutImprovement is reached.It also stops if, at any given iteration, the optimal value for the objective function is reached.
In our experiments, CRO reaches the optimal value at least twice on all of the test instances; with the worst gap of 0.62% on instance cpmp11.The average of all gaps is of 0.089%.Notice that the constructive algorithm is able to generate an optimal solution for instances cpmp02 and cpmp04.Therefore, the average number of iterations is zero.
On average, CRO is 18.24s faster than CPLEX, with an average execution time of 2.62s versus 20.86s of CPLEX.CRO is slightly slower when solving instances 11, 14 and 17.However, the differences in execution time do not exceed 1.7s.
We also compare the results from our CRO for the CPMP with the ones reported, for the same dataset, by Osman & Christofides (1994), who proposed a metaheuristic for the CPMP that   but much faster than the other methods in the group.On average, it is 381.03sfaster than Fen-Cplex, the third fastest method, and 3055.81sfaster than VNS, the slowest one.Regarding optimality, SS achieved it in only 3 instances, being the worst in the group, whereas the other methods achieved optimality in 4 or more instances, thus matching or surpassing CRO.While the average of the best gaps achieved by CRO can be considered very satisfactory, in our opinion, it is higher than the other methods, especially the ones that used an MIP solver.
Despite not being the fastest or most accurate method in this comparison, we believe that our implementation of the CRO for CPMP may have a competitive advantage on applications that must run on less capable hardware.The reason is that it requires a much smaller memory footprint (typically less than 64 MB) and runs on a single CPU/Core, whereas most hybrid solutions that employ MIP solvers, such as IBM CPLEX, require large amounts of memory (typically more than 2GB) and multicore CPUs to perform well.The high cost of licensing an MIP solver may also be a limiting factor to such hybrid solutions.Furthermore, the average gap/processing time that can be obtained with CRO may be acceptable in many scenarios.

Results for Datasets 3, 4 and 5
The last three datasets we utilize for accuracy and speed comparisons was proposed by Stefanello et al. (2015) and is comprised of 3 data sets of 5 instances each, adapted from TSP-LIB, with 318 to 724 nodes and 5 to 200 medians.We solve all problems 10 times using CRO and compare the accuracy and execution times of our results with the ones obtained Stefanello et al. (2015), who also solved those 10 times using the matheuristic IRMA, which combines a model reduction heuristic and a MIP solver (CPLEX 12.3).We compare the average gaps and average execution times obtained by CRO with the ones obtained by Phase 2 and 3 of IRMA.Since most of these problems cannot be solved to optimality by MIP solvers in a timely fashion, we report in Table 7, their objective function values for the Best Know feasible Solution (BKS), regardless of being achieved by CRO or IRMA.Whenever optimal values are available, they are marked in bold.Table 7 shows the average gaps and average execution times for IRMA's Phase 2 and Phase3, as well as for CRO.It also shows, for CRO only, the percentage of times BKS is achieved, the average and minimum objectives attained, the standard deviation, the average execution time, the average number of decompositions, on-wall collisions, syntheses and intermolecular collisions.At last, we report in Table 7 the values we set for MaxIterations, MaxIterationsWithoutImprovement and PopSize.All other CRO parameters remain constant across all instances: PopSize = 10; KELossRate = 0.8, MoleColl = 0.1, InitialKE = 1,000,000, α = 10, β = 100,000, buffer = 0, MinMol = 2, MaxMol = 100 and Minλ-InterchangeIterations = 1.We set κ0 = 1 and κ = 1 on all instances, but u724 125 and u724 200, which have both κ0 and κ = 1 set to 0.1.A 1-opt optimization is used on all ineffective collisions.
Test run results show that our CRO for the CPMP performs better when solving instances where p (number of medians) is low, more specifically, ranging from 5 to 15 medians.Within this range, it is able to achieve BKS in, at least, 30 percent of the test runs, with average gaps below 0.08%.As p increases, results become less accurate and execution times increase significantly.This behavior can be explained by the nature of the λ-interchange mechanism we use, in which every customer of a given cluster is systematically relinked to all selected stations in the vicinity of the station it is currently assigned to, looking for possible improvements in the objective function value.In the worst case scenario, for a given λ there are p( p − 1)/2 different pairs of clusters to be examined.For the tested group of datasets, the worst average gap and execution time are, respectively, 3.56% and 3621.24 for instance ali 535 150.Although instance ali 535 005 has the best gap and average execution time, BKS is achieved during the constructive phase, thus not requiring any iterations from CRO's main loop.Therefore, the best results that required CRO are obtained with instance ali 318 005, also with no gap and 9.76s execution time.
In comparison with IRMA's Phase 2, the fastest but least accurate version of IRMA, CRO achieves same or better gaps in 6 out of 15 instances and it is faster in 5 of them.The average gap of all instances is 1.17%, thus 0.45% higher than IRMA's Phase 2 of 0.72%.Phase 3 is the most accurate IRMA and comprises Phase 2. Since not all instances required Phase 3, our comparison is restricted to the instances reported by Stefanello et al. (2015).CRO is faster in 5 out of 15 instances but achieved the same or worse gaps on all instances.On average, the gap percentage difference between CRO and Phase 3 is of 1.14%.

CONCLUSIONS
This paper presents an implementation of the CRO metaheuristic for solving the capacitated p-median problem (CPMP), which we call the CRO for the CPMP.To provide intensification of the solutions stored in molecules, we use a modified λ-interchange mechanism operator on all intermolecular collisions to generate new neighborhoods.To limit the number of iterations executed by the λ-interchange mechanism, we implement a Proximity List Set, which contains lists of nodes that are in the vicinity of each selected service station.To provide diversification we use a Distance Preserving Crossover (DPX) operator on syntheses and a Half-Total Change operator on decompositions.
Computational results on benchmark data sets of the literature demonstrate that our proposed solution can be competitive on a variety of implementations, as it is does not require complex and potentially expensive MIP solvers and runs on a single processor with low memory usage, while still providing results with accuracy and speed that may be acceptable on many real-life applications.
Future developments may include the development of a parallel version of the CRO for the CPMP, which may be able to tackle larger instances, above 1000 nodes.Additionally other, more efficient, operators can be developed to support molecular collisions, thus increasing the quality of results.
CRO has been used to address a broad range of problems in both discrete and continuous domains.Lam & Li (2010) proposed a solution for the quadratic assignment problem, whereas Xu et al. (2010) implemented a parallelized version of CRO for the same problem.James et al. (2011) employed CRO to train artificial neural networks.Lam et al. (2012) extended CRO to solve continuous problems.Xu et al. (2011b) utilized it for stock portfolio selection.CRO was also used for solving task scheduling problems in grid computing by Xu et al. (2011a).

Figure 4
Figure4illustrates the aforementioned shift and interchange processes for a 1-interchange mechanism (λ = 1).A shift process occurs from Figure4(a) to Figure4(b), where a customer i is shifted by the (1, 0) operator.As a result, customer j becomes the new median.Figure4(c) and Figure4(d), show the change in the clusters after customers i and j are interchanged, which also causes the medians to change on both clusters.
is a replacement of a subset C i ⊆ C i of size |C i | ≤ λ by another subset C j ⊆ C j of size |C j | ≤ λ to get two new clusters The λ-interchange generates new neighborhoods as follows: Let C i be a cluster comprised of a number of customers assigned to a station (median), ζ i .Given a solution S = {C 1 , . . ., C i , . . ., C j , . . ., C p } with ζ i as the median of cluster C i , a λ-interchange between two given clusters C i and C j