A simple multi-wave algorithm for the uncapacitated facility location problem

Abstract

The multi-wave algorithm (Glover, 2016) integrates tabu search and strategic oscillation utilizing repeated waves (nested iterations) of constructive search or neighborhood search. We propose a simple multi-wave algorithm for solving the Uncapacitated Facility Location Problem (UFLP) to minimize the combined costs of selecting facilities to be opened and of assigning each customer to an opened facility in order to meet the customers’ demands. The objective is to minimize the overall cost including the costs of opening facilities and the costs of allocations. Our experimental tests on a standard set of benchmarks for this widely-studied class of problems show that our algorithm outperforms all previous methods.

Keywords

discrete optimization / UFLP / multi-wave optimization / strategic oscillation / tabu search

Introduction

Recently, ^{Glover (2016)} introduced a new multi-wave algorithm (MWA) that embraces both (1) iterated constructive methods and (2) iterated neighborhood search methods (see Algorithm 1). An iterated constructive search starts from a null solution and employs moves that add elements until a complete solution is constructed, and an iterated neighborhood search starts from a given solution and employs improving moves to reach a local optimum.

In this paper, we focus on iterated improvement neighborhood search (which only accepts moves that are improving) as a natural analog of iterated constructive search.

To describe the method in overview, let N(x) = the set of solutions that are neighbors of x, and M(x) = the set of moves that lead to these neighbors; i.e. M(x) = {m = m(x) ∈ N(x)}. Let CL denote a candidate list of moves extracted from M(x). For problems with large neighborhoods we assume M(x) has been screened to reduce its size using a candidate list approach as described in ^{Glover and Laguna (1997)}.

In the following, we use the term boundary solution to refer to both a completely feasible solution obtained by a constructive algorithm and to a local optimum obtained by a neighborhood search algorithm. We denote a best solution by x*, which begins as a dummy solution with f(x*) = + ∞, and f(x) denotes the objective function to be minimized, as identified below the progression from a starting solution x° to a boundary solution x will be called a wave of the solution process. MaxPass and MaxWave refer to selected iteration limits. The reference to generating a null starting solution at the start of the “Pass loop” applies to constructive search. We assume MaxPass = 1 for this form of search, for reasons elaborated later.

Algorithm 1. Multi-Wave Algorithm (Overview)

1) x* = ∅ ; f(x*) = + ∞;

2) For Pass = 1 to MaxPass do

3) Generate a starting solution x°;

4) For Wave= 1 to MaxWave do

5) x = x°;

6) While x is not a boundary solution do x = m(x) where m is a move in CL;

7) x* = argmin{f(x), f(x*)};

8) Update relevant search parameters for the next Wave;

9) EndFor

10) Update relevant diversification parameters for the next Pass;

11) EndFor

12) Return x*

The Uncapacitated Facility Location Problem (UFLP).

Our following development focuses on the design and implementation of a simple multi-wave algorithm for solving the Uncapacitated Facility Location Problem (UFLP).

A variety of facility location problems arise in supply chain management and most of them are combinatorial in nature (^{Blum and Roli, 2003}; ^{Ardjmand et al., 2014}; ^{De Corte and Sörensen, 2015}; ^{Cerrone et al., 2017}; ^{De Armas et al., 2017}). The location problems have several industrial applications. For example, determining the number of facilities to be installed and their locations has become an important issue for manufacturing and industrial companies (^{Melo et al., 2009}; ^{Chen et al., 2014}). The determination of the best facilities to be implemented, by solving location problems, can significantly reduce the management cost and increase the performance of manufacturing and industrial companies (^{Cheung et al., 2001}; ^{Melo et al., 2009}; ^{Amin and Zhang, 2013}).

The Uncapacitated Facility Location Problem (UFLP) is one of the most widely studied discrete location problems from this class, and consists of choosing a subset of facilities to be opened from a set of potential candidates, and of assigning each customer to an opened facility in order to meet the customers’ demands. The objective is to minimize the overall cost which is the sum of the costs of opening facilities and allocating customers to them. An UFLP can be defined on a bipartite graph G = (P∪ Q, A) with p + q vertices and pq arcs where p = |P| is the number customers and q= |Q| is the number of candidate facilities. The cost of opening a facility is f_i and the cost of serving customer j from facility i is c_ij. The UFLP is to find a subset O⊆ Q of facilities to open to minimize the function

f\left(O\right)={\displaystyle {\sum }_{i\in O}{{\displaystyle f}}_i}+{\displaystyle {\sum }_{j\in P}\underset{i\in O}{\min }{{\displaystyle c}}_{i,j}}.

Since the only combinatorial part of the problem is the selection of open facilities, it is natural to represent a solution by a binary vector x∈ {0, 1}^q where x_i equals 1 if facility i is open and equals 0 otherwise. We use the notation Q¹ = Q¹(x) = {i∈ Q: x_i = 1} to represent the facilities that are opened in a solution x. Hence, the optimization problem of interest may be expressed in the form

Minimizef\left(x\right)={\displaystyle \sum _{i\in Q^{\mathrm{1}}\left(x\right)}{{\displaystyle f}}_i+}{\displaystyle \sum _{j\in P}\underset{i\in Q^{\mathrm{1}}\left(x\right)}{\min }}{{\displaystyle c}}_{i,j}:x\in {\left\{\mathrm{0},\mathrm{1}\right\}}^q-\left\{\mathrm{0}\right\}.

Note that all binary vectors in {\left\{\mathrm{0},\mathrm{1}\right\}}^q are feasible solutions except for 0.

In recent decades, many methods have been proposed for solving the UFLP, including heuristic methods, exact methods and hybrid methods. Several state-of-the-art surveys on location problems and their applications can be found in ^{Hale and Moberg (2003)}, ^{Klose and Drexl (2005)}, ^{ReVelle and Eiselt (2005)}. Among exact algorithms that have been developed to solve UFLP problems optimally are the branch and bound algorithm of ^{Khumawala (1972)}, a dual-based approach by ^{Erlenkotter (1978)} and its improved version by ^{Körkel (1989)}, a three-stage algorithm by ^{Galvão and Raggi (1989)} and the branch and PEG algorithms of ^{Goldengorin et al. (2003)}. In addition, ^{Posta et al. (2014)} recently introduced a cooperative exact method for solving UFLP integrating tabu search and branch and bound with the aim of finding an optimal solution.

A variety of methods have also been proposed based on Lagrangian relaxation. Within this class, ^{Beasley (1993)} developed several Lagrangian heuristics for solving location problems including the UFLP. Later, ^{Beltran-Royo et al. (2012)} and ^{Jörnsten and Klose (2016)} used Semi-Lagrangian relaxation to solve the UFLP drawing on the notion of Semi-Lagranian relaxation introduced in ^{Beltran et al. (2006)} for solving the p-median problem.

