# Parameter Tuning for Local-Search-Based Matheuristic Methods

Authors: Guillermo Cabrera-Guerrero; Carolina Lagos; Carolina Castañeda; Franklin Johnson; Fernando Paredes; Enrique Cabrera
Journal: Complexity
 (2017)

---

## Abstract

Algorithms that aim to solve optimisation problems by combining heuristics and mathematical programming have attracted researchers’ attention. These methods, also known asmatheuristics, have been shown to perform especially well for large, complex optimisation problems that include both integer and continuous decision variables. One common strategy used by matheuristic methods to solve such optimisation problems is to divide the main optimisation problem into several subproblems. While heuristics are used to seek for promising subproblems, exact methods are used to solve them to optimality. In general, we say that both mixed integer (non)linear programming problems and combinatorial optimisation problems can be addressed using this strategy. Beside the number of parameters researchers need to adjust when using heuristic methods, additional parameters arise when using matheuristic methods. In this paper we focus on one particular parameter, which determines the size of the subproblem. We show how matheuristic performance varies as this parameter is modified. We considered a well-known NP-hard combinatorial optimisation problem, namely, the capacitated facility location problem for our experiments. Based on the obtained results, we discuss the effects of adjusting the size of subproblems that are generated when using matheuristics methods such as the one considered in this paper.


---

## Body

## 1. Introduction

Solving mixed integer programming (MIP) problems as well as combinatorial optimisation problems is, in general, a very difficult task. Although efficient exact methods have been developed to solve these problems to optimality, as the problem size increases exact methods fail to solve it within an acceptable computational time. As a consequence, nonexact methods such as heuristic and metaheuristic algorithms have been developed to find good quality solutions. In addition, hybrid strategies combining different nonexact algorithms are also promising ways to tackle complex optimisation problems. Unfortunately, algorithms that do not consider exact methods cannot give us any guarantee of optimality and, thus, we do not know how good (or bad) solutions found by these methods are.One hybrid strategy that combines nonexact methods are memetic algorithms, which are population-based metaheuristics that use an evolutionary framework integrated with local search algorithms [1]. Memetic algorithms provide lifetime learning process to refine individuals in order to improve the obtained solutions every iteration or generation; their applications have been grown significantly over the years in several NP-hard optimisation problems [2]. These algorithms are part of the paradigm of memetic computation, where the concept of meme is used to automate the knowledge transfer and reuse across problems [3]. A large number of memetic algorithms can be found in the literature. Depending on its implementation, memetic algorithms might (or might not) give guarantee of local optimality: roughly speaking, if heuristic local search algorithms are considered, then no optimality guarantee is given; if exact methods are considered as the local optimisers, then local optimality could be ensured.To overcome the situation described above, hybrid methods that combine heuristics and exact methods to solve optimisation problems have been proposed. These methods, also known as matheuristics [4], have been shown to perform better than both heuristic and mathematical programming methods when they are applied separately. The idea of combining the power of mathematical programming with flexibility of heuristics has gained attention within researchers’ community. We can found matheuristics attempting to solve problems arising in the field of logistics [5–9], health care systems [10–13], and pure mathematics [14, 15], among others. Matheuristics have been demonstrated to be very effective in solving complex optimisation problems. Some interesting surveys on matheuristics are [16, 17]. Although there is some overlapping between memetic algorithms and matheuristic ones, in this paper we have chosen to label the studied strategy as matheuristic, as we think matheuristics definition better fits the framework we are interested in.Because of their complexity, MIP problems as well as combinatorial optimisation problems are often tackled using matheuristic methods. One common strategy to solve this class of optimisation problems is to divide the main optimisation problem into several subproblems. While heuristics are used to seek for promising subproblems, exact methods are used to solve them to optimality. One advantage of this approach is that it does not depend on the (non)linearity of the resulting subproblem. Instead, it has been pointed out that it is desirable that the resulting subproblem would be convex [11]. Having a convex subproblem would allow us to solve it to optimality, and thus comparing solutions obtained at each subproblem becomes more senseful. This strategy has been successfully applied to problems arising in fields as diverse as logistics and radiation therapy.In this paper we aim to study the impact of parameter tuning on the performance of matheuristic methods as the one described above. To this end, the well-known capacitated facility location problem is used as an application of hard combinatorial optimisation problem. To the best of our knowledge, no paper has focused on parameter tuning for matheuristic methods.This paper is organised as follows: Section2 shows the general matheuristic framework we consider in this paper. Details on the algorithms that are used in this study are also shown in this section. In Section 3 the capacitated facility location problem is introduced and its mathematical model is described in Section 3.1. The experiments performed in this study are presented and the obtained results are discussed in Section 3.2. Finally, in Section 4 some conclusions are presented and the future work is outlined.


## 2. Matheuristic Methods

This section is twofold. We start by describing a general matheuristic framework that is used to solve both MIP problems and combinatorial optimisation problems and how it is different from other commonly used approaches such as memetic algorithms and other evolutionary approaches. After that, we present the local-search-based algorithms we consider in this work to perform our experiments. We finish this section by introducing the parameter we will be focused on in this study.


### 2.1. General Framework