Of particular relevance to our current work, a wide range of algorithms have been developed for UFLP using a metaheuristic framework. ^{Basu et al. (2015)} present a review on the application of metaheuristics for solving discrete location problems including the UFLP, in which they focus on four methods: scatter search, tabu search, particle swarm optimization and genetic algorithms. Specifically addressing the UFLP, ^{Ghosh (2003)} studied several neighborhood search methods including tabu search and a complete local search with memory. ^{Resende and Werneck (2006)} proposed a multistart hybrid algorithm (called HYBRID) for UFLP which is based on the algorithm introduced in ^{Resende and Werneck (2004)} for solving the p-median problem, and which combined efficient procedures such as path-relinking (used in tabu search and scatter search) together with notions from genetic algorithms. In addition to these approaches, many different variants of tabu search have been successfully applied to UFLP (^{Al-Sultan and Al-Fawzan, 1999}; ^{Michel and Van Hentenryck, 2004}; ^{Sun, 2005}; ^{Sun, 2006}) and Genetic Algorithms (^{Kratica et al., 2001}; ^{Homberger and Gehring, 2008}; ^{Tohyama et al., 2011}). Other types of approaches include those of ^{Greistorfer and Rego (2006)}) who used a simple filter-and-fan method to deal with UFLP and ^{Yigit et al. (2006)} who proposed an evolutionary simulated annealing algorithm for solving UFLP. The study of ^{Arostegui et al. (2006)} presented a comparison study between simulated annealing, tabu search and genetic algorithms for location problems, in which tabu search obtained particularly promising results. A multi-population particle swarm optimization (MPSO) algorithm was applied by ^{Wang et al. (2008)} which included a parallel implementation. ^{Tseng and Wu (2009a)} developed a novel hybrid algorithm called multiple trajectory search (MTS) and the same authors also presented a multi-start drop-add-swap heuristic to deal with UFLP in ^{Tseng and Wu (2009b)}. Other multi-start approaches for the UFLP include those of ^{Cura (2010)} and ^{Ardjmand et al. (2014)}. The latter multi-start heuristic was dubbed discrete unconscious search and obtained highly promising results. Employing a synthesis of several of the preceding ideas, De Armas et al. 2017 have proposed a method for the stochastic instance of UFLP (SUFLP), using a simheuristic (simulation optimization) approach that first embeds a fast heuristic inside the metaheuristic framework of iterated local search (ILS) for the deterministic version of UFLP and then integrates this procedure with Monte Carlo simulation techniques. Additional recent investigations of UFLP have been carried out by ^{Albareda-Sambola et al. (2017)}, ^{Atta et al. (2018)}, ^{Galli et al. (2018)}, ^{Sahman et al. (2017)}, ^{Tsuya et al. (2017)}. There also exist several recent studies of variants and extensions of UFLP, including ^{Akbaripour et al. (2017)}, ^{Chalupa and Nielsen (2018)}, ^{Han et al. (2018)}, ^{Jiang et al. (2018)}, ^{Pearce and Forbes (2018)}, ^{Todosijević et al. (2017)}, ^{An and Svensson (2017)}. Also quite recently, ^{Ortiz-Astorquiza et al. (2017a}, ^2017b) respectively provide formulations and approximation algorithms for Multilevel Uncapacitated Facility Location and a comprehensive review on multi-level facility location problems that extend several classical facility location problems.

The remainder of the paper is organized as follows. The next section presents the principles that underlie the Multi-Wave Algorithm and Section 3 exposes the general design of the algorithm. Our comparative experimental study is presented in Section 4 and concluding observations are provided in Section 5.

Multi-wave metaheuristic principles

The Multi-Wave Algorithm is based on combining tabu search and strategic oscillation to exploit three chief components: (1) a candidate list (CL) for exploring the neighborhood of the current solution; (2) an active move record (AMR) to capture the conditional effects of sequential decisions; (3) intensification and diversification strategies founded on the notion of persistent attractiveness. The following subsections present a brief description of each component, but for our present simple version of MWA we only make use of components (1) and (2). As shown in Section 4, even without component (3), our method outperforms all previous methods on the set of benchmark problems standardly used for comparative testing in the literature.

Candidate list

Let \Delta \left(m\right) to be the evaluation of the move m, i.e. the change in the objective function that results when performing m on x. We focus attention on a candidate list that contains improving moves with higher evaluations, hence CL = {m₁, m₂, …, m_s} so that D(m₁)≥...≥D(m_s)≥max{D(m): m∈ M(x) - CL}. Often the size s = |CL| is determined by a threshold T (^{Glover, 1989}, ¹⁹⁹⁵), i.e. CL = {m∈ M(x): D(m)≥T} where in the present work, the value of T is chosen as

(1)T={{\displaystyle \mathrm{\Delta }}}_{\mathit{min}}+\lambda \left({{\displaystyle \mathrm{\Delta }}}_{me\mathit{a}n}-{{\displaystyle \mathrm{\Delta }}}_{\mathit{min}}\right),

(2)T={{\displaystyle \mathrm{\Delta }}}_{mean}+\lambda \left({{\displaystyle \mathrm{\Delta }}}_{\mathit{max}}-{{\displaystyle \mathrm{\Delta }}}_{mean}\right),

and where D_min, D_mean and D_max refer to the min, mean and max evaluations D(m) over M(x), and λ∈ [0, 1] is a parameter. Note that a threshold T is also used in GRASP ^{Resende and Ribeiro (2003)} which is computed without reference to D_mean by defining T = D_min + α(D_max−D_min), relative to a parameter α∈ [0, 1]. If the midpoint of the preceding representation, for α = 0.5, is taken to be an approximation to D_mean, then a rough correspondence between this latter candidate list and those based on D_mean occurs by specifying α = (λ + 1)/2 (hence α∈ [.5, 1]) to associate it with (2).

We assume that the evaluations have been scaled and translated, if necessary, so that they are all positive. In the present work, we choose a move m from CL based on the principle of probabilistic tabu search (P-TS) (^{Glover, 1989}) by constructing a sample of size 1 from the uniform distribution with probability Pr(m) given by

(3)\mathit{Pr}\left(m\right)={\displaystyle \frac{{\displaystyle \mathrm{\Delta }}{\left(m\right)}^{\beta }}{{\displaystyle {\sum }_{m^{\prime }\in CL}{\displaystyle \mathrm{\Delta }}{\left(m^{\prime }\right)}^{\beta }}}.}

The exponent β can be selected greater than 1 to emphasize the differences between evaluations, or selected less than 1 to de-emphasize these differences. (Another early form of probabilistic tabu search generates probabilities by reference to the exponential function e^r for chosen values of the exponent r (^{Xu et al., 1997}).) The completely random choice that results for β = 0 effectively corresponds to the approach subsequently proposed by GRASP for its candidate list. More recently, probabilistic variants of GRASP have been introduced in ^{Wang et al. (2013)} and ^{Grasas et al. (2017)} which adopt a number of ideas from probabilistic tabu search.

In the context of UFLP, let x and x′ represent two binary solutions where x′ is obtained from x by a 1-flip move which changes the value of a single variable from {{\displaystyle x}}_k to {{\displaystyle x}}_k^{\text{'}}=\mathrm{1}-{{\displaystyle x}}_k, so we have x^{\prime }=x+\left(\mathrm{1}-{{\displaystyle \mathrm{2}x}}_k\right)e^k, where e^k is unit vector with all components zero except the k^th component. To evaluate a 1-flip move efficiently, for each customer j, we maintain the closest facility and the second closest facility, i.e.,

{{\displaystyle i}}_j^{\mathrm{1}}={{\displaystyle i}}_j^{\mathrm{1}}\left(x\right)=\mathit{arg}\mathit{min}\left\{{{\displaystyle c}}_{i,j}:i\in Q^{\mathrm{1}}\left(x\right)\right\},

{{\displaystyle i}}_j^{\mathrm{2}}={{\displaystyle i}}_j^{\mathrm{2}}\left(x\right)=\mathit{arg}\mathit{min}\left\{{{\displaystyle c}}_{i,j}:i\in Q^{\mathrm{1}}\left(x\right)-\left\{{{\displaystyle i}}_j^{\mathrm{1}}\left(x\right)\right\}\right\}.

Then the objective function can be rewritten as f\left(x\right)={\displaystyle {\sum }_{i\in Q^{\mathrm{1}}}{{\displaystyle f}}_i}+{\displaystyle {\sum }_{j\in P}{{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{1}}}j}. Then the objective function change produced by flipping x_k is given by {{\displaystyle \mathrm{\Delta }}}_k=f\left(x\right)-f\left(x^{\prime }\right). The objective function value after closing facility k\in Q^{\mathrm{1}} is given by

f\left(x^{\text{'}}\right)={\displaystyle \sum _{i\in Q^{\mathrm{1}}}{{\displaystyle f}}_i}-{{\displaystyle f}}_k+{\displaystyle \sum _{j\in P-P_k^{=}}{{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{1}}}}j+{\displaystyle \sum _{j\in P_k^{=}}{{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{2}}}}j,

where {{\displaystyle P}}_k^{=} is the set of customers assigned to k, i.e., {{\displaystyle P}}_k^{=}=\left\{j\in P:{{\displaystyle i}}_j^{\mathrm{1}}=k\right\}. The change {{\displaystyle \mathrm{\Delta }}}_k^{-}of closing facility k is given by

(4){{\displaystyle \mathrm{\Delta }}}_k^{-}={{\displaystyle f}}_k-{\displaystyle {\sum }_{j\in {{\displaystyle P}}_k}=\left({{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{2}}}j-{{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{1}}}j\right)}.

The objective function value after opening facility k\in Q^{\mathrm{0}}=Q^{\mathrm{0}}\left(x\right)=\left\{i\in Q:{{\displaystyle x}}_i=\mathrm{0}\right\} is given by

f\left(x^{\text{'}}\right)={\displaystyle {\sum }_{i\in Q^{\mathrm{1}}}{{\displaystyle f}}_i+{{\displaystyle f}}_k}+{\displaystyle {\sum }_{j\in P}\mathit{min}\left({{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{2}}j}-{{\displaystyle c}}_{kj}\right)}.

The change {{\displaystyle \Delta }}_k^{+} of opening facility k is given by

(5){{\displaystyle \mathrm{\Delta }}}_k^{+}=-{{\displaystyle f}}_k+{\displaystyle {\sum }_{j\in P}\mathit{max}\left\{\mathrm{0},{{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{1}}j}-{{\displaystyle c}}_{kj}\right\}}=-{{\displaystyle f}}_k+{\displaystyle {\sum }_{j\in {{\displaystyle P}}_k-{{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{1}}j}-{{\displaystyle c}}_{kj}}},

where . Hence, the change {{\displaystyle \Delta }}_k of flipping facility k can be expressed as

{{\displaystyle \mathrm{\Delta }}}_k=\{\begin{array}{lll}{{\displaystyle \mathrm{\Delta }}}_k^{-}\hfill & if\hfill & {{\displaystyle x}}_k=\mathrm{1}\hfill \\ {{\displaystyle \mathrm{\Delta }}}_k^{+}\hfill & if\hfill & {{\displaystyle x}}_k=\mathrm{0}\hfill \end{array}.

^{Michel and Van Hentenryck (2004)} and ^{Sun (2006)} provide an efficient procedure in O(p log q) time to update the evaluation {{\displaystyle \Delta }}_k for facility and warehouse location problems after each flip move.

Now we consider a swap move. Let x and x′ represent two binary solutions where x′ is obtained from x by swapping the value of two variables {{\displaystyle x}}_k and {{\displaystyle x}}_h with h\in Q^{\mathrm{1}}\left(x\right)and h\in Q^{\mathrm{0}}\left(x\right), so we have Q^{\mathrm{1}}\left(x^{\text{'}}\right)=Q^{\mathrm{1}}\left(x\right)-k+h. Then the objective function f\left(x^{\text{'}}\right)={\displaystyle {\sum }_{i\in Q^{\mathrm{1}}\left(x\right)}{{\displaystyle f}}_i-{{\displaystyle f}}_k+{{\displaystyle f}}_h}+{\displaystyle {\sum }_{j\in P}\underset{i\in Q^{\mathrm{1}}\left(x\right)-k+h}{\min }}{{\displaystyle c}}_{i,j}. Note that

Let , and P^{\mathrm{3}}=P-P^{\mathrm{1}}-P^{\mathrm{2}} be a partition of P. Then the objective function change produced by swapping x_h and x_k is given by {{\displaystyle \mathrm{\Delta }}}_{hk}=f\left(x\right)-f\left(x^{\text{'}}\right). More precisely, we have

(6){{\displaystyle \mathrm{\Delta }}}_{hk}={{\displaystyle f}}_k-{{\displaystyle f}}_h+{\displaystyle {\sum }_{j\in P}\left(\underset{i\in Q^{\mathrm{1}}\left(x\right)}{\min }{{\displaystyle c}}_{ij}-\underset{i\in Q^{\mathrm{1}}\left(x\right)-k+h}{\min }{{\displaystyle c}}_{ij}\right)}

This last evaluation can be simplified as follows

(7){{\displaystyle \mathrm{\Delta }}}_{hk}={{\displaystyle f}}_k-{{\displaystyle f}}_h+{\displaystyle {\sum }_{j\in P^{\mathrm{2}}}\left({{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{1}}j}-{{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{2}}j}\right)}+{\displaystyle {\sum }_{j\in P^{\mathrm{3}}}\left({{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{1}}j}-{{\displaystyle c}}_{hj}\right)}.

According to ^{Resende and Werneck (2007)}, the objective function change can be rewritten in the following form

(8){{\displaystyle \mathrm{\Delta }}}_{hk}={{\displaystyle \mathrm{\Delta }}}_k^{+}-{{\displaystyle \mathrm{\Delta }}}_h^{-}+{{\displaystyle \mathrm{\Delta }}}_{hk}^{\mp },

where

{{\displaystyle \mathrm{\Delta }}}_{hk}^{\mp }={\displaystyle \sum _{j\in P^{\mathrm{2}}}{{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{2}}j}-\max \left({{\displaystyle c}}_{{{\displaystyle i}}_j^{\mathrm{1}}j}-{{\displaystyle c}}_{kj}\right),}

and {{\displaystyle \mathrm{\Delta }}}_k^{+} is called gain(k) which corresponds to the decrease in solution value due to the insertion of facility k (with no associated removal). Similarly, {{\displaystyle \mathrm{\Delta }}}_h^{-} is called loss(h) corresponding to the increase in solution value due to the removal of facility h (with no associated insertion), while {{\displaystyle \mathrm{\Delta }}}_{hk}^{\mp } is called extra(k, h) corresponding to a positive correction term that accounts for the fact that the insertion of facility k and the removal of facility h may not be independent (a customer previously assigned to facility h may be reassigned to facility k). Note that i) {{\displaystyle \mathrm{\Delta }}}_{hk}^{\mp } is always positive or equal to 0; ii) If {{\displaystyle \mathrm{\Delta }}}_h^{-}is negative, then there is an improvement when closing h; iii) If {{\displaystyle \mathrm{\Delta }}}_k^{+} is positive, then there is an improvement when opening k.

In our implementation, the candidate list contains Drop, Add and Swap moves. To keep the size of the candidate list moderate, only the Swap(k, h) moves that satisfy the condition {{\displaystyle \mathrm{\Delta }}}_{hk}>\mathrm{0} and {{\displaystyle \mathrm{\Delta }}}_{hk}^{\mp }>{{\displaystyle \mathrm{\Delta }}}_h^{-} are candidates to belong to CL(x). This condition assures by Eq. (8) that we only consider Swap moves that dominate Add or Drop moves.

Conditional effect of sequential decisions and active record move

Starting from an initial solution, the algorithm repeatedly performs improving moves until a boundary solution is reached. In this sequence of decisions, the moves applied in early steps of the MWA affect the evaluation and the choice of the decision in the current step because the information underlying the current decision depends on these earlier moves. (This is the so-called conditional effect of sequential decisions.) Furthermore, the evaluation of current decisions is able to take advantage of more refined information than is available to early decisions, conditional upon decisions already made, due to the fact that the problem is now reduced, with the implication that later decisions are likely to be better in a conditional sense than earlier decisions. Consequently, as observed in ^{Glover (2000}, ²⁰¹⁶), it becomes advantageous to re-examine earlier decisions (which are likely to be bad) at intermediate stages during the search process in order to amend them and thus create paths to higher quality boundary solutions.

In its general design, the Multi-Wave Algorithm exploits this notion of the conditional effect of sequential decisions using an Active Move Record (AMR) whose form is described below. Specifically, the algorithm starts with a vertical phase by performing constructive or improving moves, and then after a given number of iterations the moves (decisions) executed earlier are re-examined a horizontal phase to see if additional improvement results by replacing them with other moves that now are evaluated to be better. An alternation between vertical and horizontal phases which are triggered by parameters of the search continues until a boundary solution is reached, whereon updates are made and the process begins again.

We refer to both constructive and improving moves as forward moves, and refer to destructive moves and moves that complement (or “cancel”) previous improving moves as reverse moves. For a constructive algorithm, the operation of dropping and reversing a move m is often very simple, consisting only of removing the move m from AMR. We will understand the reference to a reverse move \overline{m} in this case to involve just the dropping operation. We also refer both to the (arbitrary) beginning solution that launches an improving phase of neighborhood search and to the null construction that launches a constructive phase, as the initial solution.

More precisely, the Active Move Record identifies the forward moves selected so far that have not been cancelled by reverse moves. It is convenient to represent this record by AMR = (m₁, m₂, …, m_r), where r = |AMR|, and where the indexing sequence indicates that move m_i was executed before move m_i+1 for i∈ [1, r - 1]. It should be noted that M(x) and CL can include moves that reverse moves in AMR.

At any stage, the effect of applying the moves in AMR to the initial solution yields a current solution x and in the case of a constructive algorithm, we assume that the objective function value f(x) can be evaluated for a partial solution as well as a complete solution. We first sketch the general form of the Multi-Wave algorithm and then describe the ways that AMR and CL are used within it.

General design of the multi-wave algorithm

Each wave is divided into alternating vertical and horizontal components, where vertical steps consist of forward moves that enlarge the AMR and horizontal steps consist of combined forward and reverse moves designed to leave the cardinality of AMR unchanged. Here we handle the horizontal phase by the simpler approach of dropping moves from AMR by reversing them, and then adding an equal number of forward moves to AMR. The key steps for carrying this out involve: (1) deciding when a horizontal phase is to be launched, (2) selecting the number of moves to be executed during this phase, and (3) identifying the particular moves drawn from AMR to drop and the associated forward moves from CL that make up the phase.

We index the successive intervention points by k = 1, 2, …, k*. The parameter v_k denotes the number of vertical steps (forward moves) executed before launching intervention k, which therefore occurs when |AMR| = v_k. Finally, h_k denotes the number of horizontal steps executed during this intervention and d_k denotes the number of moves dropped on each of these steps. The values of these parameters are interrelated and we give simple rules to determine them. The condition for reaching an intervention point and executing a horizontal phase may be expressed as follows.

If |AMR| = v_k, perform h_k horizontal steps each consisting of dropping (reversing) d_k moves from the start of AMR and sequentially adding d_k new moves to AMR from CL (where CL has access to the moves dropped from AMR). This sequential selection implies that the move evaluations may change from one step to the next. If no forward move exists to replace a given reverse move (as where no improving moves currently exist), a boundary solution is reached and the wave terminates.

In short, all steps of adding moves to AMR during the horizontal phase have exactly the same form as the steps of choosing and executing forward moves during the vertical phase.

Algorithm 2 describes in detail the main steps used in the proposed Multi-Wave Algorithm.

Figure 1 shows the evolution of the AMR in the Vertical Phase when which consists of adding moves to AMR until it reaches the size r={{\displaystyle v}}_k. Figure 2 illustrates the steps of the Horizontal Phase that performs d_k Drop / Add moves h_k times. Note that the Concluding Horizontal Phase can be also illustrated by Fig. 2 where the process iterates h₀ times by dropping d₀ moves from AMR and adding moves until reaching a boundary solution.

Algorithm 2. Multi-Wave Algorithm

1) x = ∅ ; x* = x; f(x*) = ∞;

2) For Pass = 1 to MaxPas do

3) Generate a starting solution x° (null, randomly or by post-analysis diversification);

4) For Wave= 1 to MaxWave do

5) x = x°; AMR = ∅; k = 1;

6) While x is not a boundary solution do

7) If |AMR| = v_k then // Horizontal Phase

8) For h = 1 to h_k do

9) Reverse the first d_k moves and remove them from AMR