Equations (1a) to (1f) show the general form of MIP problems. Hereafter we will refer to this problem as the MIP(x,y) problem or themain problem.(1a)MIPx,y:minx,yfx,y(1b)s.t.gjx,y≤bjforj=1,…,m,(1c)gj-x≤bj-forj=1,…,m-,(1d)gj¯¯y≤bj̿forj=1,…,m̿,(1e)y∈0,1o,(1f)x∈Rn,where f(x,y) is an objective function, m is the number of inequality constraints on x and y, m- is the number of inequality constraints on x, m̿ is the number of inequality constraints on y, o is the number of binary (≥0) decision variables y, and n is the number of continuous decision variables. Combinatorial optimisation problems, as the one we consider in this paper, can be easily obtained by either removing the continuous decision variable from the model or making it integer (i.e., x∈Z≥o).Although there exist a number of exact algorithms that can find an optimal solution for the MIP(x,y), as the size of the problem increases, exact methods fail or take too long. Because of this, heuristic methods are used to obtain good quality solutions of the problem within an acceptable time. Heuristic methods cannot guarantee optimality though.During the last two decades, the idea of combining heuristic methods and mathematical programming has received much more attention. Exploiting the advantages of each method appears to be a senseful strategy to overcome their inherent drawbacks. Several strategies have been proposed to combine heuristics and exact methods to solve optimisation problems such as the MIP(x,y) problem. For instance, Chen and Ting [18] and Lagos et al. [8] combine the well-known ant colony optimisation (ACO) algorithm and Lagrangian relaxation method to solve the single source capacitated facility location problem and a distribution network design problem, respectively. In these articles, Lagrangian multipliers are updated using ACO algorithm. Another strategy to combine heuristics and exact methods is to let heuristics seek for subproblems of MIP(x,y) which, in turn, are solved to optimality by some exact method. One alternative to obtain subproblems of MIP(x,y) is to add a set of additional constraints on a subset of binary decision variables. These constraints are of the form yi=0, with i∈I and I being the set of index that are restricted in subproblem MIPI(x,y) (see (1a) to (1f)). The portion of binary decision variables yi that are set to 0 is denoted by α (i.e., #I=α×o), with o being the number of binary variables yi. Then, the obtained subproblem, which we call MIPI(x,y), is(2a)MIPIx,y:minx,yfx,y(2b)s.t.gjx,y≤bjforj=1,…,m,(2c)gj-x≤bj-forj=1,…,m¯,(2d)gj̿y≤bj̿forj=1,…,m̿,(2e)yi=0i∈I,(2f)y∈0,1o,(2g)x∈Rn.In this paper we assume that a constraint on theith resource of vector y of the form yi=0, that is, i∈I, means that such a resource is not available for subproblem MIPI(x,y). We can note that, as the number of constrained binary decision variables yi associated with I increases, that is, as (1-α) increases, the subproblem MIPI(x,y) becomes smaller. Similarly, as the number of constrained binary decision variables yi associated with I gets smaller, that is, (1-α) decreases, the subproblem MIPI(x,y) gets larger, as there are more available resources. It is also easy to note that, as (1-α) increases, we can obtain optimal solutions of the corresponding subproblem MIPI(x,y) relatively faster. However, quality of the obtained solutions is usually impaired as the solution space is restricted by the additional constraints associated with the set I. Similarly, as (1-α) gets smaller, obtaining optimal solutions of the associated subproblems might take longer but the quality of the obtained solutions is, in general, greatly improved. Finally, when α=0 we have that subproblem MIPI(x,y) is identical to the main problem, MIP(x,y). We assume that optimal solution of subproblem MIPI(x,y), (x^,y^), is also feasible for the MIP(x,y) problem. Moreover, we can note that there must be a minimal set I, for which an optimal solution (x^,y^) of the associated subproblem MIPI(x,y) is also an optimal solution (x∗,y∗) of the main problem MIP(x,y).In some cases, the value ofα might be predefined by the problem that is being solved. For instance, the problem of finding the best beam angle configuration for radiation delivery in cancer treatment (beam angle optimisation problem) usually sets the number of beams to be used in a beam angle configuration (see Cabrera-Guerrero et al. [11], Li et al. [19], and Li et al. [20]). This definition is made by the treatment planner and it does not take into account the algorithm performance but clinical aspects. Unlike this kind of problems, there are many other problems where the value of α is not predefined and, then, setting it to an efficient value is important for the algorithm performance. Figure 1 shows the interaction between the heuristic and the exact method.Figure 1
Interaction between a heuristic method and an exact method.One distinctive feature of matheuristics is that there exists an interaction between the heuristic method and the exact method. In the method depicted in Figure1 we have that, on one hand, the heuristic method influences the solver by passing onto it a set of constraints, I, that defines the subproblem to be solved by the exact method. On the other hand, we have that the solution obtained by the solver is returned to the heuristic method such that it can be used to influence the selection of the elements in the next set of constraints; that is, the solution of a subproblem might be used by the heuristic to obtain some useful information to generate next subproblems.As mentioned above, there is one parameter that is not part of the set of parameters of the heuristic method nor part of the set of parameters of the exact method. This parameter, which we callα, only comes up after these two methods are posed together. Then, in this paper we are interested in how the choice of α can modify the performance of the proposed strategy in terms of the quality of the obtained solutions of the main problem. Since many different matheuristic frameworks might be proposed to solve the MIP(x,y) problem, we restrict this study to local-search-based matheuristics. Thus, in this study, three local-search-based matheuristic algorithms are implemented and their results compared. Additionally, we implement a very simple method which we call “blind algorithm” as a baseline for this study. Next sections explain these algorithms.


### 2.2. Local Search Algorithms

As mentioned in previous sections, the aim of this paper is not to provide a “state-of-the-art” algorithm to solve MIP problems but, instead, to study the effect that changing the size of subproblems has on the quality of the obtained solutions when using local-search-based matheuristic methods as the one described in Section2.1. Thus, three local search algorithms are implemented, namely, steepest descent (SD), next descent (ND), and tabu search (TS). Local search algorithms need a neighbourhood N to be defined by means of a neighbourhood movement. We define the same neighbourhood movement for the all three methods.Since the local search algorithms move on the subproblem space, that is, the local search algorithms look for promising subproblemsMIPI(x,y), or, equivalently, look for promising sets I, then we need to define a neighbourhood within the same search space. Thus, the neighbour set of I, which we denote by N(I) or simply I′, corresponds to all resources yi such that they either are not considered in the optimal solution of MIPI(x,y) or are actually part of I. Since the number of resources yi that meet these criteria might be above or below the number of elements set I′ must contain according to the parameter α that is being considered, we randomly select among the variables that meet these criteria until the neighbour set I′ is completed. This neighbourhood movement ensures that those resources yi that are part of the optimal solution of subproblem MIPI(x,y), that is, yi∉I and yi=1, will continue to be available for sets I′∈N(I). Thus, the set of resources yi that can potentially be part of sets I′∈N(I) is defined as follows:(3)I′∈NI=yi:yi∈Ioryi∉I,y^i=0,where (x^,y^) is the optimal solution of MIPI(x,y). The initial set I local search methods start with is set randomly for the all three algorithms implemented in this paper.


#### 2.2.1. Steepest Descent Algorithm

The steepest descent algorithm starts with an initialI0 for which its associated subproblem, MIPI0(x,y), is labelled as the current subproblem. For this current subproblem, its entire neighbourhood N(I0) is generated, and the neighbour that leads to the best objective function value is selected. If the best neighbour subproblem is better than the current one, then the best neighbour is set as the new current solution. If the best neighbour is not better than the current subproblem, then the algorithm stops and the optimal solution of the current subproblem is returned as a locally optimal solution of the main MIP(x,y) problem. Algorithm 1 shows the pseudocode for the steepest descent algorithm implemented in this paper.Algorithm 1:Steepest descent algorithm.
Input:  α (portion of the resources included in I)
Output:   (x^,y^) (locally optimal solution)
(1) begin
(2)  k=0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (x^,y^)= solve MIPIk(x,y);
(5)  repeat
(6)   localOptimum = true;
(7)   k=k+1;
(8)   foreach  Ik-1′∈N(Ik-1)  do
(9)    (x′,y′)= solve MIPIk-1′(x,y);
(10)    If  MIPx′,y′<MIPx^,y^  then
(11)      (x^,y^)=(x′,y′);
(12)      Ik=Ik-1′;
(13)      localOptimum = false;
(14)  until localOptimum;
(15)  return(x^,y^);Although the steepest descent algorithm can be considered, in general, a deterministic algorithm, in the sense that, given an initial solution, it converges to the same local optima, in our implementation the algorithm does not visit all possible neighbours and, therefore, it becomes a stochastic local search. Since only one local optimum is generated at each run, we repeat the algorithm until the time limit is reached. Same is done for both ND and TS algorithms we introduce next.


#### 2.2.2. Next Descent Algorithm

As mentioned in Section2.2.1, the steepest descent might take too long to converge if the size of the neighbourhood is too large or solving an individual subproblem is too time-consuming. Thus, we include in this paper a local search algorithm called next descent that aims to converge faster than the steepest descent, without a major impact on the solution quality. Algorithm 2 shows the next descent method.Algorithm 2:Next descent algorithm.
Input: α (portion of the resources included in I)
Output:  (x^,y^) (locally optimal solution)
(1) begin
(2)  k=0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (x^,y^)= solve MIPIk(x,y);
(5)   repeat
(6)   localOptimum = true;
(7)   k=k+1;
(8)   foreach  Ik-1′∈N(Ik-1)do
(9)     (x′,y′)= solve MIPIk-1′(x,y);
(10)     if  MIP(x′,y′)<MIP(x^,y^)  then
(11)      (x^,y^)=(x′,y′);
(12)      Ik=Ik-1′;
(13)      localOptimum = false;
(14)      break;
(15)  until localOptimum;
(16)  return (x^,y^);Just as in the steepest descent algorithm, the next descent algorithm starts with an initial solution, which is labelled as the current solution. Then, a random element from the neighbourhood of the current solution is selected and solved, and the objective function value of its optimal solution is compared to the objective function value of the current solution. Like in the SD algorithm, if the neighbour solution is not better than the current solution, another randomly selected setI′ from the neighbourhood N(Ik-1) is generated and compared to the current solution. Unlike the SD algorithm, in the ND algorithm, if the neighbour is better than the current solution, then the neighbour is labelled as the new current solution and no other I′∈N(Ik-1) is visited. The algorithm repeats these steps until the entire neighbourhood has been computed with no neighbour resulting in a better solution than the current one, in which case the algorithm stops.


#### 2.2.3. Tabu Search Algorithm

Unlike the algorithms described in the previous sections, tabu search is a local search technique guided by the use of adaptive or flexible memory structures [5], and thus, it is inherently a stochastic local search algorithm. The variety of the tools and search principles introduced and described in [21] are such that the TS can be considered as the seed of a general framework for modern heuristic search [22]. We include TS as it has been applied to several combinatorial optimisation problems (see, e.g., [5, 23–27]) including, of course, mixed integer programming problems as the one we consider in this study.As in the algorithms introduced above, TS also starts with an initial set of constraintsI0 for which its associated subproblem MIPI0(x,y) is solved. Then, as in the SD algorithm, the “entire” neighbourhood of the current solution is computed. As explained before, since the neighbourhood N has a stochastic component, it is actually not possible to generate the entire neighbourhood of a solution; however, we use the term “entire” to stress the fact that all thegenerated neighbours of N are solved. This is different from the ND algorithm, where not all the generated neighbours are necessarily solved. The number of generated neighbours, ns, is, as in the algorithms above, set equal to 10. Moreover, TS implements a list calledtabu list that aims to avoid cycles during the search. Each time a neighbourhood is generated, the warehouses removed from I(k-1) are marked as tabu. Once we have solved the subproblems MIPIk′(x,y), with Ik′∈N(Ik), its optimal solutions are ranked and the one with the best objective function value is chosen as candidate solution for the next iteration. If the candidate solution was generated using a movement within the tabu list, it should be discarded, and the next best neighbour of the list should be chosen as the new candidate solution. However, there is one exception to this rule: if the candidate solution is within the tabu list but its objective function value is better than the best objective function found so far by the algorithm, then the so calledaspiration criterion is invoked and the candidate solution is set as the new current solution and passed on to the next iteration.Unlike the algorithms introduced above, TS does not require next current solution to be better than the previous one. This means that it is able to avoid local optimal by choosing solutions that are more expensive as they allows the algorithm to visit other (hopefully promising) areas in the search space. However, in case the algorithm cannot make any improvement after a predefined number of iterations, a diversification mechanism is used to get out from low-quality neighbourhoods and “jump” to other neighbourhoods. The diversification mechanism implemented here is a restart method, which set the current solution to a randomly generated solution without losing the best solution found so far. Termination criterion implemented here is the time limit.The tabu search implemented in this paper is as follows.As Algorithm3 shows, tabu search requires the following parameters:(i)
Time limit: total time the algorithm will perform.(ii)
DiversBound: total number of iteration without improvements on the best solution before diversification criterion (restart method) is applied.(iii)
TabuListSize(ts): size oftabuList. Number of iterations for which a specific movement remains banned.Algorithm 3:Tabu search algorithm.
Input: α (portion of the resources included in Ik)
Output:  (x^,y^) (approximately optimal solution)
(1) begin
(2)  k=0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (xc,yc)= solve MIPk(x,y);   // set current sol
(5)  (x^,y^)=(xc,yc); //set best sol
(6)  Y=∅; // tabu list initially empty
(7)  noImprovementCounter=0;
(8)  repeat
(9)    C=∅;
(10)   k=k+1;
(11)   foreachIk-1′∈N(Ik-1′)   do
(12)    (xn,yn)= solve MIPIk-1′(x,y);
(13)    C=C∪xn,yn;
(14)   SortC;
(15)   foreach  {(xn,yn)}∈C  do
(16)    if  isTabu((xn,yn))  then
(17)     if  MIP(xn,yn)<MIP(x^,y^)  then
(18)      Y=updateTabuList(Y,(xc,yc),(xn,yn),ts);
(19)      (xc,yc)=(xn,yn);
(20)      (x^,y^)=(xn,yn);
(21)      Ik=Ik-1n;
(22)      noImprovementCounter=0;
(23)      break;
(24)    else
(25)     Y=updateTabuList(Y,(xc,yc),(xn,yn));
(26)     (xc,yc)=(xn,yn);
(27)     Ik=Ik-1n  if  MIP(xn,yn)<MIP(x^,y^)  then
(28)       (x^,y^)=(xn,yn);
(29)       noImprovementCounter=0;
(30)     else
(31)       noImprovementCounter++;
(32)     break;
(33)   if  noImprovementCounter≥diversBound  then
(34)    Ik=selectBinaryVariablesRandomly(α);
(35)  until  time limit is reached;
(36)  return  (x^,y^);


#### 2.2.4. Blind Algorithm

We finally implement a very simple heuristic method we callblind as it moves randomly at each iteration. Pseudocode of this method is presented in Algorithm 4.Algorithm 4:Matheuristic method using the blind algorithm.
// To find an (approximately) optimal solution for MIP
(
x
,
y
)
Input:  α (portion of the resources included in Ik)
Output:   (x^,y^) (approximately optimal solution)
(1) begin
(2)  k = 0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (x^,y^) = solve (MIPIk(x,y));
(5)  while  !stopCriterion  do
(6)   k=k+1;
(7)   Ik=selectBinaryVariablesRandomly(α);
(8)   (xc,yc) = solve (MIPIk(x,y));
(9)   if  MIP(xc,yc)<MIP(x^,y^)  then
(10)    (x^,y^) = xc,yc;
(11)    Ik=Ik-1
(12)  return (x^,y^);As we can see, no additional intelligence is added to the blind algorithm. It is just a random search that, after a predefined number of iterations (or any other “stop criterion”), returns the best solution it found during its search. Thus, we can consider this algorithm as a baseline of this study.We apply the all four algorithms described in this section to two prostate cases. Details on this case and the obtained results are presented in the next section.


### 2.3. Subproblem Sizing

All the algorithms above perform very different as the value of the input parameterα varies. On the one hand, setting α to a very small value provokes that either the solver fails to solve the subproblem because of lack of memory or it takes too long to find an optimal solution of the subproblem. On the other hand, setting α close to the total number of binary decision variables provokes that the algorithm fails to find a solution for the generated subproblems as there is no feasible solution to it (i.e., the subproblems are too restrictive). Thus, we have to find a value of α such that exact methods can solve the obtained subproblems within an acceptable time. Further, α should allow local search methods to iterate as much as needed. Solving the obtained subproblems within few seconds is critical to matheuristic methods as they usually need several iterations before to converge to a good quality solution. This is especially true for matheuristics that consider local search or population-based heuristic methods. Then, there is a trade-off between the quality of the solution of subproblems and the time that is needed to generate such solutions. Therefore, finding a value of α that gives us a good compromise between these two aspects is critical for the overall performance of the matheuristic methods explained before. It is interesting that the problem of findingefficient values of α might be seen as a multiobjective optimisation problem.In next section we explain the experiments that we perform to study how the choice ofα hits local-search-based matheuristics performance. Based on the results, we draft some guidelines to set value of parameter α at the end of next section.


## 2.1. General Framework

Equations (1a) to (1f) show the general form of MIP problems. Hereafter we will refer to this problem as the MIP(x,y) problem or themain problem.(1a)MIPx,y:minx,yfx,y(1b)s.t.gjx,y≤bjforj=1,…,m,(1c)gj-x≤bj-forj=1,…,m-,(1d)gj¯¯y≤bj̿forj=1,…,m̿,(1e)y∈0,1o,(1f)x∈Rn,where f(x,y) is an objective function, m is the number of inequality constraints on x and y, m- is the number of inequality constraints on x, m̿ is the number of inequality constraints on y, o is the number of binary (≥0) decision variables y, and n is the number of continuous decision variables. Combinatorial optimisation problems, as the one we consider in this paper, can be easily obtained by either removing the continuous decision variable from the model or making it integer (i.e., x∈Z≥o).Although there exist a number of exact algorithms that can find an optimal solution for the MIP(x,y), as the size of the problem increases, exact methods fail or take too long. Because of this, heuristic methods are used to obtain good quality solutions of the problem within an acceptable time. Heuristic methods cannot guarantee optimality though.During the last two decades, the idea of combining heuristic methods and mathematical programming has received much more attention. Exploiting the advantages of each method appears to be a senseful strategy to overcome their inherent drawbacks. Several strategies have been proposed to combine heuristics and exact methods to solve optimisation problems such as the MIP(x,y) problem. For instance, Chen and Ting [18] and Lagos et al. [8] combine the well-known ant colony optimisation (ACO) algorithm and Lagrangian relaxation method to solve the single source capacitated facility location problem and a distribution network design problem, respectively. In these articles, Lagrangian multipliers are updated using ACO algorithm. Another strategy to combine heuristics and exact methods is to let heuristics seek for subproblems of MIP(x,y) which, in turn, are solved to optimality by some exact method. One alternative to obtain subproblems of MIP(x,y) is to add a set of additional constraints on a subset of binary decision variables. These constraints are of the form yi=0, with i∈I and I being the set of index that are restricted in subproblem MIPI(x,y) (see (1a) to (1f)). The portion of binary decision variables yi that are set to 0 is denoted by α (i.e., #I=α×o), with o being the number of binary variables yi. Then, the obtained subproblem, which we call MIPI(x,y), is(2a)MIPIx,y:minx,yfx,y(2b)s.t.gjx,y≤bjforj=1,…,m,(2c)gj-x≤bj-forj=1,…,m¯,(2d)gj̿y≤bj̿forj=1,…,m̿,(2e)yi=0i∈I,(2f)y∈0,1o,(2g)x∈Rn.In this paper we assume that a constraint on theith resource of vector y of the form yi=0, that is, i∈I, means that such a resource is not available for subproblem MIPI(x,y). We can note that, as the number of constrained binary decision variables yi associated with I increases, that is, as (1-α) increases, the subproblem MIPI(x,y) becomes smaller. Similarly, as the number of constrained binary decision variables yi associated with I gets smaller, that is, (1-α) decreases, the subproblem MIPI(x,y) gets larger, as there are more available resources. It is also easy to note that, as (1-α) increases, we can obtain optimal solutions of the corresponding subproblem MIPI(x,y) relatively faster. However, quality of the obtained solutions is usually impaired as the solution space is restricted by the additional constraints associated with the set I. Similarly, as (1-α) gets smaller, obtaining optimal solutions of the associated subproblems might take longer but the quality of the obtained solutions is, in general, greatly improved. Finally, when α=0 we have that subproblem MIPI(x,y) is identical to the main problem, MIP(x,y). We assume that optimal solution of subproblem MIPI(x,y), (x^,y^), is also feasible for the MIP(x,y) problem. Moreover, we can note that there must be a minimal set I, for which an optimal solution (x^,y^) of the associated subproblem MIPI(x,y) is also an optimal solution (x∗,y∗) of the main problem MIP(x,y).In some cases, the value ofα might be predefined by the problem that is being solved. For instance, the problem of finding the best beam angle configuration for radiation delivery in cancer treatment (beam angle optimisation problem) usually sets the number of beams to be used in a beam angle configuration (see Cabrera-Guerrero et al. [11], Li et al. [19], and Li et al. [20]). This definition is made by the treatment planner and it does not take into account the algorithm performance but clinical aspects. Unlike this kind of problems, there are many other problems where the value of α is not predefined and, then, setting it to an efficient value is important for the algorithm performance. Figure 1 shows the interaction between the heuristic and the exact method.Figure 1
Interaction between a heuristic method and an exact method.One distinctive feature of matheuristics is that there exists an interaction between the heuristic method and the exact method. In the method depicted in Figure1 we have that, on one hand, the heuristic method influences the solver by passing onto it a set of constraints, I, that defines the subproblem to be solved by the exact method. On the other hand, we have that the solution obtained by the solver is returned to the heuristic method such that it can be used to influence the selection of the elements in the next set of constraints; that is, the solution of a subproblem might be used by the heuristic to obtain some useful information to generate next subproblems.As mentioned above, there is one parameter that is not part of the set of parameters of the heuristic method nor part of the set of parameters of the exact method. This parameter, which we callα, only comes up after these two methods are posed together. Then, in this paper we are interested in how the choice of α can modify the performance of the proposed strategy in terms of the quality of the obtained solutions of the main problem. Since many different matheuristic frameworks might be proposed to solve the MIP(x,y) problem, we restrict this study to local-search-based matheuristics. Thus, in this study, three local-search-based matheuristic algorithms are implemented and their results compared. Additionally, we implement a very simple method which we call “blind algorithm” as a baseline for this study. Next sections explain these algorithms.


## 2.2. Local Search Algorithms

As mentioned in previous sections, the aim of this paper is not to provide a “state-of-the-art” algorithm to solve MIP problems but, instead, to study the effect that changing the size of subproblems has on the quality of the obtained solutions when using local-search-based matheuristic methods as the one described in Section2.1. Thus, three local search algorithms are implemented, namely, steepest descent (SD), next descent (ND), and tabu search (TS). Local search algorithms need a neighbourhood N to be defined by means of a neighbourhood movement. We define the same neighbourhood movement for the all three methods.Since the local search algorithms move on the subproblem space, that is, the local search algorithms look for promising subproblemsMIPI(x,y), or, equivalently, look for promising sets I, then we need to define a neighbourhood within the same search space. Thus, the neighbour set of I, which we denote by N(I) or simply I′, corresponds to all resources yi such that they either are not considered in the optimal solution of MIPI(x,y) or are actually part of I. Since the number of resources yi that meet these criteria might be above or below the number of elements set I′ must contain according to the parameter α that is being considered, we randomly select among the variables that meet these criteria until the neighbour set I′ is completed. This neighbourhood movement ensures that those resources yi that are part of the optimal solution of subproblem MIPI(x,y), that is, yi∉I and yi=1, will continue to be available for sets I′∈N(I). Thus, the set of resources yi that can potentially be part of sets I′∈N(I) is defined as follows:(3)I′∈NI=yi:yi∈Ioryi∉I,y^i=0,where (x^,y^) is the optimal solution of MIPI(x,y). The initial set I local search methods start with is set randomly for the all three algorithms implemented in this paper.


### 2.2.1. Steepest Descent Algorithm

The steepest descent algorithm starts with an initialI0 for which its associated subproblem, MIPI0(x,y), is labelled as the current subproblem. For this current subproblem, its entire neighbourhood N(I0) is generated, and the neighbour that leads to the best objective function value is selected. If the best neighbour subproblem is better than the current one, then the best neighbour is set as the new current solution. If the best neighbour is not better than the current subproblem, then the algorithm stops and the optimal solution of the current subproblem is returned as a locally optimal solution of the main MIP(x,y) problem. Algorithm 1 shows the pseudocode for the steepest descent algorithm implemented in this paper.Algorithm 1:Steepest descent algorithm.
Input:  α (portion of the resources included in I)
Output:   (x^,y^) (locally optimal solution)
(1) begin
(2)  k=0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (x^,y^)= solve MIPIk(x,y);
(5)  repeat
(6)   localOptimum = true;
(7)   k=k+1;
(8)   foreach  Ik-1′∈N(Ik-1)  do
(9)    (x′,y′)= solve MIPIk-1′(x,y);
(10)    If  MIPx′,y′<MIPx^,y^  then
(11)      (x^,y^)=(x′,y′);
(12)      Ik=Ik-1′;
(13)      localOptimum = false;
(14)  until localOptimum;
(15)  return(x^,y^);Although the steepest descent algorithm can be considered, in general, a deterministic algorithm, in the sense that, given an initial solution, it converges to the same local optima, in our implementation the algorithm does not visit all possible neighbours and, therefore, it becomes a stochastic local search. Since only one local optimum is generated at each run, we repeat the algorithm until the time limit is reached. Same is done for both ND and TS algorithms we introduce next.


### 2.2.2. Next Descent Algorithm

As mentioned in Section2.2.1, the steepest descent might take too long to converge if the size of the neighbourhood is too large or solving an individual subproblem is too time-consuming. Thus, we include in this paper a local search algorithm called next descent that aims to converge faster than the steepest descent, without a major impact on the solution quality. Algorithm 2 shows the next descent method.Algorithm 2:Next descent algorithm.
Input: α (portion of the resources included in I)
Output:  (x^,y^) (locally optimal solution)
(1) begin
(2)  k=0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (x^,y^)= solve MIPIk(x,y);
(5)   repeat
(6)   localOptimum = true;
(7)   k=k+1;
(8)   foreach  Ik-1′∈N(Ik-1)do
(9)     (x′,y′)= solve MIPIk-1′(x,y);
(10)     if  MIP(x′,y′)<MIP(x^,y^)  then
(11)      (x^,y^)=(x′,y′);
(12)      Ik=Ik-1′;
(13)      localOptimum = false;
(14)      break;
(15)  until localOptimum;
(16)  return (x^,y^);Just as in the steepest descent algorithm, the next descent algorithm starts with an initial solution, which is labelled as the current solution. Then, a random element from the neighbourhood of the current solution is selected and solved, and the objective function value of its optimal solution is compared to the objective function value of the current solution. Like in the SD algorithm, if the neighbour solution is not better than the current solution, another randomly selected setI′ from the neighbourhood N(Ik-1) is generated and compared to the current solution. Unlike the SD algorithm, in the ND algorithm, if the neighbour is better than the current solution, then the neighbour is labelled as the new current solution and no other I′∈N(Ik-1) is visited. The algorithm repeats these steps until the entire neighbourhood has been computed with no neighbour resulting in a better solution than the current one, in which case the algorithm stops.


### 2.2.3. Tabu Search Algorithm

Unlike the algorithms described in the previous sections, tabu search is a local search technique guided by the use of adaptive or flexible memory structures [5], and thus, it is inherently a stochastic local search algorithm. The variety of the tools and search principles introduced and described in [21] are such that the TS can be considered as the seed of a general framework for modern heuristic search [22]. We include TS as it has been applied to several combinatorial optimisation problems (see, e.g., [5, 23–27]) including, of course, mixed integer programming problems as the one we consider in this study.As in the algorithms introduced above, TS also starts with an initial set of constraintsI0 for which its associated subproblem MIPI0(x,y) is solved. Then, as in the SD algorithm, the “entire” neighbourhood of the current solution is computed. As explained before, since the neighbourhood N has a stochastic component, it is actually not possible to generate the entire neighbourhood of a solution; however, we use the term “entire” to stress the fact that all thegenerated neighbours of N are solved. This is different from the ND algorithm, where not all the generated neighbours are necessarily solved. The number of generated neighbours, ns, is, as in the algorithms above, set equal to 10. Moreover, TS implements a list calledtabu list that aims to avoid cycles during the search. Each time a neighbourhood is generated, the warehouses removed from I(k-1) are marked as tabu. Once we have solved the subproblems MIPIk′(x,y), with Ik′∈N(Ik), its optimal solutions are ranked and the one with the best objective function value is chosen as candidate solution for the next iteration. If the candidate solution was generated using a movement within the tabu list, it should be discarded, and the next best neighbour of the list should be chosen as the new candidate solution. However, there is one exception to this rule: if the candidate solution is within the tabu list but its objective function value is better than the best objective function found so far by the algorithm, then the so calledaspiration criterion is invoked and the candidate solution is set as the new current solution and passed on to the next iteration.Unlike the algorithms introduced above, TS does not require next current solution to be better than the previous one. This means that it is able to avoid local optimal by choosing solutions that are more expensive as they allows the algorithm to visit other (hopefully promising) areas in the search space. However, in case the algorithm cannot make any improvement after a predefined number of iterations, a diversification mechanism is used to get out from low-quality neighbourhoods and “jump” to other neighbourhoods. The diversification mechanism implemented here is a restart method, which set the current solution to a randomly generated solution without losing the best solution found so far. Termination criterion implemented here is the time limit.The tabu search implemented in this paper is as follows.As Algorithm3 shows, tabu search requires the following parameters:(i)
Time limit: total time the algorithm will perform.(ii)
DiversBound: total number of iteration without improvements on the best solution before diversification criterion (restart method) is applied.(iii)
TabuListSize(ts): size oftabuList. Number of iterations for which a specific movement remains banned.Algorithm 3:Tabu search algorithm.
Input: α (portion of the resources included in Ik)
Output:  (x^,y^) (approximately optimal solution)
(1) begin
(2)  k=0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (xc,yc)= solve MIPk(x,y);   // set current sol
(5)  (x^,y^)=(xc,yc); //set best sol
(6)  Y=∅; // tabu list initially empty
(7)  noImprovementCounter=0;
(8)  repeat
(9)    C=∅;
(10)   k=k+1;
(11)   foreachIk-1′∈N(Ik-1′)   do
(12)    (xn,yn)= solve MIPIk-1′(x,y);
(13)    C=C∪xn,yn;
(14)   SortC;
(15)   foreach  {(xn,yn)}∈C  do
(16)    if  isTabu((xn,yn))  then
(17)     if  MIP(xn,yn)<MIP(x^,y^)  then
(18)      Y=updateTabuList(Y,(xc,yc),(xn,yn),ts);
(19)      (xc,yc)=(xn,yn);
(20)      (x^,y^)=(xn,yn);
(21)      Ik=Ik-1n;
(22)      noImprovementCounter=0;
(23)      break;
(24)    else
(25)     Y=updateTabuList(Y,(xc,yc),(xn,yn));
(26)     (xc,yc)=(xn,yn);
(27)     Ik=Ik-1n  if  MIP(xn,yn)<MIP(x^,y^)  then
(28)       (x^,y^)=(xn,yn);
(29)       noImprovementCounter=0;
(30)     else
(31)       noImprovementCounter++;
(32)     break;
(33)   if  noImprovementCounter≥diversBound  then
(34)    Ik=selectBinaryVariablesRandomly(α);
(35)  until  time limit is reached;
(36)  return  (x^,y^);


### 2.2.4. Blind Algorithm

We finally implement a very simple heuristic method we callblind as it moves randomly at each iteration. Pseudocode of this method is presented in Algorithm 4.Algorithm 4:Matheuristic method using the blind algorithm.
// To find an (approximately) optimal solution for MIP
(
x
,
y
)
Input:  α (portion of the resources included in Ik)
Output:   (x^,y^) (approximately optimal solution)
(1) begin
(2)  k = 0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (x^,y^) = solve (MIPIk(x,y));
(5)  while  !stopCriterion  do
(6)   k=k+1;
(7)   Ik=selectBinaryVariablesRandomly(α);
(8)   (xc,yc) = solve (MIPIk(x,y));
(9)   if  MIP(xc,yc)<MIP(x^,y^)  then
(10)    (x^,y^) = xc,yc;
(11)    Ik=Ik-1
(12)  return (x^,y^);As we can see, no additional intelligence is added to the blind algorithm. It is just a random search that, after a predefined number of iterations (or any other “stop criterion”), returns the best solution it found during its search. Thus, we can consider this algorithm as a baseline of this study.We apply the all four algorithms described in this section to two prostate cases. Details on this case and the obtained results are presented in the next section.


## 2.2.1. Steepest Descent Algorithm

The steepest descent algorithm starts with an initialI0 for which its associated subproblem, MIPI0(x,y), is labelled as the current subproblem. For this current subproblem, its entire neighbourhood N(I0) is generated, and the neighbour that leads to the best objective function value is selected. If the best neighbour subproblem is better than the current one, then the best neighbour is set as the new current solution. If the best neighbour is not better than the current subproblem, then the algorithm stops and the optimal solution of the current subproblem is returned as a locally optimal solution of the main MIP(x,y) problem. Algorithm 1 shows the pseudocode for the steepest descent algorithm implemented in this paper.Algorithm 1:Steepest descent algorithm.
Input:  α (portion of the resources included in I)
Output:   (x^,y^) (locally optimal solution)
(1) begin
(2)  k=0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (x^,y^)= solve MIPIk(x,y);
(5)  repeat
(6)   localOptimum = true;
(7)   k=k+1;
(8)   foreach  Ik-1′∈N(Ik-1)  do
(9)    (x′,y′)= solve MIPIk-1′(x,y);
(10)    If  MIPx′,y′<MIPx^,y^  then
(11)      (x^,y^)=(x′,y′);
(12)      Ik=Ik-1′;
(13)      localOptimum = false;
(14)  until localOptimum;
(15)  return(x^,y^);Although the steepest descent algorithm can be considered, in general, a deterministic algorithm, in the sense that, given an initial solution, it converges to the same local optima, in our implementation the algorithm does not visit all possible neighbours and, therefore, it becomes a stochastic local search. Since only one local optimum is generated at each run, we repeat the algorithm until the time limit is reached. Same is done for both ND and TS algorithms we introduce next.


## 2.2.2. Next Descent Algorithm

As mentioned in Section2.2.1, the steepest descent might take too long to converge if the size of the neighbourhood is too large or solving an individual subproblem is too time-consuming. Thus, we include in this paper a local search algorithm called next descent that aims to converge faster than the steepest descent, without a major impact on the solution quality. Algorithm 2 shows the next descent method.Algorithm 2:Next descent algorithm.
Input: α (portion of the resources included in I)
Output:  (x^,y^) (locally optimal solution)
(1) begin
(2)  k=0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (x^,y^)= solve MIPIk(x,y);
(5)   repeat
(6)   localOptimum = true;
(7)   k=k+1;
(8)   foreach  Ik-1′∈N(Ik-1)do
(9)     (x′,y′)= solve MIPIk-1′(x,y);
(10)     if  MIP(x′,y′)<MIP(x^,y^)  then
(11)      (x^,y^)=(x′,y′);
(12)      Ik=Ik-1′;
(13)      localOptimum = false;
(14)      break;
(15)  until localOptimum;
(16)  return (x^,y^);Just as in the steepest descent algorithm, the next descent algorithm starts with an initial solution, which is labelled as the current solution. Then, a random element from the neighbourhood of the current solution is selected and solved, and the objective function value of its optimal solution is compared to the objective function value of the current solution. Like in the SD algorithm, if the neighbour solution is not better than the current solution, another randomly selected setI′ from the neighbourhood N(Ik-1) is generated and compared to the current solution. Unlike the SD algorithm, in the ND algorithm, if the neighbour is better than the current solution, then the neighbour is labelled as the new current solution and no other I′∈N(Ik-1) is visited. The algorithm repeats these steps until the entire neighbourhood has been computed with no neighbour resulting in a better solution than the current one, in which case the algorithm stops.


## 2.2.3. Tabu Search Algorithm

Unlike the algorithms described in the previous sections, tabu search is a local search technique guided by the use of adaptive or flexible memory structures [5], and thus, it is inherently a stochastic local search algorithm. The variety of the tools and search principles introduced and described in [21] are such that the TS can be considered as the seed of a general framework for modern heuristic search [22]. We include TS as it has been applied to several combinatorial optimisation problems (see, e.g., [5, 23–27]) including, of course, mixed integer programming problems as the one we consider in this study.As in the algorithms introduced above, TS also starts with an initial set of constraintsI0 for which its associated subproblem MIPI0(x,y) is solved. Then, as in the SD algorithm, the “entire” neighbourhood of the current solution is computed. As explained before, since the neighbourhood N has a stochastic component, it is actually not possible to generate the entire neighbourhood of a solution; however, we use the term “entire” to stress the fact that all thegenerated neighbours of N are solved. This is different from the ND algorithm, where not all the generated neighbours are necessarily solved. The number of generated neighbours, ns, is, as in the algorithms above, set equal to 10. Moreover, TS implements a list calledtabu list that aims to avoid cycles during the search. Each time a neighbourhood is generated, the warehouses removed from I(k-1) are marked as tabu. Once we have solved the subproblems MIPIk′(x,y), with Ik′∈N(Ik), its optimal solutions are ranked and the one with the best objective function value is chosen as candidate solution for the next iteration. If the candidate solution was generated using a movement within the tabu list, it should be discarded, and the next best neighbour of the list should be chosen as the new candidate solution. However, there is one exception to this rule: if the candidate solution is within the tabu list but its objective function value is better than the best objective function found so far by the algorithm, then the so calledaspiration criterion is invoked and the candidate solution is set as the new current solution and passed on to the next iteration.Unlike the algorithms introduced above, TS does not require next current solution to be better than the previous one. This means that it is able to avoid local optimal by choosing solutions that are more expensive as they allows the algorithm to visit other (hopefully promising) areas in the search space. However, in case the algorithm cannot make any improvement after a predefined number of iterations, a diversification mechanism is used to get out from low-quality neighbourhoods and “jump” to other neighbourhoods. The diversification mechanism implemented here is a restart method, which set the current solution to a randomly generated solution without losing the best solution found so far. Termination criterion implemented here is the time limit.The tabu search implemented in this paper is as follows.As Algorithm3 shows, tabu search requires the following parameters:(i)
Time limit: total time the algorithm will perform.(ii)
DiversBound: total number of iteration without improvements on the best solution before diversification criterion (restart method) is applied.(iii)
TabuListSize(ts): size oftabuList. Number of iterations for which a specific movement remains banned.Algorithm 3:Tabu search algorithm.
Input: α (portion of the resources included in Ik)
Output:  (x^,y^) (approximately optimal solution)
(1) begin
(2)  k=0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (xc,yc)= solve MIPk(x,y);   // set current sol
(5)  (x^,y^)=(xc,yc); //set best sol
(6)  Y=∅; // tabu list initially empty
(7)  noImprovementCounter=0;
(8)  repeat
(9)    C=∅;
(10)   k=k+1;
(11)   foreachIk-1′∈N(Ik-1′)   do
(12)    (xn,yn)= solve MIPIk-1′(x,y);
(13)    C=C∪xn,yn;
(14)   SortC;
(15)   foreach  {(xn,yn)}∈C  do
(16)    if  isTabu((xn,yn))  then
(17)     if  MIP(xn,yn)<MIP(x^,y^)  then
(18)      Y=updateTabuList(Y,(xc,yc),(xn,yn),ts);
(19)      (xc,yc)=(xn,yn);
(20)      (x^,y^)=(xn,yn);
(21)      Ik=Ik-1n;
(22)      noImprovementCounter=0;
(23)      break;
(24)    else
(25)     Y=updateTabuList(Y,(xc,yc),(xn,yn));
(26)     (xc,yc)=(xn,yn);
(27)     Ik=Ik-1n  if  MIP(xn,yn)<MIP(x^,y^)  then
(28)       (x^,y^)=(xn,yn);
(29)       noImprovementCounter=0;
(30)     else
(31)       noImprovementCounter++;
(32)     break;
(33)   if  noImprovementCounter≥diversBound  then
(34)    Ik=selectBinaryVariablesRandomly(α);
(35)  until  time limit is reached;
(36)  return  (x^,y^);


## 2.2.4. Blind Algorithm

We finally implement a very simple heuristic method we callblind as it moves randomly at each iteration. Pseudocode of this method is presented in Algorithm 4.Algorithm 4:Matheuristic method using the blind algorithm.
// To find an (approximately) optimal solution for MIP
(
x
,
y
)
Input:  α (portion of the resources included in Ik)
Output:   (x^,y^) (approximately optimal solution)
(1) begin
(2)  k = 0;
(3)  Ik=selectBinaryVariablesRandomly(α);
(4)  (x^,y^) = solve (MIPIk(x,y));
(5)  while  !stopCriterion  do
(6)   k=k+1;
(7)   Ik=selectBinaryVariablesRandomly(α);
(8)   (xc,yc) = solve (MIPIk(x,y));
(9)   if  MIP(xc,yc)<MIP(x^,y^)  then
(10)    (x^,y^) = xc,yc;
(11)    Ik=Ik-1
(12)  return (x^,y^);As we can see, no additional intelligence is added to the blind algorithm. It is just a random search that, after a predefined number of iterations (or any other “stop criterion”), returns the best solution it found during its search. Thus, we can consider this algorithm as a baseline of this study.We apply the all four algorithms described in this section to two prostate cases. Details on this case and the obtained results are presented in the next section.


## 2.3. Subproblem Sizing

All the algorithms above perform very different as the value of the input parameterα varies. On the one hand, setting α to a very small value provokes that either the solver fails to solve the subproblem because of lack of memory or it takes too long to find an optimal solution of the subproblem. On the other hand, setting α close to the total number of binary decision variables provokes that the algorithm fails to find a solution for the generated subproblems as there is no feasible solution to it (i.e., the subproblems are too restrictive). Thus, we have to find a value of α such that exact methods can solve the obtained subproblems within an acceptable time. Further, α should allow local search methods to iterate as much as needed. Solving the obtained subproblems within few seconds is critical to matheuristic methods as they usually need several iterations before to converge to a good quality solution. This is especially true for matheuristics that consider local search or population-based heuristic methods. Then, there is a trade-off between the quality of the solution of subproblems and the time that is needed to generate such solutions. Therefore, finding a value of α that gives us a good compromise between these two aspects is critical for the overall performance of the matheuristic methods explained before. It is interesting that the problem of findingefficient values of α might be seen as a multiobjective optimisation problem.In next section we explain the experiments that we perform to study how the choice ofα hits local-search-based matheuristics performance. Based on the results, we draft some guidelines to set value of parameter α at the end of next section.


## 3. Computational Experiments

This section starts briefly introducing the problem we consider in this paper. Then, the experiments performed here are presented and their results are discussed.


### 3.1. The Capacitated Facility Location Problem

The capacitated facility location problem (CFLP) is a well-known problem in combinatorial optimisation. The CFLP has been shown to be NP-hard [28]. The problem consists of selecting specific sites at which to install plants, warehouses, and distribution centres while assigning customers to service facilities and interconnecting facilities using flow assignment decisions [23]. In this study we consider the CFLP problem to evaluate the performance of a very simple matheuristic algorithm. We consider a two-level supply chain in which a single plant serves a set of warehouses, which in turn serve a set of end customers or retailers. Figure 2 shows the basic configuration of our supply chain. Thus, the goal is to find a set of locations that serves the entire set of customers in an optimal way. As Figure 2 shows, each customer (or cluster) is served only by one warehouse.Figure 2
Distribution network structure considered in this study. It consists of one central plant, a set of potential warehouses, and a set of customers or retailers [23].The optimisation model considers theinstallation cost (i.e., the cost associated with opening a specific warehouse) and transportation orallocation cost (i.e., the cost of transporting one item from a warehouse to a customer). The mathematical model for the CFLP is(4a)CFLPx,y:minx,y∑i=1nfiyi+∑i=1n∑j=1mcijdjxij(4b)s.t.∑j=1mdjxij≤Iicapyi∀i=1,…,n,(4c)∑i=1nxij=1∀j=1,…,m,(4d)xij≤yi∀i=1,…,n;∀j=1,…,m,(4e)yi,xij∈0,1∀i=1,…,n;∀j=1,…,m.Equation (4a) is the total system cost. The first term is the fixed setup and operating cost when opening warehouses. The second term is the daily transport cost between warehouse and customers which depends on the customer demand d and distance cij between warehouse i and customer j. Inequality (4b) ensures that total demand of warehouse i will never be greater than its capacity Iicap. Equation (4c) ensures that customers are served by only one warehouse. Equation (4d) makes sure that customers are only allocated to available warehouses. Finally, (4e) states integrality (0-1) for the binary variables xij and yi. Other versions of the CFLP relax this constraints by making x∈Rn×m with 0≤xij≤1. In that case, the CFLP would be just as the MIP(x,y) problem.


### 3.2. Experiments

In this paper three benchmarks for the CFLP are considered. The first benchmark corresponds to problem setsA,  B, and C from the OR Library [29]. Instances in problem sets A, B, and C consider 1,000 clients and 100 warehouses. Warehouses capacity for instances A1, A2, A3, and A4 are equal to 8,000; 10,000; 12,000; and 14,000, respectively. Warehouses capacity for instances B1, B2, B3, and B4 are equal to 5,000; 6,000; 7,000; and 8,000, respectively. Warehouses capacity for instances C1, C2, C3, and C4 are equal to 5,000; 5,750; 6,500; and 7,250, respectively. Table 1 shows these instances and their corresponding optimal values obtained by the MILP solver. Column t shows the time needed by the MILP solver to reach the optimal solution.Table 1
Instances for the first benchmark used in this study (OR Library).
OR Library
Instances
|
I
|
×
|
J
|
I
c
a
p
f
(
x
∗
)
t (sec)
A
1
1
,
000
×
100
8
,
000
19,240,822.45
1,745
A
2
1
,
000
×
100
10
,
000
18,438,046.54
1,778
A
3
1,000
×
100
12
,
000
17,765,201.95
199
A
4
1
,
000
×
100
14
,
000
17,160,439.01
8
B
1
1
,
000
×
100
5
,
000
13,656,379.58
136
B
2
1
,
000
×
100
6
,
000
13,361,927.45
151
B
3
1
,
000
×
100
7
,
000
13,198,556.43
256
B
4
1
,
000
×
100
8
,
000
13,082,516.50
27
C
1
1
,
000
×
100
5
,
000
11,646,596.97
155
C
2
1
,
000
×
100
5
,
750
11,570,340.29
70
C
3
1
,
000
×
100
6
,
500
11,518,743.74
31
C
4
1
,
000
×
100
7
,
250
11,505,767.39
24The second benchmark is a set of instances where clients and warehouses are uniformly distributed over an imaginary square of100×100[distanceunits]2 (see Figure 3(a)). We call this set DS1U. The number of clients considered in instances belonging to set DS1U ranges from 500 to 1000 (500, 600, 700, 800, 900, and 1,000) while the number of warehouses considered varies between 500 and 1,000. Thus, set DS1U consists of 12 problem classes (500×500,500×600,…,500×900,500×1000,1000×500,…,1000×1000). For each problem class, 10 instances are randomly generated using the procedure proposed in [30] and that was also used in [5, 23]. We do this in order to minimise any instance dependant effect. Table 2 shows the average values for each class of problems.Table 2
Instances for theDS1U and DS1C benchmarks.
Instances
D
S
1
_
U
D
S
1
_
C
|
I
|
×
|
J
|
f
(
x
∗
)
t
a
v
g
t
m
i
n
t
m
a
x
f
(
x
∗
)
t
a
v
g
t
m
i
n
t
m
a
x
(1)
500
×
500
120670.3
170
95
287
123900.7
191
99
410
(2)
500
×
1000
110534.5
267
194
356
113527.1
262
201
410
(3)
600
×
500
149539.9
261
105
424
149343.2
287
120
522
(4)
600
×
1000
141847.7
399
309
568
137451.4
331
255
501
(5)
700
×
500
182434.9
427
128
674
181149.3
292
96
645
(6)
700
×
1000
161882.7
726
263
1515
161872.9
596
306
1673
(7)
800
×
500
209319.5
455
227
748
207036.9
343
131
648
(8)
800
×
1000
193647.7
761
399
1679
188859.7
1040
590
1828
(9)
900
×
500
243403.7
391
179
681
236287.7
428
178
1024
(10)
900
×
1000
211302.4
1228
672
1787
215410.7
1158
507
1564
(11)
1000
×
500
270449.9
448
189
1484
265805.9
407
154
642
(12)
1000
×
1000
243625.0
1187
670
1791
239519.8
1419
588
2653Figure 3
Example of the instances considered in this study.
(a)
Clients distributed in clusters
(b)
Clients uniformly distributedFinally, a third benchmark consisting on clients that are organised in clusters is considered. We call this benchmarkDS1C (see Figure 3(b)). Instances in DS1C are generated very similar to the ones in DS1U. The only difference is that clients in DS1C are not uniformly distributed and, thus, those warehouses that are within (or very close to) a cluster have an installation cost slightly higher than those warehouses that are far away from the clusters.Table2 shows the results obtained by the MIP solver from Gurobi when solving each instance of DS1C. As we can see, the MIP solver is able to solve all the instances to optimality. Further, the solver finds the optimal solution for almost all the instances in less than 1,800 secs. Columns tavg, tmin, and tmax in Table 2 show the average time and both minimum and maximum times, respectively. This is because we solve 10 different instances for each instance class. As mentioned before, we do this to avoid any instance dependent effect.After we have solved the problem using the MIP solver, we apply the all three local-search-based matheuristics and the blind algorithm to each instance and allow them to run for 2,000 secs. The proposed algorithms solve each instance 10 times for each value of(1-α). Table 3 shows the results obtained by the local-search-based matheuristic methods for instances A, B, and C for four different values of (1-α). As expected, the blind algorithm consistently obtains GAP values much higher than the other three methods. Further, local-search-based algorithms are able to find the optimal solution for the majority of the instances. While the average GAP for the blind algorithm is 1.65%, the average GAP for the SD, ND, and TS methods is 0.33%, 0.36%, and 0.15%, respectively. As expected, the best average value is obtained by the TS algorithm while difference between SD and ND algorithms is negligible. Regarding the time needed by the algorithms to converge, the blind algorithm is the one that takes longer (770 secs). As mentioned before, the ND algorithm is, in average, the fastest one with 299 seconds, while the SD algorithm almost doubles this time with 547 seconds, in average. TS algorithm takes 371 seconds, in average, before convergence.Table 3
Results obtained by the local search methods for the OR Library’s instances.
Inst
(
1
-
α
)
Blind
SD
ND
TS
GAP
Time
GAP
Time
GAP
Time
GAP
Time
A
1
0.1
5,22%
622
9,65%
248
9,45%
1424
6,73%
785
0.3
1,63%
1124
0,00%
1240
0,00%
489
0,00%
736
0.5
0,32%
664
0,00%
1095
0,00%
274
0,24%
700
0.7
0,10%
745
0,00%
366
0,00%
138
0,00%
1476
A
2
0.1
4,31%
463
0,00%
833
0,00%
1415
0,00%
866
0.3
1,22%
893
0,00%
366
0,19%
1290
0,19%
615
0.5
0,52%
987
1,83%
52
0,09%
956
0,00%
745
0.7
0,11%
282
0,09%
378
0,00%
1562
0,00%
743
A
3
0.1
1,81%
841
0,00%
341
0,31%
65
0,00%
686
0.3
0,23%
503
0,00%
602
0,00%
693
0,00%
1083
0.5
0,12%
880
0,00%
431
0,00%
442
0,00%
659
0.7
0,00%
272
0,00%
151
0,00%
211
0,00%
211
A
4
0.1
1,56%
990
0,00%
25
0,00%
84
0,00%
84
0.3
0,00%
411
0,00%
8
0,00%
31
0,00%
31
0.5
0,00%
111
0,00%
2
0,00%
14
0,00%
14
0.7
0,00%
23
0,00%
1
0,00%
7
0,00%
7
B
1
0.1
5,22%
622
0,03%
639
0,23%
960
0,01%
319
0.3
1,63%
1124
0,01%
588
0,08%
436
0,03%
1452
0.5
0,32%
664
0,00%
707
0,03%
216
0,00%
605
0.7
0,10%
745
0,01%
347
0,00%
1630
0,00%
1192
B
2
0.1
4,31%
463
1,44%
806
2,75%
222
0,09%
714
0.3
1,22%
893
0,00%
135
0,00%
69
0,00%
1165
0.5
0,52%
987
0,00%
186
0,00%
196
0,05%
530
0.7
0,11%
282
0,00%
112
0,05%
315
0,00%
261
B
3
0.1
1,81%
841
0,40%
403
0,43%
1270
0,30%
142
0.3
0,23%
503
0,00%
96
0,00%
181
0,00%
267
0.5
0,12%
880
0,00%
588
0,00%
236
0,00%
564
0.7
0,00%
272
0,00%
179
0,00%
316
0,00%
316
B
4
0.1
1,56%
990
0,00%
1309
0,32%
144
0,32%
144
0.3
0,00%
411
0,00%
24
0,00%
360
0,00%
360
0.5
0,00%
111
0,00%
57
0,00%
381
0,00%
381
0.7
0,00%
23
0,00%
28
0,00%
349
0,00%
349
C
1
0.1
8,09%
712
9,70%
1220
8,14%
92
10,27%
133
0.3
3,20%
919
0,00%
657
0,00%
486
0,00%
394
0.5
1,99%
1105
0,00%
173
0,00%
496
0,00%
1469
0.7
0,48%
863
0,00%
239
0,00%
355
0,00%
247
C
2
0.1
5,74%
937
1,42%
1254
2,32%
1397
0,70%
84
0.3
1,64%
934
0,00%
136
0,00%
107
0,00%
172
0.5
0,35%
562
0,00%
1225
0,00%
248
0,00%
537
0.7
0,02%
766
0,00%
1294
0,00%
164
0,00%
169
C
3
0.1
5,17%
979
1,11%
1254
1,54%
819
0,20%
585
0.3
1,04%
993
0,00%
73
0,00%
38
0,00%
272
0.5
0,22%
906
0,00%
587
0,00%
64
0,00%
1387
0.7
0,04%
740
0,00%
551
0,00%
149
0,00%
779
C
4
0.1
4,86%
663
1,35%
153
0,34%
523
0,85%
313
0.3
0,65%
711
0,03%
7
0,03%
35
0,03%
107
0.5
0,08%
821
0,00%
18
0,03%
11
0,03%
13
0.7
0,02%
228
0,03%
8
0,03%
22
0,00%
29Figures4(a), 5(a), and 6(a) show the evolution of the GAP, as parameter (1-α) increases. As we expected, the larger the value of (1-α) is set to, the smaller the GAP. Also, it is clear that the blind algorithms obtain higher GAP values than the other three algorithms considered in this study. As we mentioned before, the blind algorithm works as our baseline algorithm. It is interesting to note that, for all algorithms, the worst performance is obtained for (1-α)=0.1. This is mainly because for this value not enough facilities are available and the algorithms open a set of facilities randomly so the problem has a feasible solution. However, this repairing process is completely random and, thus, algorithms performance is impaired.Figure 4
Average results obtained by the local search algorithms for each value of(1-α) for ORLib  A instances.
(a)
Evolution of the GAP forORLib  A instances
(b)
Evolution of the time forORLib  A instancesFigure 5
Average results obtained by the local search algorithms for each value of(1-α) for ORLib  B instances.
(a)
Evolution of the GAP forORLib  B instances
(b)
Evolution of the time forORLib  B instancesFigure 6
Average results obtained by the local search algorithms for each value of(1-α) for ORLib  C instances.
(a)
Evolution of the GAP forORLib  C  instances
(b)
Evolution of the time forORLib    C  instancesSince the problems from the OR Library are only medium size instances, local-search-based matheuristics consistently find the optimal solution for almost all instances for(1-α)={0.3,0.5,0.7}. In fact, even the blind algorithm finds solutions that are very close to the optimal ones when (1-α)=0.7.Figures4(b), 5(b), and 6(b) show the time needed by our algorithms to find its best solution, as parameter (1-α) increases. Unlike we expected, the time needed to find the best solution does not increase as the parameter (1-α) gets larger. In fact, both algorithms, SD and blind, converge faster when (1-α=0.7), the larger value we tried in our experiments. This can be explained because of the problems features. We note that optimal solutions for all the problems in the benchmark of the OR Library need not too many warehouses to be open. Thus, when a large portion of the potential warehouses are available (as for the (1-α)=0.7 case) the algorithm is likely to find such optimal solutions in early iterations. In fact, for the vast majority of the experiments, optimal solution is found within the first 2 to 3 iterations.We now move on theDS1U and the DS1C benchmarks. As we noted before, these sets of instances are much larger and harder to solve than the problems in the OR Library we discussed above. Further, as the optimal solutions of these instances do not require too many warehouses to be open, the repairing procedure applied to the OR Library instances for (1-α)=0.1 does not have a great impact on the results. Thus, local-search-based matheuristics outperform the results obtained by the blind algorithm for all values of parameter α. Figures 7(a) and 7(b) show the GAP values obtained by the all four algorithms for instances DS1U and DS1C, respectively.Figure 7
Average GAP values obtained by the local search algorithms for each value of(1-α) for DS1 instances.
(a)
Evolution of the GAP forDS1U instances
(b)
Evolution of the GAP forDS1C instancesJust as in the OR Library instances, as the parameter(1-α) gets larger, the GAP approximates to 0 for the all three local-search-based matheuristic proposed in this paper. While for the DS1U instances both the TS and the ND algorithms reach GAP values very close to 0 for (1-α)=0.5 and (1-α)=0.7, for the DS1C instances the best value obtained by the TS and the ND algorithms is 0.71 and 0.61, respectively, for (1-α)=0.7. Notably, the ND algorithm performs slightly better than the TS for the DS1C instance when (1-α)=0.5, obtaining a best average value of 0.66 against the 1.42 obtained by the TS algorithm. In spite of that, we can note that the differences between the GAP obtained by both the ND and the TS algorithms when using (1-α)=0.5 and (1-α)=0.7 are negligible, just as in the OR Library instances.Figures8(a) and 8(b) shows the time needed by the all four algorithms to find their best solution for instances DS1U and DS1C, respectively. As expected, as the parameter (1-α) gets larger, longer times are needed by the algorithms to converge. This is important as one would like to use a value for (1-α) such that minimum GAP values are reached within few minutes.Figure 8
Average times needed by the local search algorithms for each value of(1-α) for DS1 instances.
(a)
Evolution of the time forDS1U instances
(b)
Evolution of the time forDS1C instancesResults in Figures7 and 8 show, clearly, that there is a compromise between the quality of the solution found by the algorithms and the time they need to do so. Further, we know that this compromise can be managed by adjusting the value of parameter (1-α). For the case of instances DS1U and DS1C, such a value should be between 0.4 and 0.6. Moreover, for the case of the instances of the OR Library, parameter (1-α) should be set around 0.3. It is interesting to note that while the size of instances DS1U and DS1C is equivalent, instances of the OR Library are much smaller.We can also note that instances that include clusters tend to take longer to converge. This is especially true as parameter(1-α) gets larger. Further, for these instances, as the parameter (1-α) gets larger, less iterations are performed by algorithms although each of these iterations takes much longer. As mentioned before, these should be taken into account when using population-based algorithms within the framework presented in this study as such kind of heuristic algorithms needs to perform several iterations before to converge to good quality solutions.


## 3.1. The Capacitated Facility Location Problem

The capacitated facility location problem (CFLP) is a well-known problem in combinatorial optimisation. The CFLP has been shown to be NP-hard [28]. The problem consists of selecting specific sites at which to install plants, warehouses, and distribution centres while assigning customers to service facilities and interconnecting facilities using flow assignment decisions [23]. In this study we consider the CFLP problem to evaluate the performance of a very simple matheuristic algorithm. We consider a two-level supply chain in which a single plant serves a set of warehouses, which in turn serve a set of end customers or retailers. Figure 2 shows the basic configuration of our supply chain. Thus, the goal is to find a set of locations that serves the entire set of customers in an optimal way. As Figure 2 shows, each customer (or cluster) is served only by one warehouse.Figure 2
Distribution network structure considered in this study. It consists of one central plant, a set of potential warehouses, and a set of customers or retailers [23].The optimisation model considers theinstallation cost (i.e., the cost associated with opening a specific warehouse) and transportation orallocation cost (i.e., the cost of transporting one item from a warehouse to a customer). The mathematical model for the CFLP is(4a)CFLPx,y:minx,y∑i=1nfiyi+∑i=1n∑j=1mcijdjxij(4b)s.t.∑j=1mdjxij≤Iicapyi∀i=1,…,n,(4c)∑i=1nxij=1∀j=1,…,m,(4d)xij≤yi∀i=1,…,n;∀j=1,…,m,(4e)yi,xij∈0,1∀i=1,…,n;∀j=1,…,m.Equation (4a) is the total system cost. The first term is the fixed setup and operating cost when opening warehouses. The second term is the daily transport cost between warehouse and customers which depends on the customer demand d and distance cij between warehouse i and customer j. Inequality (4b) ensures that total demand of warehouse i will never be greater than its capacity Iicap. Equation (4c) ensures that customers are served by only one warehouse. Equation (4d) makes sure that customers are only allocated to available warehouses. Finally, (4e) states integrality (0-1) for the binary variables xij and yi. Other versions of the CFLP relax this constraints by making x∈Rn×m with 0≤xij≤1. In that case, the CFLP would be just as the MIP(x,y) problem.


## 3.2. Experiments

In this paper three benchmarks for the CFLP are considered. The first benchmark corresponds to problem setsA,  B, and C from the OR Library [29]. Instances in problem sets A, B, and C consider 1,000 clients and 100 warehouses. Warehouses capacity for instances A1, A2, A3, and A4 are equal to 8,000; 10,000; 12,000; and 14,000, respectively. Warehouses capacity for instances B1, B2, B3, and B4 are equal to 5,000; 6,000; 7,000; and 8,000, respectively. Warehouses capacity for instances C1, C2, C3, and C4 are equal to 5,000; 5,750; 6,500; and 7,250, respectively. Table 1 shows these instances and their corresponding optimal values obtained by the MILP solver. Column t shows the time needed by the MILP solver to reach the optimal solution.Table 1
Instances for the first benchmark used in this study (OR Library).
OR Library
Instances
|
I
|
×
|
J
|
I
c
a
p
f
(
x
∗
)
t (sec)
A
1
1
,
000
×
100
8
,
000
19,240,822.45
1,745
A
2
1
,
000
×
100
10
,
000
18,438,046.54
1,778
A
3
1,000
×
100
12
,
000
17,765,201.95
199
A
4
1
,
000
×
100
14
,
000
17,160,439.01
8
B
1
1
,
000
×
100
5
,
000
13,656,379.58
136
B
2
1
,
000
×
100
6
,
000
13,361,927.45
151
B
3
1
,
000
×
100
7
,
000
13,198,556.43
256
B
4
1
,
000
×
100
8
,
000
13,082,516.50
27
C
1
1
,
000
×
100
5
,
000
11,646,596.97
155
C
2
1
,
000
×
100
5
,
750
11,570,340.29
70
C
3
1
,
000
×
100
6
,
500
11,518,743.74
31
C
4
1
,
000
×
100
7
,
250
11,505,767.39
24The second benchmark is a set of instances where clients and warehouses are uniformly distributed over an imaginary square of100×100[distanceunits]2 (see Figure 3(a)). We call this set DS1U. The number of clients considered in instances belonging to set DS1U ranges from 500 to 1000 (500, 600, 700, 800, 900, and 1,000) while the number of warehouses considered varies between 500 and 1,000. Thus, set DS1U consists of 12 problem classes (500×500,500×600,…,500×900,500×1000,1000×500,…,1000×1000). For each problem class, 10 instances are randomly generated using the procedure proposed in [30] and that was also used in [5, 23]. We do this in order to minimise any instance dependant effect. Table 2 shows the average values for each class of problems.Table 2
Instances for theDS1U and DS1C benchmarks.
Instances
D
S
1
_
U
D
S
1
_
C
|
I
|
×
|
J
|
f
(
x
∗
)
t
a
v
g
t
m
i
n
t
m
a
x
f
(
x
∗
)
t
a
v
g
t
m
i
n
t
m
a
x
(1)
500
×
500
120670.3
170
95
287
123900.7
191
99
410
(2)
500
×
1000
110534.5
267
194
356
113527.1
262
201
410
(3)
600
×
500
149539.9
261
105
424
149343.2
287
120
522
(4)
600
×
1000
141847.7
399
309
568
137451.4
331
255
501
(5)
700
×
500
182434.9
427
128
674
181149.3
292
96
645
(6)
700
×
1000
161882.7
726
263
1515
161872.9
596
306
1673
(7)
800
×
500
209319.5
455
227
748
207036.9
343
131
648
(8)
800
×
1000
193647.7
761
399
1679
188859.7
1040
590
1828
(9)
900
×
500
243403.7
391
179
681
236287.7
428
178
1024
(10)
900
×
1000
211302.4
1228
672
1787
215410.7
1158
507
1564
(11)
1000
×
500
270449.9
448
189
1484
265805.9
407
154
642
(12)
1000
×
1000
243625.0
1187
670
1791
239519.8
1419
588
2653Figure 3
Example of the instances considered in this study.
(a)
Clients distributed in clusters
(b)
Clients uniformly distributedFinally, a third benchmark consisting on clients that are organised in clusters is considered. We call this benchmarkDS1C (see Figure 3(b)). Instances in DS1C are generated very similar to the ones in DS1U. The only difference is that clients in DS1C are not uniformly distributed and, thus, those warehouses that are within (or very close to) a cluster have an installation cost slightly higher than those warehouses that are far away from the clusters.Table2 shows the results obtained by the MIP solver from Gurobi when solving each instance of DS1C. As we can see, the MIP solver is able to solve all the instances to optimality. Further, the solver finds the optimal solution for almost all the instances in less than 1,800 secs. Columns tavg, tmin, and tmax in Table 2 show the average time and both minimum and maximum times, respectively. This is because we solve 10 different instances for each instance class. As mentioned before, we do this to avoid any instance dependent effect.After we have solved the problem using the MIP solver, we apply the all three local-search-based matheuristics and the blind algorithm to each instance and allow them to run for 2,000 secs. The proposed algorithms solve each instance 10 times for each value of(1-α). Table 3 shows the results obtained by the local-search-based matheuristic methods for instances A, B, and C for four different values of (1-α). As expected, the blind algorithm consistently obtains GAP values much higher than the other three methods. Further, local-search-based algorithms are able to find the optimal solution for the majority of the instances. While the average GAP for the blind algorithm is 1.65%, the average GAP for the SD, ND, and TS methods is 0.33%, 0.36%, and 0.15%, respectively. As expected, the best average value is obtained by the TS algorithm while difference between SD and ND algorithms is negligible. Regarding the time needed by the algorithms to converge, the blind algorithm is the one that takes longer (770 secs). As mentioned before, the ND algorithm is, in average, the fastest one with 299 seconds, while the SD algorithm almost doubles this time with 547 seconds, in average. TS algorithm takes 371 seconds, in average, before convergence.Table 3
Results obtained by the local search methods for the OR Library’s instances.
Inst
(
1
-
α
)
Blind
SD
ND
TS
GAP
Time
GAP
Time
GAP
Time
GAP
Time
A
1
0.1
5,22%
622
9,65%
248
9,45%
1424
6,73%
785
0.3
1,63%
1124
0,00%
1240
0,00%
489
0,00%
736
0.5
0,32%
664
0,00%
1095
0,00%
274
0,24%
700
0.7
0,10%
745
0,00%
366
0,00%
138
0,00%
1476
A
2
0.1
4,31%
463
0,00%
833
0,00%
1415
0,00%
866
0.3
1,22%
893
0,00%
366
0,19%
1290
0,19%
615
0.5
0,52%
987
1,83%
52
0,09%
956
0,00%
745
0.7
0,11%
282
0,09%
378
0,00%
1562
0,00%
743
A
3
0.1
1,81%
841
0,00%
341
0,31%
65
0,00%
686
0.3
0,23%
503
0,00%
602
0,00%
693
0,00%
1083
0.5
0,12%
880
0,00%
431
0,00%
442
0,00%
659
0.7
0,00%
272
0,00%
151
0,00%
211
0,00%
211
A
4
0.1
1,56%
990
0,00%
25
0,00%
84
0,00%
84
0.3
0,00%
411
0,00%
8
0,00%
31
0,00%
31
0.5
0,00%
111
0,00%
2
0,00%
14
0,00%
14
0.7
0,00%
23
0,00%
1
0,00%
7
0,00%
7
B
1
0.1
5,22%
622
0,03%
639
0,23%
960
0,01%
319
0.3
1,63%
1124
0,01%
588
0,08%
436
0,03%
1452
0.5
0,32%
664
0,00%
707
0,03%
216
0,00%
605
0.7
0,10%
745
0,01%
347
0,00%
1630
0,00%
1192
B
2
0.1
4,31%
463
1,44%
806
2,75%
222
0,09%
714
0.3
1,22%
893
0,00%
135
0,00%
69
0,00%
1165
0.5
0,52%
987
0,00%
186
0,00%
196
0,05%
530
0.7
0,11%
282
0,00%
112
0,05%
315
0,00%
261
B
3
0.1
1,81%
841
0,40%
403
0,43%
1270
0,30%
142
0.3
0,23%
503
0,00%
96
0,00%
181
0,00%
267
0.5
0,12%
880
0,00%
588
0,00%
236
0,00%
564
0.7
0,00%
272
0,00%
179
0,00%
316
0,00%
316
B
4
0.1
1,56%
990
0,00%
1309
0,32%
144
0,32%
144
0.3
0,00%
411
0,00%
24
0,00%
360
0,00%
360
0.5
0,00%
111
0,00%
57
0,00%
381
0,00%
381
0.7
0,00%
23
0,00%
28
0,00%
349
0,00%
349
C
1
0.1
8,09%
712
9,70%
1220
8,14%
92
10,27%
133
0.3
3,20%
919
0,00%
657
0,00%
486
0,00%
394
0.5
1,99%
1105
0,00%
173
0,00%
496
0,00%
1469
0.7
0,48%
863
0,00%
239
0,00%
355
0,00%
247
C
2
0.1
5,74%
937
1,42%
1254
2,32%
1397
0,70%
84
0.3
1,64%
934
0,00%
136
0,00%
107
0,00%
172
0.5
0,35%
562
0,00%
1225
0,00%
248
0,00%
537
0.7
0,02%
766
0,00%
1294
0,00%
164
0,00%
169
C
3
0.1
5,17%
979
1,11%
1254
1,54%
819
0,20%
585
0.3
1,04%
993
0,00%
73
0,00%
38
0,00%
272
0.5
0,22%
906
0,00%
587
0,00%
64
0,00%
1387
0.7
0,04%
740
0,00%
551
0,00%
149
0,00%
779
C
4
0.1
4,86%
663
1,35%
153
0,34%
523
0,85%
313
0.3
0,65%
711
0,03%
7
0,03%
35
0,03%
107
0.5
0,08%
821
0,00%
18
0,03%
11
0,03%
13
0.7
0,02%
228
0,03%
8
0,03%
22
0,00%
29Figures4(a), 5(a), and 6(a) show the evolution of the GAP, as parameter (1-α) increases. As we expected, the larger the value of (1-α) is set to, the smaller the GAP. Also, it is clear that the blind algorithms obtain higher GAP values than the other three algorithms considered in this study. As we mentioned before, the blind algorithm works as our baseline algorithm. It is interesting to note that, for all algorithms, the worst performance is obtained for (1-α)=0.1. This is mainly because for this value not enough facilities are available and the algorithms open a set of facilities randomly so the problem has a feasible solution. However, this repairing process is completely random and, thus, algorithms performance is impaired.Figure 4
Average results obtained by the local search algorithms for each value of(1-α) for ORLib  A instances.
(a)
Evolution of the GAP forORLib  A instances
(b)
Evolution of the time forORLib  A instancesFigure 5
Average results obtained by the local search algorithms for each value of(1-α) for ORLib  B instances.
(a)
Evolution of the GAP forORLib  B instances
(b)
Evolution of the time forORLib  B instancesFigure 6
Average results obtained by the local search algorithms for each value of(1-α) for ORLib  C instances.
(a)
Evolution of the GAP forORLib  C  instances
(b)
Evolution of the time forORLib    C  instancesSince the problems from the OR Library are only medium size instances, local-search-based matheuristics consistently find the optimal solution for almost all instances for(1-α)={0.3,0.5,0.7}. In fact, even the blind algorithm finds solutions that are very close to the optimal ones when (1-α)=0.7.Figures4(b), 5(b), and 6(b) show the time needed by our algorithms to find its best solution, as parameter (1-α) increases. Unlike we expected, the time needed to find the best solution does not increase as the parameter (1-α) gets larger. In fact, both algorithms, SD and blind, converge faster when (1-α=0.7), the larger value we tried in our experiments. This can be explained because of the problems features. We note that optimal solutions for all the problems in the benchmark of the OR Library need not too many warehouses to be open. Thus, when a large portion of the potential warehouses are available (as for the (1-α)=0.7 case) the algorithm is likely to find such optimal solutions in early iterations. In fact, for the vast majority of the experiments, optimal solution is found within the first 2 to 3 iterations.We now move on theDS1U and the DS1C benchmarks. As we noted before, these sets of instances are much larger and harder to solve than the problems in the OR Library we discussed above. Further, as the optimal solutions of these instances do not require too many warehouses to be open, the repairing procedure applied to the OR Library instances for (1-α)=0.1 does not have a great impact on the results. Thus, local-search-based matheuristics outperform the results obtained by the blind algorithm for all values of parameter α. Figures 7(a) and 7(b) show the GAP values obtained by the all four algorithms for instances DS1U and DS1C, respectively.Figure 7
Average GAP values obtained by the local search algorithms for each value of(1-α) for DS1 instances.
(a)
Evolution of the GAP forDS1U instances
(b)
Evolution of the GAP forDS1C instancesJust as in the OR Library instances, as the parameter(1-α) gets larger, the GAP approximates to 0 for the all three local-search-based matheuristic proposed in this paper. While for the DS1U instances both the TS and the ND algorithms reach GAP values very close to 0 for (1-α)=0.5 and (1-α)=0.7, for the DS1C instances the best value obtained by the TS and the ND algorithms is 0.71 and 0.61, respectively, for (1-α)=0.7. Notably, the ND algorithm performs slightly better than the TS for the DS1C instance when (1-α)=0.5, obtaining a best average value of 0.66 against the 1.42 obtained by the TS algorithm. In spite of that, we can note that the differences between the GAP obtained by both the ND and the TS algorithms when using (1-α)=0.5 and (1-α)=0.7 are negligible, just as in the OR Library instances.Figures8(a) and 8(b) shows the time needed by the all four algorithms to find their best solution for instances DS1U and DS1C, respectively. As expected, as the parameter (1-α) gets larger, longer times are needed by the algorithms to converge. This is important as one would like to use a value for (1-α) such that minimum GAP values are reached within few minutes.Figure 8
Average times needed by the local search algorithms for each value of(1-α) for DS1 instances.
(a)
Evolution of the time forDS1U instances
(b)
Evolution of the time forDS1C instancesResults in Figures7 and 8 show, clearly, that there is a compromise between the quality of the solution found by the algorithms and the time they need to do so. Further, we know that this compromise can be managed by adjusting the value of parameter (1-α). For the case of instances DS1U and DS1C, such a value should be between 0.4 and 0.6. Moreover, for the case of the instances of the OR Library, parameter (1-α) should be set around 0.3. It is interesting to note that while the size of instances DS1U and DS1C is equivalent, instances of the OR Library are much smaller.We can also note that instances that include clusters tend to take longer to converge. This is especially true as parameter(1-α) gets larger. Further, for these instances, as the parameter (1-α) gets larger, less iterations are performed by algorithms although each of these iterations takes much longer. As mentioned before, these should be taken into account when using population-based algorithms within the framework presented in this study as such kind of heuristic algorithms needs to perform several iterations before to converge to good quality solutions.


## 4. Conclusions and Future Work

In this paper we show the impact of parameter tuning on a local-search-based matheuristic framework for solving mixed integer (non)linear programming problems. In particular, matheuristics that combine local search methods and a MIP solver are tested. In this study, we focus on the size of the subproblem generated by the local search method that is passed on to the MIP solver. As expected, the size of the subproblems that are solved in turn by the matheuristic method has a big impact on the behaviour of the matheuristic and, consequently, on its obtained results: as the size of the subproblem increases (i.e., more integer/binary decision variables are considered) the results obtained by the MIP solver are closer to the optimal solution. The time required by the algorithms tested in this paper to find its best solution also increases as the subproblem gets larger. Further, as the subproblem gets larger, fewer iterations can be performed within the allowed time. This is important as other heuristics such as evolutionary algorithms and swarm intelligence, where many iterations are needed before converging to a good quality solution, might be not able to deal with large subproblems. We also note that the improvement in the GAP values after certain value of parameter(1-α) is negligible and that this value depends to some extent on the size of the problem: for medium size instances such as the ones in the OR Library, parameter (1-α) should be set to a value around 0.3 while for larger instances (DS1U and DS1C) it should be set to a value between 0.4 and 0.6. The specific values will depend on both the accuracy level requested by users and the time available to perform the algorithm. Therefore, the challenge when designing matheuristic frameworks as the one presented in this paper is to find a value for parameter (1-α) that allows the heuristic algorithm to perform as many iterations as needed and that provides a good compromise between solution quality and run times. Although the values provided here are only valid for the problem and the algorithms considered in this paper, we think that the obtained results can be used as a guide by other researchers using similar frameworks and/or dealing with similar problems.As a future work, strategies such as evolutionary algorithms and swarm intelligence will be tested within the matheuristic framework considering the results obtained in this study. We expect that intelligent methods such as the ones named before greatly improve the results obtained by the local search methods considered in this study. Moreover, the matheuristic framework used in this paper might also be applied to other MILP and MINLP problems such as, for instance, the beam angle optimisation problem in radiation therapy for cancer treatment.



---

*Source: 1702506-2017-12-31.xml*