10) Perform d_k forward moves and add them to AMR

11) EndFor

12) k← k + 1;

13) Else // Vertical Phase

14) Perform one forward move and add it to AMR

15) EndIf

16) EndWhile

17) x* = argmin{f(x), f(x*)};

18) For h = 1 to h₀ do // Concluding Horizontal Phase

19) Reverse the first d₀ moves and remove them from AMR

20) While x is not a boundary solution do

21) Perform a forward move and add it to AMR

22) EndWhile

23) x* = argmin{f(x), f(x*)};

24) EndFor

25) Update v_k, h_k and d_k for all k;

26) EndFor

27) EndFor

28) Returnx*

A variety of rules are possible for determining the parameters of the Multi-Wave Algorithm, and particularly those of the horizontal phase. We focus here on rules that are easy to execute provided in subsequent sections.

Determining the intervention parameter values v_k and k*

The v_k values which identify the values of |AMR| at which interventions occur can be conveniently established by using an estimate r_e of the final |AMR| value upon reaching a boundary solution. Such an estimate can be determined by executing a preliminary wave without any interventions (as insured by making v₁ large) and setting r_e = |AMR| when the wave terminates. This r_e estimate can be updated after each subsequent wave to equal the final |AMR| value obtained during this wave or the average of the final |AMR| values obtained so far.

Given r_e, a desired spacing between successive interventions can be selected based on either:

(1) The separation δ_k between successive values v_k at which successive interventions occur, which yields v_k = v_k-1 + δ_k for k = 0 to k*, where by convention v₀ = 0. The value δ_k should satisfy δ_k≥d_k + 1. In the simplest case, we select a constant value δ = δ_k≥2, which therefore determines k* = ⌊r_e / δ⌋. This approach is useful for densely populating the values v_k over the estimate r_e, hence is suited to choosing k* relatively large (in comparison to r_e).

(2) The chosen total number of interventions k*. In the simplest case, we divide the interval [1, r_e] into equal parts and we consider values of k* ranging from 1 to a selected upper limit (5 for small value), then for any value of k*, let v_k = [k r_e/(k* + 1)] or v_k = [(k + 1)r_e/(k* + 2)] for each k = 1 to k*. This approach is useful for sparsely populating the values v_k over the estimate r_e, and it is suited to choosing k* relatively small.

Note that with the first technique (1), the value of d determines the total number of interventions k*. We have tested the simplest version where δ_k = δ is a constant value, but this simplest rule worked well only for some instances of UFLP. Hence we implement a new rule where δ_k increases after each η waves η∈ {1, …, 4} starting with small value of δ₁∈ {2, 3} (i.e., δ_k+h = δ_k + τ with τ∈ {1, …, 6}). To implement the second rule (2), based on the estimation of r_e, we keep a record at each pass of the average number of moves performed up to the current wave_. Accompanying this, we increase the number of interventions k* by 1 at each wave (k* = k* + 1), starting with a small value (in our implementation, k* = 5). Hence, we have v_k = v_k−1 + d where \delta =\frac{{{\displaystyle \tilde{r}}}_e}{k*}.

Determining the number of drop moves d_k

The value d_k identifying the number of moves to be dropped during the kth intervention point should be chosen, as previously noted, so that d_k<δ_k = v_k- v_k-1. Again, an equal spacing option gives a possible starting point, giving each d_k a value d_k = d>1 for all k. Alternatively, d_k may be selected randomly at each intervention point to lie in the interval [1, d]. Another option is to let d_k be larger on the first step of the horizontal phase and then decrease it to a smaller value on subsequent steps. In the latter case, two alternatives of interest are to decrease d_k immediately to 1 and to decrease it by a single unit at each step until it reaches 1. In our implementation for UFLP we determined the number of drop moves by setting d_k = 2 for all k, adopting the equal spacing option because of its simplicity and effectiveness.

Determining the number of horizontal steps h_k

We determine the number of horizontal steps h_k as a function of the number of dropped moves d_k. To assure that every move in AMR is dropped on one of these steps performed on the kth intervention, we select h_k = ⌈|AMR|/ d_k⌉ in our current implementation, where |AMR| refers here to the size of the active move record when the kth intervention is launched. In particular, each time a move or a block of moves is selected from AMR: (a) the moves dropped from AMR are made available to become members of CL; (b) the rule of probabilistic tabu search is employed for successively choosing moves from CL to be executed; and (c) M(x) and hence CL are updated after each of these moves. It is reasonable to expect that the horizontal interventions will contribute to an intensification/diversification tradeoff in a way that will enable the P-TS choice rules to focus primarily on intensification, as by using relatively large λ and b values. If λ is chosen sufficiently large, e.g., in the interval [.98, 1.0], then the value of b will be largely irrelevant. In our implementation, at each pass, the value of b starts from 1.5 and ends at 0.5 and at each wave decreases by stepwize 1/MaxWave, while the value of λ starts from 0.8 and ends at 0 and at each wave decreases by stepwize 0.8/MaxWave. Starting from large λ and b values allows intensification of the search while slowly reducing those values at each pass permits an intensification/diversification tradeoff.

Concluding horizontal phase

The goal behind the concluding horizontal phase is to iron out imperfections produced after the last intervention. This procedure is applied to the AMR produced by the current wave. The difference between the concluding horizontal phase and the horizontal phase performed during the wave is: after dropping d₀ moves from AMR and reversing them, the concluding horizontal phase performs forward moves until a boundary solution is reached, rather than stopping after executing d₀ forward moves, as done in the horizontal phase of the wave. From this fact, the initial values of h₀ and d₀ are chosen to depend on the value |AMR|. Choosing a large value of d₀ can make a significant change in the solution x and leads us to start a new search from scratch, as in a simple multi-start search, and this will be done h₀ times. On the other hand, choosing a small value of d₀ allows us to iron out early imperfections, but if |AMR| is large then h₀ will receive a large value, and this will negatively affect the search procedure and its CPU-time. To achieve the goal of the concluding horizontal phase without inhibiting the effectiveness of the search process, we fix h₀ and d₀ by setting d₀ = 2 and h₀ = min(⌈|AMR|/d₀⌉, \overline{h}), where \overline{h} is the maximum number of iterations we imposed on h₀. Consequently, if ⌈|AMR|/d₀⌉>\overline{h} the procedure will not examine all the moves in AMR and in such a case only early moves in AMR are reversed. We are motivated to use this technique based on the notion of the conditional effect of sequential decisions.

Computational results

Our simple version of the multi-wave algorithm was coded in JAVA and all experiments were conducted on a computer with Intel Xeon E3-1505M CPU @ 2.80 GHz and 16 GB RAM. To test the performance of our MWA, we use four of the most widely referenced and challenging benchmark data sets of UFLP from the literature including: The benchmark of ^{Ghosh (2003)} (GHOSH), the data set of ^{Kratica et al. (2001)} (M*), the large-scale data set used in ^{Barahona and Chudak (1999)} (MED) and the test problems from the OR-Library of ^{Beasley (1990)} (ORLIB). These data sets can be uploaded from http://resources.mpi-inf.mpg.de/departments/d1/projects/benchmarks/UflLib/. In addition, the MWA is compared against nine of the most efficient methods developed for the UFLP. These state-of-the-art methods are:

1) LAG: The Lagrange-method from ^{Sun (2006)} which is implemented in ^{Sun (2005)}.

2) CLM: The Complete Local search with Memory of ^{Ghosh (2003)}.

3) GTS: The Tabu Search of Ghosh (2003).

4) HYB: The Hybrid multistart method of ^{Resende and Werneck (2006)}.

5) TS: The Tabu Search of ^{Sun (2006)}.

6) MDAS: The Multi-start Drop Add Swap heuristic of ^{Tseng and Wu (2009b)}.

7) MTS: Multiple Trajectory Search of ^{Tseng and Wu (2009a)}.

8) US: The discrete Unconscious Search of ^{Ardjmand et al. (2014)}.

9) ILS: The Iterated Local Search method of ^{De Armas et al. (2017)}.

Our experiments are divided into two parts. In the first part, we present a comparison between the proposed MWA and two multi-start algorithms we have implemented with the aim to identify performance differences between the MWA and classical multi-start approaches. In the second part of our experimentation, the MWA results are compared to the outcomes from the 9 state-of-the-art algorithms identified above. In all experiments, we report CPU-times for the MWA in two columns: Column All which refers to the full CPU-time of the execution and column Find which refers to the CPU-time needed to find the best solution in the execution. In what follows, first we give a brief description of the benchmark data sets used, then we present the two parts of the experiments.

Benchmarks for UFLP

We use the following four benchmark data sets to evaluate the MWA:

• GHOSH: In this data set from ^{Ghosh (2003)}, the number of customers p and the number of facilities q are set to be (p, q) ∈ {(250, 250), (500, 500), (750, 750)} with 30 problem instances of each size. These instances are divided into three groups A, B and C in which the facility setup costs f_i, are chosen uniformly from intervals [100, 200], [1000, 2000] and [10000, 20000] respectively for low, medium and high fixed cost values. In all instances, the assignment cost c_ij is chosen at random from [1000, 2000]. In addition, each set of test problems of the same size and from the same group are further divided into two categories, where in 5 instances c_ij is symmetric (i.e., c_ij = c_ji) and in the other 5 instances c_ij is asymmetric.

• MED: This benchmark contains 18 large-scale problems, originally introduced for the p-median problem by ^{Ahn et al. (1988)} and used later by ^{Barahona and Chudak (1999)} in the context of UFLP. We name each instance as q-y where the number of facilities q∈ {500, 1000, 1500, 2000, 2500, 3000} (each facility is a candidate for opening and represents a customer at the same time) and the parameter y∈ {10, 100, 1000}. Each combination of q and y is included to yield a total of 18 instances. The parameter y is used to calculate the setup cost {{\displaystyle f}}_i=\frac{\sqrt{q}}{y}, which is the same for all facilities.

• M*: The third benchmark is the data set of ^{Kratica et al. (2001)}. According to the authors, these problems are generated to closely approximate real-life cases. In these problems the number of customers is equal to the number of facilities and varies from 100 to 2000 (i.e., p = q∈{100, 200, 300, 500, 1000, 2000}).

• ORLIB: The fourth test problem set is taken from the OR-Library of ^{Beasley (1990)}. These instances were initially proposed for the capacitated facility location problem and then used for the uncapacitated facility location problem case by ignoring the capacity of facilities. This data set contains 15 instances, in which the numbers of facilities and customers vary from ((q, p) = (16, 50) to (q, p) = (100, 1000)).

Comparison between basic multi-start and the MWA

In this subsection, we present a comparison between two basic multi-start methods and the MWA. The aim behind this comparison is to identify the contribution provided by techniques used in our simple MWA that differentiate it from multi-start methods: (1) the use of the AMR and the division of a wave into horizontal and vertical phases in order to exploit the conditional effects of the sequential decisions, and (2) the manner in which we exploit the candidate list in order to select the most appropriate moves to be applied at each step of the algorithm. To achieve this aim, we compare the results returned by the MWA to those obtained by the two basic multi-start methods using the same CPU times.

The basic multi-start algorithms we implemented may be viewed as simple instances of the GRASP approach and consist of two sequential phases: (i) a Greedy Randomized Procedure (GRP) in which a new starting solution is generated, and (ii) a local search which improves the solution returned from the first phase. These steps are performed until the time limit stopping criterion is met. In traditional GRASP methods the GRP procedures are used to construct complete and/or feasible solutions but in our case, a UFLP-solution with only one facility is a complete and feasible solution. In fact, in our multi-start methods, the GRP procedure is used to find starting solutions with good tradeoff between the assignment cost and facilities opening cost. To do this, the GRP starts from a solution with one facility (i.e., with low opening cost even that generated a large transportation cost). To establish a good tradeoff between the opening cost and the assignment cost (that minimizes the overall cost) we search for the best number of facilities to be opened by starting with one facility and continue to add facilities using the greedy randomized technique until no improving add moves are available. Thus, the method performs improving add moves which both improve the solution quality and converge to a good tradeoff between the assignment and opening costs. After that, the method performs a local search procedure in which the Add, Drop and Swap moves are considered.

In this study, the two versions of the multi-start approach are examined: MultiS-FI, in which the local search adopts the first improvement choice strategy and MultiS-BI, in which the local search adopts the best improvement choice strategy. Table 1 shows the results of the three algorithms MultiS-FI, MultiS-BI and MWA: the first column presents the instances, then in columns from two to seven, for the two multi-start methods, we present the average solution cost, the best solution cost and the average number of iterations out of 10 runs. In columns eight and nine we show the MWA results, and finally, in column, ten we present the CPU times of the methods. These results are obtained as follows: First, we run the MWA algorithm on these instances by setting the MaxPas= 10 and MaxWave= 10 (so 100 waves are performed). Then we run the two multi-start methods using the same amount of CPU-times consumed by the MWA.

As shown in Table 1, the MWA algorithm clearly proves more effective than the two multi-start methods. The MWA obtains higher quality best solutions for all MED instances and in most instances, the average solution cost obtained by the MWA is better than the best solution cost obtained by the multi-start methods. In addition, we observe that the average number of iterations for MultiS-FI equals 1170 and for MultiS-BI equals 928. So, on average, around 1000 solutions are investigated by these methods in one execution, whereas the 100 solutions generated by executing only 100 waves of the MWA turn out to include much better outcomes. Hence, we conclude that the new techniques used in our simple MWA method allow it to uncover better regions in the search space where high-quality solutions are available.

Comparison with state-of-the-art algorithms

We now present the results of the MWA on each benchmark and compare these results to those obtained by the most effective methods from the literature.

Results on GHOSH test problems

As previously noted, this benchmark contains hard test problems. We compare the results of the MWA on these problems to the results of: TS, CLM, GTS, HYB, MTS, US and MDAS. Table 2 presents the solution costs (objective function values) obtained by each method and Table 3 presents the CPU-times of each method. The results for TS, GTS, HYB, MTS, CLM and US are taken from ^{Sun (2006)} and ^{Ardjmand et al. (2014)}. The results of MDAS are not reported in ^{Ardjmand et al. (2014)} and we take them directly from the original paper of ^{Tseng and Wu (2009b)}.

As shown in Table 2, the MWA was able to reach all best-known solutions (which were found by at least one of the most effective methods in the literature) and to obtain a new best known solution for the instance 750-asym-A. By contrast, none of the previously proposed methods was able to find all the best solutions previously known, disregarding the new best solution we have found. Table 3 additionally shows that our simple MWA has a competitive computing time compared to the other methods. The MWA performs particularly well on the type C problems where it is able to reach the best solution of two large-scale instances in less than 2 s (750-sym-C and 750-asym-C) and obtain best solutions for two of the medium size instances in less than 1 s (500-sym-C and 500-asym-C). In addition, on average, our algorithm needs only half the CPU-time required by the fastest method to obtain its best results (comparing the averages of these CPU-times).

Results on MED test problems

The second comparison is made on the MED test problems. In Table 4, we compare the four methods that have results on this data set: ILS, MDAS, HYB and MWA. In the first and the second columns we present the instance name and the optimal solution found by the solver Gurobi, then for each method, we report the average solutions cost and the average CPU-times out of 30 runs for ILS and HYB, 50 runs for MDAS and 10 runs for the MWA. The results of HYB and ILS are taken from ^{Resende and Werneck (2006)} and ^{De Armas et al. (2017)} and the results of MDAS are taken from ^{Tseng and Wu (2009b)}. Note that the results for MDAS are not reported in ^{De Armas et al. (2017)}.

Comparing MWA to ILS and MDAS, the MWA obtains better solutions for all problem instances with the exception of one problem in which MDAS and MWA obtain the same optimal solutions over all executions (500-10). Compared to ILS, the MWA method finds better solutions to all problem instances and requires smaller average CPU-time (778,23 versus 931,41). Although the average CPU-time for MDAS is slightly smaller than the average CPU-time for MWA, we observe that in all instances where q∈ {500, 1000, 1500, 2000, 2500, 3000}-y = 10 the CPU-times of MWA are better than those of MDAS. Interestingly, for the indicated problems from the MED test set, in contrast to the problem from the GHOSH test set, the HYB method is much closer to matching the MWA performance overall – obtaining better results than MWA in 7 instances, while the MWA method obtains better results than HYB in 10 instances (and in one instance the two methods obtain the same average solution value). The two methods require approximately the same average CPU-times for this problem set.

Table 5 compares the best solutions for problems from the MED test set, comparing MWA to ILS, where MWA executes 10 runs and ILS executes 30 runs. Here we present only the results for ILS and MWA because, unfortunately, the best solutions returned by MDAS and HYB are not reported for these particular problems. Table 5 shows that the best solutions found by MWA are better than the best solutions found by ILS in 15 out of 18 instances, and the two methods found the same optimal solutions for 3 instances.

From Table 4, we observe that: (1) For all instance q-10 (q∈ {500, 1000, 1500, 2000, 2500, 3000}) the MWA performs very well by reaching the best average solution cost in all executions while requiring less CPU-time than the state-of-the-art methods from the literature. (2) For instances q-1000 (q∈ {500, 1000, 1500, 2000, 2500, 3000}), the MWA provides high-quality solutions and outperforms the two state-of-the-art methods ILS and MDAS. We observe also that the CPU-times of MWA in these instances are slightly worse than the other methods. Finally, Table 5 shows that the MWA was able, out of 10 runs, to find the optimal solutions for 12 out of 18 large-scale instances (with up to 3000 customers and 3000 facilities).

Results on M* and ORLIB test problems

The results of MWA versus the literature results on M* and ORLIB benchmarks are presented in Table 6 and Table 7 respectively. In Table 6 (for the M* test problems) we compare MWA to TS, LAG, US and HYB and for each instance show the solution costs and CPU-times reported by each method. In Table 7 (for the ORLIB test problems) we present the name of the instance and the optimal solution in columns one and two, then present the CPU-times of TS, LAG and HYB in columns three, four, and five respectively. The CPU-times of MWA are presented in columns six and seven. These tables show that the M* and ORLIB problems were easy to solve.

Table 6 shows that TS, US, and HYB, along with our MWA approach, were all able to reach every best-known solution for the M* problems in competitive CPU-times compared to the other state-of-the-art methods. The ORLIB problems turn out to be even easier to solve, since all state-of-the-art methods, without exception, are able to obtain optimal solutions to all problems. For this reason, Table 7 shows only the CPU-times. It may be noted that our MWA method requires significantly less average computation time for these problems than the other methods, demonstrating its efficiency in solving these problems as well as in the case of the more complex problems, where in addition it obtains better solution results.

Conclusions

Our simple version of the multi-wave algorithm for the uncapacitated facility location problem proves highly effective, as determined by computational tests comparing our method to nine state-of-the-art methods from the literature. The efficacy of our MWA algorithm derives from strategies for exploiting the conditional effects of sequential decisions, embodied in two alternating phases, a vertical phase that successively augments or improves a current solution, and a horizontal phase that modifies the current solution by replacing moves stored in an active move record (AMR). The alternation of vertical and horizontal phases continues until reaching a boundary solution, after which parameters are updated based on performance history and the method repeats. Our findings motivate an investigation of more advanced strategies of the general multi-wave algorithm that additionally exploit other conditional effects, including particularly the effect of persistent attractiveness.

Our findings suggest that the ideas we have incorporated in our multi-wave implementation for UFLP may likewise prove valuable in application to other combinatorial optimization problems. It has been beyond the scope of this study to investigate the full set of strategies available to the multi-wave algorithm, and consequently, our findings also motivate an exploration of other multi-wave components.

References

[1] Ahn S, Cooper C, Cornuéjols G, Frieze A (1988). Probabilistic analysis of a relaxation for the k-median problem. Mathematics of Operations Research, 13(1): 1–31

[2] Akbaripour H, Masehian E, Roostaei A (2017). Landscape analysis and scatter search metaheuristic for solving the uncapacitated single allocation hub location problem. International Journal of Industrial and Systems Engineering, 26(4): 425–459

[3] Al-Sultan K S, Al-Fawzan M A (1999). A tabu search approach to the uncapacitated facility location problem. Annals of Operations Research, 86: 91–103

[4] Albareda-Sambola M, Fernández E, Saldanha-da-Gama F (2017). Heuristic solutions to the facility location problem with general bernoulli demands. INFORMS Journal on Computing, 29(4): 737–753

[5] Amin S H, Zhang G (2013). A multi-objective facility location model for closed-loop supply chain network under uncertain demand and return. Applied Mathematical Modelling, 37(6): 4165–4176

[6] An H C, Svensson O (2017). Recent developments in approximation algorithms for facility location and clustering problems. In: Fukunaga T, Kawarabayashi K, eds. Combinatorial Optimization and Graph Algorithms, 1–19. Singapore: Springer

[7] Ardjmand E, Park N, Weckman G, Amin-Naseri M R (2014). The discrete unconscious search and its application to uncapacitated facility location problem. Computers & Industrial Engineering, 73: 32–40

[8] Arostegui M A Jr, Kadipasaoglu S N, Khumawala B M (2006). An empirical comparison of tabu search, simulated annealing, and genetic algorithms for facilities location problems. International Journal of Production Economics, 103(2): 742–754

[9] Atta S, Mahapatra P R S, Mukhopadhyay A (2018). Deterministic and randomized heuristic algorithms for uncapacitated facility location problem. Information and Decision Sciences, 205–216. Singapore: Springer

[10] Barahona F, Chudak F (1999). Near-optimal Solutions to Large Scale Facility Location Problems. Technical Report RC21606, IBM, USA

[11] Basu S, Sharma M, Ghosh P S (2015). Metaheuristic applications on discrete facility location problems: A survey. OPSEARCH, 52(3): 530–561

[12] Beasley J E (1990). OR-Library: Distributing test problems by electronic mail. Journal of the Operational Research Society, 41(11): 1069–1072

[13] Beasley J E (1993). Lagrangean heuristics for location problems. European Journal of Operational Research, 65(3): 383–399

[14] Beltran C, Tadonki C, Vial J Ph (2006). Solving the p-median problem with a semi-Lagrangian relaxation. Computational Optimization and Applications, 35(2): 239–260

[15] Beltran-Royo C, Vial J P, Alonso-Ayuso A (2012). Semi-lagrangian relaxation applied to the uncapacitated facility location problem. Computational Optimization and Applications, 51(1): 387–409

[16] Blum C, Roli A (2003). Metaheuristics in combinatorial optimization: Overview and conceptual comparison. ACM Computing Surveys, 35(3): 268–308

[17] Cerrone C, Cerulli R, Golden B (2017). Carousel greedy: A generalized greedy algorithm with applications in optimization and statistics. Computers & Operations Research, 85: 97–112

[18] Chalupa D, Nielsen P (2018). Instance scale, numerical properties and design of metaheuristics: A study for the facility location problem. arXiv preprint arXiv:1801.03419

[19] Chen L, Olhager J, Tang O (2014). Manufacturing facility location and sustainability: A literature review and research agenda. International Journal of Production Economics, 149: 154–163

[20] Cheung B S, Langevin A, Villeneuve B (2001). High performing evolutionary techniques for solving complex location problems in industrial system design. Journal of Intelligent Manufacturing, 12(5): 455–466

[21] Cura T (2010). A parallel local search approach to solving the uncapacitated warehouse location problem. Computers & Industrial Engineering, 59(4): 1000–1009

[22] De Armas J, Juan A A, Marques J M, Pedroso J P (2017). Solving the deterministic and stochastic uncapacitated facility location problem: From a heuristic to a simheuristic. Journal of the Operational Research Society, 68(10): 1161–1176

[23] De Corte A, Sörensen K (2015). An iterated local search algorithm for water distribution network design optimization. Networks, 67(3): 187–198

[24] Erlenkotter D (1978). A dual-based procedure for uncapacitated facility location: General solution procedures and computational experience. Operations Research, 26(6): 992–1009

[25] Galli L, Letchford A N, Miller S J (2018). New valid inequalities and facets for the simple plant location problem. European Journal of Operational Research, 269(3): 824–833

[26] Galvão R D, Raggi L A (1989). A method for solving to optimality uncapacitated location problems. Annals of Operations Research, 18(1): 225–244

[27] Ghosh D (2003). Neighborhood search heuristics for the uncapacitated facility location problem. European Journal of Operational Research, 150(1): 150–162

[28] Glover F (1989). Tabu search—Part I. ORSA Journal on Computing, 1(3): 190–206

[29] Glover F (1995). Tabu thresholding: Improved search by nonmonotonic trajectories. ORSA Journal on Computing, 7(4): 426–442

[30] Glover F (2000). Multi-start and strategic oscillation methods—principles to exploit adaptive memory. In: Laguna M, Gonzales-Velarde J L, eds. Computing Tools for Modeling, Optimization and Simulation: Interfaces in Computer Science and Operations Research, 1–24. Berlin: Kluwer Academic Publishers

[31] Glover F (2016). The multi-wave algorithm for metaheuristic optimization. Journal of Heuristics, 22(3): 331–358

[32] Glover F, Laguna M (1997). Tabu Search. Boston: Kluwer Academic Publishers

[33] Goldengorin B, Ghosh D, Sierksma G (2003). Branch and peg algorithms for the simple plant location problem. Computers & Operations Research, 30(7): 967–981

[34] Grasas A, Juan A A, Faulin J, de Armas J, Ramalhinho H (2017). Biased randomization of heuristics using skewed probability distributions: A survey and some applications. Computers & Industrial Engineering, 110: 216–228

[35] Greistorfer P, Rego C (2006). A simple filter-and-fan approach to the facility location problem. Computers & Operations Research, 33(9): 2590–2601

[36] Hale T S, Moberg C R (2003). Location science research: A review. Annals of Operations Research, 123(1/4): 21–35

[37] Han L, Xu D, Du D, Zhang D (2018). A local search approximation algorithm for the uniform capacitated k-facility location problem. Journal of Combinatorial Optimization, 35(2): 409–423

[38] Homberger J, Gehring H (2008). A two-level parallel genetic algorithm for the uncapacitated warehouse location problem. In: Proceedings of Hawaii international conference on system sciences

[39] Jiang Y, Xu D, Du D, Wu C, Zhang D (2018). An approximation algorithm for soft capacitated k-facility location problem. Journal of Combinatorial Optimization, 35(2): 493–511

[40] Jörnsten K, Klose A (2016). An improved Lagrangian relaxation and dual ascent approach to facility location problems. Computational Management Science, 13(3): 317–348

[41] Khumawala B M (1972). An efficient branch and bound algorithm for the warehouse location problem. Management Science, 18(12): 718–731

[42] Klose A, Drexl A (2005). Facility location models for distribution system design. European Journal of Operational Research, 162(1): 4–29

[43] Körkel M (1989). On the exact solution of large-scale simple plant location problems. European Journal of Operational Research, 39(2): 157–173

[44] Kratica J, Tošic D, Filipović V, Ljubić I (2001). Solving the simple plant location problem by genetic algorithm. Operations Research, 35(1): 127–142

[45] Melo M T, Nickel S, Saldanha-Da-Gama F (2009). Facility location and supply chain management–A review. European Journal of Operational Research, 196(2): 401–412

[46] Michel L, Van Hentenryck P (2004). A simple tabu search for warehouse location. European Journal of Operational Research, 157(3): 576–591

[47] Ortiz-Astorquiza C, Contreras I, Laporte G (2017a). Formulations and approximation algorithms for multilevel uncapacitated facility location. INFORMS Journal on Computing, 29(4): 767–779

[48] Ortiz-Astorquiza C, Contrerasn I, Laporte G (2017b). Multi-level facility location problems. European Journal of Operational Research, 267(3): 791–805

[49] Pearce R H, Forbes M (2018). Disaggregated Benders decomposition and branch-and-cut for solving the budget-constrained dynamic uncapacitated facility location and network design problem. European Journal of Operational Research, 270(1): 78–88

[50] Posta M, Ferland J A, Michelon P (2014). An exact cooperative method for the uncapacitated facility location problem. Mathematical Programming Computation, 6(3): 199–231

[51] Resende M G, Werneck R F (2004). A hybrid heuristic for the p-median problem. Journal of Heuristics, 10(1): 59–88

[52] Resende M G C, Ribeiro C C (2003). Greedy randomized adaptive search procedures. In: Glover F, Kochenberger G, eds. Handbook of Metaheuristics, 219–249. Berlin: Kluwer Academic Publishers

[53] Resende M G C, Werneck R F (2006). A hybrid multistart heuristic for the uncapacitated facility location problem. European Journal of Operational Research, 174(1): 54–68

[54] Resende M G C, Werneck R F (2007). A fast swap-based local search procedure for location problems. Annals of Operations Research, 150(1): 205–230

[55] ReVelle C S, Eiselt H A, ReVelle C S (2005). Location analysis: A synthesis and survey. European Journal of Operational Research, 165(1): 1–19

[56] Sahman M A, Altun A A, Dündar A O (2017). The binary differential search algorithm approach for solving uncapacitated facility location problems. Journal of Computational and Theoretical Nanoscience, 14(1): 670–684

[57] Sun M (2005). A tabu search heuristic procedure for the uncapacitated facility location problem. In: Rego C, Alidaee B, eds. Metaheuristic Optimization via Memory and Evolution: Tabu Search and Scatter Search, 191–211. Boston: Kluwer Academic Publishers

[58] Sun M (2006). Solving the uncapacitated facility location problem using tabu search. Computers & Operations Research, 33(9): 2563–2589

[59] Todosijević R, Urošević D, Mladenović N, Hanafi S (2017). A general variable neighborhood search for solving the uncapacitated r-allocation p-hub median problem. Optimization Letters, 11(6): 1109–1121

[60] Tohyama H, Ida K, Matsueda J (2011). A genetic algorithm for the uncapacitated facility location problem. Electronics and Communications in Japan, 94(5): 47–54

[61] Tseng L Y, Wu C S (2009a). Multiple trajectory search for uncapacitated facility location problems. In: Proceedings of the Second International Joint Conference on Computational Sciences and Optimization, Sanya, China. 2: 965–968

[62] Tseng L Y, Wu C S (2009b). The multistart drop-add-swap heuristic for the uncapacitated facility location problem. In: Proceedings of the 6th International Conference on Informatics in Control, Automation and Robotics, Intelligent Control Systems and Optimization, Milan, Italy. 21–28

[63] Tsuya K, Takaya M, Yamamura A (2017). Application of the firefly algorithm to the uncapacitated facility location problem. Journal of Intelligent & Fuzzy Systems, 32(4): 3201–3208

[64] Wang D, Wu C H, Ip A, Wang D, Yan Y (2008). Parallel multipopulation particle swarm optimization algorithm for the uncapacitated facility location problem using OpenMP. In: Proceedings of 2008 IEEE Congress on Evolutionary Computation, IEEE World Congress on Computational Intelligence. 1214–1218

[65] Wang Y, Lu Z, Glover F, Hao J K (2013). Probabilistic GRASP-tabu search algorithms for the UBQP problem. Computers & Operations Research, 40(12): 3100–3107

[66] Xu J, Chiu S Y, Glover F (1997). Probabilistic tabu search for telecommunications network design. Combinatorial Optimization: Theory and Practice, 1(1): 69–94

[67] Yigit V, Aydin M E, Turkbey O (2006). Solving large-scale uncapacitated facility location problems with evolutionary simulated annealing. International Journal of Production Research, 44(22): 4773–4791

RIGHTS & PERMISSIONS

The Author(s) 2018. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)