algorithm

Computers & Industrial Engineering 110 (2017) 216–228

Contents lists available at ScienceDirect

Computers & Industrial Engineering

journal homepage: www.elsevier.com/locate/caie

Biased randomization of heuristics using skewed probability distributions: A survey and some applications

http://dx.doi.org/10.1016/j.cie.2017.06.019 0360-8352/� 2017 Elsevier Ltd. All rights reserved.

⇑ Corresponding author. E-mail addresses: [email protected] (A. Grasas), [email protected] (A.A. Juan), javier.

[email protected] (J. Faulin), [email protected] (J. de Armas), helena. [email protected] (H. Ramalhinho).

Alex Grasas a, Angel A. Juan b,⇑, Javier Faulin c, Jesica de Armas b,d, Helena Ramalhinho d a Institute for Transformational Leadership, Spain b Open University of Catalonia, IN3 – Dept. of Computer Science, Spain c Public University of Navarre, Dept. of Statistics and OR, Spain d Univ. Pompeu Fabra, Dept. of Economics and Business, Spain

a r t i c l e i n f o a b s t r a c t

Article history: Received 15 June 2016 Received in revised form 24 January 2017 Accepted 10 June 2017 Available online 13 June 2017

Keywords: Heuristics Biased randomization Real-time decision making Combinatorial optimization Logistics Transportation Production

Randomized heuristics are widely used to solve large scale combinatorial optimization problems. Among the plethora of randomized heuristics, this paper reviews those that contain biased-randomized proce- dures (BRPs). A BRP is a procedure to select the next constructive ‘movement’ from a list of candidates in which their elements have different probabilities based on some criteria (e.g., ranking, priority rule, heuristic value, etc.). The main idea behind biased randomization is the introduction of a slight modifi- cation in the greedy constructive behavior that provides a certain degree of randomness while maintain- ing the logic behind the heuristic. BRPs can be categorized into two main groups according to how choice probabilities are computed: (i) BRPs using an empirical bias function; and (ii) BRPs using a skewed the- oretical probability distribution. This paper analyzes the second group and illustrates, throughout a series of numerical experiments, how these BRPs can benefit from parallel computing in order to significantly outperform heuristics and even simple metaheuristic approaches, thus providing reasonably good solu- tions in ‘real time’ to different problems in the areas of transportation, logistics, and scheduling.

� 2017 Elsevier Ltd. All rights reserved.

1. Introduction

A number of complex decision-making processes in real-life transportation, logistics, and production systems can be modeled as combinatorial optimization problems (Faulin, Juan, Grasman, & Fry, 2012). Among many others, some typical examples include: vehicle routing problems (VRP) (Toth & Vigo, 2014), arc routing problems (Corberán & Laporte, 2014), facility location problems (Chan, 2011), or scheduling problems (Pinedo, 2012). All these problems are NP-hard in nature, meaning that the space of poten- tial solutions grows very fast (exponential explosion) as the instance size increases. Therefore, using exact methods is not always the most efficient strategy, especially when the size of the problem instance is large and reasonably good decisions are needed in negligible computing times. Under these circumstances, heuristic-based approaches constitute an excellent alternative to exact methods (Talbi, 2009). Accordingly, a large number of heuris- tic and metaheuristic algorithms have been developed during the

last decades to solve large scale combinatorial optimization prob- lems and, eventually, support intelligent decision-making pro- cesses in a myriad of fields, including transportation, logistics, production, finance, telecommunication, Internet computing, health care, etc.

A constructive heuristic is a computational method that employs an iterative process to generate a feasible solution of rea- sonable quality. At each iteration of the solution-building process, the next ‘movement’ is selected from a list of potential candidates that has been sorted according to some criteria. Pure greedy heuris- tics always select the next ‘most promising’ movement. As a result, these heuristics are expected to generate a reasonably good solu- tion once the entire list is traversed. Notice, however, that this is a somewhat myopic behavior, since the heuristic selects the next movement without considering how the current selection will affect subsequent decisions as the list is processed downwards. Even worse, this property results in a deterministic procedure, i.e., the same solution is obtained every time the algorithm is run. Examples of such methods are the nearest neighbor for travel- ing salesman problems (Lawler, Lenstra, Rinnooy Kan, & Shmoys, 1985), the shortest processing time dispatching rule for scheduling problems (Pinedo & Chao, 1999), or the savings algorithm for VRPs (Clarke & Wright, 1964). Although these methods are easy to

A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228 217

implement and can be run almost instantaneously, the real-time solutions they provide are usually far from being optimal. To improve the quality of these heuristic solutions –and as far as more time is available–, different types of local search methods can be used to explore the solution neighborhood (Aarts & Lenstra, 1997). Typically, the neighbor selection is based on a certain logic that tries to take advantage of the specific characteristics of the optimization problem being considered. This usually leads to local optimal solutions. As in the construction phase, if the neighbor chosen is always the next ‘most promising’ movement according to some criteria, the resulting searching process will be determin- istic too.

Randomization techniques are frequently used to escape from this local optimality trap and improve the overall quality of the solution. These techniques can be incorporated either in the con- struction phase and/or the local search. Randomization allows exploring alternative solutions by selecting an element other than the ‘most promising option on the short run’. This leads to different outputs each time the entire procedure is executed. Since running a heuristic might take only a few seconds –or even less in a modern computer if the heuristic is correctly implemented and the instance size is not extremely large–, one can execute it several times, either in sequential mode or in parallel mode by using dif- ferent threads, and then select the best of the stochastic outputs. Countless metaheuristic algorithms include uniform randomiza- tion in their procedures. However, a uniform randomization of the list of candidate elements destroys the logic behind the heuris- tic greedy behavior. In order to maintain this logic, the randomiza- tion can be biased (i.e., oriented) so that higher probabilities are given to the most promising candidates. Thus, the main idea behind biased randomization is the introduction of a slight modi- fication in the greedy constructive behavior that provides a certain degree of randomness while maintaining the main logic behind the heuristic. In a seminal paper on the Monte Carlo method, King (1953) already emphasized the enormous improvement of biasing probabilities on sampling efficiency. Different methods to bias the randomization have been used in multiple contexts thereafter (Fig. 1). Among them, this paper pays special attention to the ones that use skewed (non-symmetric) theoretical probability distribu- tions in order to introduce an appropriate bias in the process of selecting elements from the list during the constructive and/or local search stages. Some skewed theoretical distributions, such as the geometric or the decreasing triangular ones, offer at least two advantages over using empirical distributions: (i) they contain at most one simple parameter, which can be easily set; and (ii)

Biased Randomized Procedures (BRPs)

Empirical Bias Functions

Skewed Probability Distribu

Fig. 1. A classification of biase

they can be sampled using well-known analytical expressions, which from a computational perspective is typically faster than other sampling techniques involving the use of loops.

In particular, the main contributions of this paper are: (i) to pro- vide a review of the most relevant biased randomized procedures (BRPs) used in the literature to solve combinatorial optimization problems; (ii) to provide a general framework for BRPs that use a skewed theoretical probability distribution to bias the selection of the next movement during the constructive and/or local search processes; and (iii) to illustrate, throughout a series of numerical experiments, how these BRPs can significantly outperform heuris- tics, and even simple metaheuristic approaches, thus providing reasonably good solutions in ‘real time’ (e.g., one or two seconds) to different transportation, logistics, and scheduling problems.

The remainder of this paper is structured as follows: Section 2 introduces the concept of randomized algorithms; Section 3 pre- sents different BRPs that use empirical bias functions; Section 4 provides a general framework for BRPs with a skewed theoretical probability distribution, and discusses the advantages of this approach over the one based on empirical bias functions; Section 5 analyzes different applications of BRPs to the fields of logistics, transportation, and scheduling; Section 6 describes a series of com- putational experiments that contribute to illustrate and quantify the potential of BRPs; finally, Section 7 summarizes the main con- tributions of the paper.

2. Randomized algorithms

There is an enormous body of literature that study probabilistic or randomized algorithms and a review of that is far beyond the scope of this paper. The reader is referred to Collet and Rennard (2006) for a review, and to Clerc (2015) for a vast discussion about the stochastic aspects of optimization. The focus of this paper is in the subset of randomized algorithms that include some type of bias in any of their random processes. A randomized algorithm uses random bits to make random choices during its execution. Unlike deterministic algorithms, different solutions are obtained every time the procedure is executed. The most successful approaches to solve large combinatorial problems take advantage of this fea- ture to perform several iterations and collect the best overall out- put. These approaches are commonly known as multi-start methods (Martí, Resende, & Ribeiro, 2013). In general terms, they all contain two differentiated phases: a construction process and a local or neighborhood search. The former diversifies the search for solu- tions while the latter intensifies this search. These two phases

tions (geometric, decreasing triangular, etc.)

Biased Random Sampling (BRS)

Parameterized BRS

Probabilistic Tabu Search

Ant Colony Systems

Heuristic Biased Stochastic Sampling (HBSS)

Value Biased Stochastic Sampling (VBSS)

d randomized procedures.

218 A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228

are repeated until a stopping criterion is satisfied. Note that the randomized procedure can be applied at either phase because there is always a discrete choice that has to be made.

Many randomized procedures found in the literature rely on uniform randomization, that is, they use the uniform probability distribution when selecting an element, neighbor, or solution. These could be categorized as uniformly-randomized algorithms. The main drawback of such approaches is that they do not benefit from the heuristic ‘common sense’: if candidate elements are ranked according to their ‘goodness’ on a given criterion, choosing one via a uniform random process fades away the advantages of the sorting. This is partially overcome by partially-randomized algo- rithms, that use uniform randomization but on a subset of candi- dates. The greedy randomized adaptive search procedure (GRASP) is the most representative and commonly used algorithm of this type. It was initially proposed by Feo and Resende (1995) and extensively used in multiple applications (Resende & Ribeiro, 2010). As a multi-start method, each GRASP iteration is composed of a construction phase and a local search. In the construction phase, all candidate elements are sorted according to a greedy evaluation function. The ‘best next’ elements, i.e., those whose addition represents the highest improvement on the objective function, are added to a restricted candidate list (RCL). The next element is chosen randomly (using a uniform distribution) from this RCL and added to the partial solution. This is performed itera- tively until there are no more candidates. A local search is then applied until a local optimal is reached. Since only the elements in the RCL can be included in the partial solution, this method can be seen as a multi-start method with a partially random con- struction heuristic. A similar idea is the window random sampling by Valls, Quintanilla, and Ballestı́n (2003), used in a resource- constrained project scheduling problem. A window parameter is defined as the maximum difference allowed between the order of the candidate and the minimal order. In the field of genetic algo- rithms, there also exists heuristics with randomization partially guided according to some criteria. This is the case, for instance, of the biased random key genetic algorithm (Gonçalves & Resende, 2011). In this genetic algorithm, the population is parti- tioned into two groups: elite and non-elite individuals. When this population is evolved to obtain the next generation, some of the children are obtained by the process of mating two randomly selected individuals, one from the elite group and one from the non-elite group. Mating is done using parameterized uniform crossover, that is, the genes from the elite parents have larger probability of being selected. This way, the randomization is par- tially biased to favor elite parent’s characteristics over the non- elite parent’s ones.

3. BRPs using empirical bias functions

A biased (oriented) randomized procedure aims at selecting the next element while capturing the best of two realms: exploitation and exploration. On the one hand, the procedure favors the most promising candidates to exploit the solution space; on the other hand, it introduces a weighted randomness degree to explore this solution space. To determine the balance between either one, a BRP may use a bias function. A bias function, biasð�Þ, is a function that assigns a non-negative weight to all elements in the candidate list. These weights are then normalized to obtain an empirical proba- bility distribution that will define the set of probabilities. Note that two extreme cases can be constructed by giving the same weight to all elements (uniform distribution), or by giving a positive weight to the top element and zero weight to the rest (greedy). All other weight allocations provide intermediate bias configurations. This section presents different BRPs that include some sort of bias func-

tion to build the empirical probability distribution. The bias func- tion is therefore an algorithm input that must be designed considering both the problem characteristics as well as the respon- siveness of the chosen heuristic –in terms of performance– to dif- ferent types of bias.

3.1. Biased random sampling

Biased random sampling (BRS) was one of the earliest BRPs employed in the literature. In these early heuristics, some specified criteria were used to bias the choice of randomly generated solu- tions. To the best of our knowledge, the first heuristics that incor- porated BRS were used to solve assembly line balance problems (Arcus, 1965; Tonge, 1965), production scheduling problems (Giffler, Thompson, & Van Ness, 1963; Heller & Logemann, 1962), location problems (Mabert & Whybark, 1977; Nugent, Vollmann, & Ruml, 1968), and an inventory management problem (Berry, Marcus, & Williams, 1977).

3.2. Parameterized biased random sampling

Parameterized BRS is a randomized method in which the prob- ability values to select the next candidates are biased according to priority rules. Numerous priority rule-based heuristics have been designed to tackle resource-constrained project scheduling prob- lems. For this vastly studied problem, Cooper (1976) presented the first BRS scheme that used nine different priority rules as weighting factor to bias the probability of choosing an activity. This probability was calculated by dividing the activity priority value by the sum of the priority values of all activities in the candidate list. Later, Drexl (1991) introduced regret based biased random sampling. This sampling technique uses the priority values indirectly via regret values. The regret of a job is the difference between its pri- ority value and the lowest overall priority value. Probabilities are then calculated using a parameter that controls the bias degree. Schirmer and Riesenberg (1997) proposed BRS variants, dubbed the normalized BRS and the modified regret based BRS, to cope with some of the drawbacks of the existing sampling approaches. The authors stated that, in general, they always outperformed uniform random sampling approaches. In a similar line of research, Valls et al. (2003) developed the b-BRS method. The b parameter establishes the probability of choosing the activity on top of the priority-rule based list. Lastly, Coelho and Tavares (2003) designed the global BRS method. Unlike previous sampling schemes, this one perturbs the priority values by adding a random value between 0 and 1. Activities are then scheduled in the order defined by the modified priority list. The reader is referred to Kolisch and Hartmann (1999) for a summary of some of these sampling tech- niques in the resource-constrained project scheduling problem.

3.3. Probabilistic Tabu search

For search heuristics, Glover (1989) introduced the first big family of metaheuristics, the probabilistic tabu search (PTS), which incorporated a BRP –usually in the move acceptance function. Tabu search (Glover, 1990) is a ‘‘higher level heuristic designed to guide other methods to escape the trap of local optimality”. PTS is an extension of tabu search that includes a skewed probabilistic ele- ment within the search. Biasing is a way to control the diversity in the search, and can be achieved by considering: (i) move attrac- tiveness (i.e., the change in the objective function); (ii) tabu status (i.e., tenure on a tabu list); and/or (iii) aspiration level (i.e., the objective function value in relation to a historical standard). The probabilistic nature of the approach can be a substitute for mem- ory when it is unclear how memory can be used to enhance the result. Some years later, Løkketangen and Glover (1996) adapted

A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228 219

PTS to zero-one mixed integer programming problems with prob- abilistic measures that were both effective and easy to implement.

3.4. Ant systems and ant colony systems

Another family of probabilistic algorithms that uses a BRP is the ant system (Dorigo, Maniezzo, & Colorni, 1996) and its subsequent variant, the ant colony system (Dorigo & Gambardella, 1997). Inspired by the behavior of real ants, these algorithms mimic the pheromone trails that insects use to establish the shortest paths. The ant system was first applied to the traveling salesman prob- lem. To complete a tour, the cities visited are chosen probabilisti- cally via Monte Carlo sampling. The probabilities in the state transition are biased using the so-called random-proportional rule to favor shorter edges with a greater amount of pheromone. The ant colony system includes three main variants, one of which is in the state transition. The modified transition follows the pseudo-random-proportional rule. This rule adds a previous step: it randomly selects a number uniformly distributed between 0 and 1; if it is below a given threshold the best edge is selected, otherwise an edge is selected according to the random- proportional rule. By calibrating the threshold, the algorithm determines the relative importance of exploitation (best next edge) versus exploration (biased random edge).

3.5. Reactive GRASP

In the previous section, GRASP was introduced as one of the most well-known partially-randomized algorithms. A variant of GRASP that includes a BRP is the so-called reactive GRASP, first proposed by Prais and Ribeiro (2000). Unlike the original GRASP, the size of the restricted candidate list in a reactive GRASP is not fixed but self-adjusted according to the quality of the solutions found during the search. A BRP is used when selecting the restric- tiveness, or size, of the candidate list (i.e., the parameter a). The algorithm starts with a discrete set of predetermined list sizes, ai. The probability of choosing a given ai from this set is drawn ini- tially from a uniform distribution. As the algorithm advances, these probabilities are biased using information collected during the search. One possible biasing strategy is to use the average values of the solutions obtained to recompute the probabilities of the dif- ferent a’s. The values that lead to better solutions will be more fre- quently used in the construction phase of the GRASP procedure.

3.6. Heuristic-biased stochastic sampling

Bresina (1996) devised a BRP called heuristic-biased stochastic sampling (HBSS) to solve scheduling problems and other con- strained optimization problems. The motivating idea again was to bias the probability of choosing the next partial solution. To avoid a complete random exploration, HBSS uses the search heuris- tic to guide this exploration. The guidance degree is determined by a given bias function. In the search, the elements of the candidate list are ranked according to the heuristic, and the bias function assigns a weight to each element. These weights are then normal- ized to be transformed into probabilities. Thus, if ri denotes the rank of the element ei, the probability of choosing ei is given by:

PðeiÞ ¼ biasðriÞ

P jbiasðrjÞ

Usually, when the heuristic accuracy is high a stronger bias (weight) is set to increase the probabilities of selecting better solu- tions. On the contrary, when the heuristic accuracy is lower a weaker bias is set to widen the exploration of the solution space.

Bresina, Drummond, Swanson, and Edgington (1994) experi- mented with the following bias functions in the telescope observa- tion scheduling problem:

� Logarithmic: biasðrÞ ¼ log�1ðr þ 1Þ � Linear: biasðrÞ ¼ 1=r � Polynomial (n = 2, 3, 4): biasðrÞ ¼ r�n � Exponential: biasðrÞ ¼ e�r � Uniform: biasðrÞ ¼ 1

The best performing bias functions for the particular problem the authors analyzed were the exponential and the second degree polynomial, which were the two functions in the middle in terms of bias strength. The HBSS approach encompasses a family of search algorithms that can be modulated via a bias function rang- ing from a greedy search to a uniform random search.

3.7. Value-biased stochastic sampling

Following a similar reasoning as in the HBSS, Cicirello and Smith (2005) proposed the value-biased stochastic sampling (VBSS) as a search heuristic. In HBSS, the bias function gives weight to the can- didates solely based on their rank, ignoring completely their heuristic values. Alternatively, VBSS not only considers rank but it also incorporates the ‘‘discriminatory power inherent in the heuristic”. Thus, according to the VBSS approach, the probability of choosing element ei is given by:

PðeiÞ ¼ biasðheuristicðeiÞÞ

P jbiasðheuristicðejÞÞ

The resulting probabilities from both BRPs (VBSS and HBSS), may differ considerably when the choices of the candidate list have very different heuristic values. Cicirello and Smith (2005) provided an illustration of such a case. The authors also tested their approach in the weighted tardiness scheduling problem with sequence-dependent setups. In their experiments, the VBSS approach was able to outperform the HBSS approach.

4. BRPs with skewed theoretical probability distributions

This section presents a general framework for BRPs that use a skewed theoretical probability distribution to bias the probabilities of the candidate elements. The BRPs considered in the previous section share a common feature: they all rely on some kind of bias function to define the choice probabilities. Using a bias function, each element in the list was assigned a different weight based on some criteria (e.g., ranking, priority rule, heuristic value), and then an empirical probability distribution was built. A random number drawn from this distribution determined the choice of the next ele- ment. Alternatively, instead of using an empirically-constructed probability distribution, one could resort directly to a theoretical probability distribution that is already skewed or non-symmetric by definition. Examples of such distributions are the geometric or the decreasing triangular. In particular, in most of our previous work we have used the geometric distribution, since it only has one parameter which determines its specific shape. Also, this parameter varies between 0 and 1, thus facilitating its setting in most practical applications. As values of this parameter get closer to 0, the more uniform-randomized the selection process will be. On the contrary, as these values get closer to 1, the more greedy the selection will be (Juan, Faulin, Ruiz, Barrios, & Caballé, 2010). This type of distributions provides a natural bias for the candidate elements in the list. For instance, Fig. 2 compares the probabilities of selecting each of the twenty-five elements of a given sorted list

Fig. 2. Use of skewed distributions to introduce randomness.

220 A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228

using: a geometric distribution with parameter 0.2 (left part), and a discretized decreasing triangular distribution (right part).

Most heuristics iteratively perform a construction phase fol- lowed by a neighborhood search. As discussed above, BRPs can be employed in any of the two stages. Regardless of the stage, there is always a discrete choice that has to be made from a list of poten- tial candidates. Potential candidates could be neighbor solutions, edges in traveling salesman problems and VRPs, jobs or machines in scheduling problems, or tasks or resources in resource- constrained project scheduling problems, for example. The list of candidates is sorted according to a problem-specific criterion, and probabilities are assigned to each element according to a skewed probability distribution. Fig. 3 displays a general pseudo- code for a BRP with a skewed distribution. The procedure requires the following inputs: (i) the list L of potential candidates, which is sorted according to the criteria provided by the heuristic; (ii) the seed s for the (uniform) pseudo-random number generator; (iii) the skewed probability distribution PD; and (iv) the parametric values p of this distribution. The procedure returns the selected element l from the sorted list.

To the best of our knowledge, the first heuristic algorithms including BRPs with a skewed theoretical probability distribution were introduced in Juan et al. (2010, 2011) in order to solve the VRP. The authors developed a hybrid algorithm that combined the classical savings heuristic (Clarke & Wright, 1964) with Monte Carlo simulation. At each step of the solution-construction process, eligible edges were sorted in a non-increasing savings list. Edges enjoyed (quasi-) exponentially diminishing probabilities that were

Procedure BRP(L, s, PD, p) 01 using seed s, generate ps 02 using , generate random 03 l select the -th element o 04 return l End

Fig. 3. Pseudocode to select the next el

variable and based upon a geometric distribution. Consequently, the next element was selected by a guided random sampling. The procedure continued until there were no more edges to be selected.

Even for large-size instances, heuristics with this type of BRPs can generate a large number of promising solutions in a few sec- onds, with some of these solutions outperforming the one pro- vided by the original heuristic (Fig. 4). Moreover, just by employing different threads the computation can be trivially parallelized by assigning a different seed to each thread, i.e., the BRP allows to generate solutions in real-time that outper- form the one provided by the original heuristic. Notice that this might be especially interesting in online optimization problems, where only one or two seconds are allowed before taking a deci- sion –which invalidates the use of time-consuming local search processes. This is the case, for instance, of the mobile cloud com- puting application analyzed in Mazza, Pagès-Bernaus, Tarchi, Juan, and Corazza (2017), where mobile users require immediate assignment but, at the same time, some intelligence must be incorporated into the assignment process to optimize the use of shared telecommunication and computing resources. Also, in De Armas, Cadarso, Juan, and Faulin (2016), the authors propose the use of BRPs for fast generation of crew rostering plans in air- line companies under realistic conditions. As discussed in Juan, Cáceres-Cruz, González-Martín, Riera, and Barrios (2014), another interesting application of these BRPs is the quick gener- ation of a diversified population of starting ‘promising’ solutions for metaheuristics.

eudo-random number in [0,1) variate from distribution PD(p) f the sorted list L

ement using a skewed distribution.

Fig. 4. Using skewed distributions to generate alternative solutions.

A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228 221

From a computational perspective, the generation of random variates coming from probability distributions other than the uni- form is always more time-consuming than the generation of pseudo-random numbers. Fortunately, there exist analytical expressions that allow fast generation of random observations from most theoretical distributions. Fig. 5 shows an example of Java code that makes use of these analytical expressions to effi- ciently generate random positions in a sorted list, with these ran- dom positions following either a geometric probability distribution or a triangular distribution. Since a large number of such random positions is usually required during the BRP execution, not using such analytical expressions might significantly increase the com- putational time required by the algorithm. This is, in our opinion, one of the main advantages of using skewed theoretical probability distributions over empirical bias functions, since the latter typi-

private static int getRandomPosGeom(doubl { // Returns random position between 0 a int index = (int) (Math.log(rng.nextDo pos = pos % n; return pos; }

private static int getRandomPosTriang(Ran { // Returns random position between 0 a pos = (int) (n * (1 - Math.sqrt(rng.ne return pos; }

Fig. 5. Efficient generation of random po

cally requires more computational time and also more parameter fine-tuning processes than the former.

5. Areas of application

This section illustrates the use of BRPs with skewed theoretical probability distributions through some examples of applications to different combinatorial optimization problems.

5.1. Vehicle routing problem

In a VRP, a set of customers’ demands must be satisfied by a fleet of capacitated vehicles that typically depart from a central depot. Moving vehicles between any two nodes (customers or depots) in the map has a distance-based cost. The goal is to find

e beta, Random rng, int n) nd n-1 based on Geometric(beta) uble()) / Math.log(1 - beta));

dom rng, int n) nd n-1 based on decreasing Triangular xtDouble())));

sitions using skewed distributions.

222 A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228

the set of vehicle routes that minimizes the delivery cost while serving all demands and taking into consideration the vehicle capacity constraints. One popular procedure for solving this prob- lem is the aforementioned savings heuristic. In the savings heuris- tic, an initial dummy solution is constructed by sending a virtual vehicle from the depot to each customer. Then, the list of edges connecting each pair of nodes is considered. This list is sorted according to the ‘savings’ criterion that would be obtained by using the corresponding edge to merge two routes in the dummy solu- tion. Thus, merging edges associated with higher savings are located at the top of the list, while edges with lower savings are located at its bottom. At this point, the sorted list of edges is tra- versed from the top to the bottom, and new route merges are car- ried out whenever the corresponding edge can be used to merge the two routes it connects without violating any constraint.

The diversification of the savings heuristic is rather old, reach- ing many variants very soon (Toth & Vigo, 2014). As far as we know, the first randomization of the savings heuristic was done by Buxey (1979), who made a random selection of one shortlist of savings according to a probability distribution built taking the savings themselves as weights. Afterwards, Fernández de Córdoba, Garcı́a Raffi, and Sanchis (1998), Fernández de Córdoba, García-Raffi, Mayado, and Sanchis (2000) developed two proce- dures using randomization to solve a real version of the Capaci- tated VRP and the Rural Postman Problem. Subsequently, the ALGACEA-1 method for the Capacitated VRP included the control of the randomization using an entropy function (Faulin & Juan, 2008).

Nevertheless, as stated above, the first implementations of a BRP with a skewed theoretical distribution was carried out by Juan et al. (2010). Later, Juan et al. (2011) improved the algorithm by incorporating some splitting and cache (memory-based) tech- niques. This conceptual idea was generalized in Juan, Faulin, Ferrer, Lourenço, and Barrios (2013) to the multi-start biased ran- domization of classical heuristics with adaptive local search (MIRHA). This solution approach was able to handle, in an efficient way, realistic vehicle routing problems under more complex sce- narios dominated by non-smooth/non-convex objective functions and non-convex regions.

BRPs of this type were also used in some VRP extensions, namely, the heterogeneous fleet VRP (Juan, Faulin, Cáceres-Cruz, Barrios, & Martinez, 2014), the heterogeneous fleet VRP with multi-trips (Cáceres-Cruz, Grasas, Ramalhinho, & Juan, 2014; Grasas, Cáceres-Cruz, Lourenço, Juan, & Roca, 2013), the VRP with asymmetric costs and heterogeneous fleets (Herrero, Rodríguez, Cáceres-Cruz, & Juan, 2014), the VRP with multiple driving ranges – i.e., heterogeneous fleet with respect to maximum route lengths – (Juan, Goentzel, & Bektas�, 2014), the multi-depot VRP with a lim- ited number of identical vehicles per depot (Juan, Pascual, Guimarans, & Barrios, 2015), and the two-dimensional loading capacitated VRP with homogeneous fleet (Dominguez, Juan, & Faulin, 2014), with heterogeneous fleet (Dominguez, Juan, Barrios, Faulin, & Agustin, 2016), with heterogeneous fleet and sequential loading and items rotation (Dominguez, Juan, De, & Ouelhadj, 2016), and with backhauls (Dominguez, Guimarans, Juan, & Nuez, 2016).

A similar ‘savings-based’ heuristic, called SHARP, was developed by González et al. (2012) for solving the arc routing problem. The arc routing problem is similar to the previously described VRP, but it differs in several details: first, demands are not located on the nodes, but on the edges connecting these nodes; second, only some nodes are directly connected among them (i.e., the underly- ing graph or network connecting the nodes in the problem is not complete). Again, the SHARP heuristic makes use of a dummy ini- tial solution and a sorted list of connecting edges to merge those routes that provide the highest possible savings at each step with-

out violating any problem constraint (e.g., vehicle capacity). The edges are selected with biased probabilities according to a geomet- ric distribution with a parameter randomly selected between 0.10 and 0.25.

5.2. Scheduling problem

Another well-known optimization problem is the permutation flow-shop problem (PFSP). This is a problem frequently encoun- tered in production processes, where a sequence of jobs or tasks has to be processed in a set of machines. Each job requires a given time to be processed by each machine, and the goal here is to find the permutation of jobs that minimizes the makespan, i.e., the total time necessary to complete the processing of all the jobs in all the machines. The NEH heuristic (Nawaz, Enscore, & Ham, 1983) is probably the best well-known heuristic for solving this problem. In the NEH heuristic, the list of jobs is sorted according to the total time each job would require to be processed by all the machines if it were the only job in the set. Then, the sorted list of jobs is tra- versed from top to bottom, and a new emerging solution (permu- tation of jobs) is constructed by locating each new job extracted from the list in the position that minimizes the makespan of the jobs considered so far. Juan, Lourenço, Mateo, Luo, and Castellà (2014) employed a BRP with a discretized version of the decreasing triangular distribution during the solution–construction process to select the jobs. This way, eligible jobs were assigned linearly diminishing probabilities according to their corresponding total processing time. In Juan, Barrios, Vallada, Riera, and Jorba (2014), the former BRP was combined with simulation in order to deal with the PFSP with stochastic processing times.

5.3. Facility location problem

The facility location problem (FLP), sometimes referred to as the location-allocation problem, consists of deciding the location of facilities and allocating demand points to one or multiple facilities (Reese, 2006). The objectives can be manifold: minimizing the cost of serving all customers (p-median problem), minimizing the long- est distance between any customer and its assigned facility (p- center problem), minimizing the sum of fixed setup costs and vari- able costs of serving the customers (uncapacitated facility location problem), minimizing a total cost that is a function of the distance and flow between the facilities plus the fixed cost of placing a facil- ity (quadratic assignment problem), among others.

Cabrera, Gonzalez-Martin, Juan, Marquès, and Grasman (2014) modeled a telecommunications problem as an uncapacitated facil- ity location problem, in which web-servers (facilities) needed to be placed in a distributed network to provide some service to a given set of customers. The authors developed a probabilistic algorithm that combined an iterated local search framework with a BRP with a geometric distribution. The use of a biased distribution led to shorter convergence times than those of a uniform distribution.

Similarly, De Armas, Juan, Marques, and Pedroso (2017) propose a new heuristic for the uncapacitated FLP, and then extend this heuristic to a BRP. They show the efficiency of this approach in solving very large-scale instances in low computing times and, then, they extend the BRP into a simheuristic (Juan, Faulin, Grasman, Rabe, & Figueira, 2015) able to deal with the stochastic version of the problem.

6. Computational experiments

In order to provide some empirical evidences on the use of BRP techniques to enhance classical heuristics and quantify the gains obtained, a series of experiments have been developed for five

A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228 223

well-known combinatorial optimization problems. The problems, heuristics, and benchmarks selected, as well as the results achieved, are described and analyzed in the next subsections.

6.1. Selected problems and heuristics

Five well-known optimization problems have been chosen to illustrate the improvements that can be reached by the introduc- tion of BRPs in constructive heuristics: the VRP, the ARP, the PFSP, the uncapacitated FLP (UFLP), and the 2D strip packing problem (2DSPP).

For the first three problems the following heuristics have been selected: the savings heuristic (Clarke & Wright, 1964) for the VRP, the SHARP heuristic (González et al., 2012) for the ARP, and the NEH heuristic (Nawaz et al., 1983) for the PFSP. These three heuristics make use of a sorted list that is traversed from the top to the bottom. Thus, at each iteration the next element in the list is chosen without knowing how this selection will condition future decisions during the solution building process. To avoid this greedy behavior, we use a skewed probability distribution to select the next element from the list.

Regarding the UFLP, this is a location problem which involves locating an undetermined number of facilities to minimize the sum of the setup costs of these facilities and the costs of serving the customers from these facilities. It is assumed that there is no limit on the number of customers that can be served from each sin- gle facility. In order to solve this problem, we have used the con- structive heuristic proposed in De Armas et al. (2017). This heuristic works as follows: for a given instance, a scenario with all facilities open is considered; then, the marginal savings or loses obtained when each facility is closed in this ‘‘all-open” scenario are computed. This way, we obtain a list of possible closures that can be sorted by the savings value. Afterwards, starting from the ‘‘all- open” scenario, the savings list is traversed from the beginning, and the next closure is performed as far as it reduces the total cost. After each closure, the savings/losses associated with closing each open facility are updated to take into account the new scenario, and the list is re-sorted accordingly. Obviously, since this heuristic also uses a dynamically-sorted saving list, a BRP technique can be introduced using a skewed probability distribution.

Finally, the 2DSPP –also referred to as the Open Dimension Problem (Wäscher, Haußner, & Schumann, 2007) – involves pack- ing items into a single bin or strip of fixed width and infinite height, with the objective of minimizing the total height of the packing within the strip. For this problem we have selected the ‘best-fit decreasing height decreasing width’ heuristic (Mumford- Valenzuela, Wang, & Vick, 2001; Ntene & van Vuuren, 2009). This heuristic is a variation of the ‘first fit decreasing height’ heuristic (Coffman, Garey, & Tarjan, 1980). Initially, all rectangles to be packed are sorted by decreasing height (or decreasing width in case of rectangles with equivalent height). Again, a BRP can be applied regarding this sorted list of rectangles to be processed.

6.2. A first experiment regarding parallelization and computing times

We have implemented in Java 8 the previously described heuristics and their corresponding multi-start biased-randomized versions for each of the five optimization problems. For all exper- iments we have used a geometric distribution with a beta param- eter adapted to each problem. A series of classical benchmarks were then run on a desktop computer (Intel Core i5 CPU @2.7 GHz with 8 GB on OS X). Each instance was run 100 times with different seeds as parallel agents, so that each agent is run- ning a certain time. Different time steps have been taken as refer- ences to compare the quality of the solutions. In order to compare the heuristic value, h, and the best value obtained with the biased-

randomized version, rh, the percentage gap between both solu- tions, computed as gap = (rh � h)/h, has been used.

More specifically, the benchmark used for the VRP was the clas- sical Kelly instances (Golden, Wasil, Kelly, & Chao, 1998). This benchmark, available at http://neo.lcc.uma.es/vrp/vrp- instances/capacitated-vrp-instances/, is composed of 20 large- scale instances, using from 200 customers to 480. Some instances have restrictions on the maximum length of every route. The benchmark used for the ARP was the classical Egl instances (Li & Eglese, 1996). This set of instances, available at http://logistik. bwl.uni-mainz.de/benchmarks.php, was constructed using as underlying graph regions of the road network of the county of Lan- cashire, UK. Costs and demands are proportional to the length of the edges, except for non-required edges that have zero demand. The Taillard benchmark is the most used benchmark in the litera- ture for the PFSP (Taillard, 1993). This set of instances, which is available at http://mistic.heig-vd.ch/taillard/problemes.dir/ordon- nancement.dir/ordonnancement.html, is composed of 120 differ- ent instances ranging from 20 jobs and 5 machines to 500 jobs and 20 machines. For the UFLP we have selected the Fpp17 bench- mark, available at http://www.math.nsc.ru/LBRT/k5/Kochetov/ bench.html. It is a set of medium-sized instances introduced by Kochetov and Ivanenko (2003). It consists of 30 instances with 307 customers and 307 facilities. Finally, the Zdf benchmark (Leung & Zhang, 2011; Zhang, Wei, Leung, & Chen, 2013) has been used for the 2DSPP. This benchmark, available at http://paginas.fe. up.pt/~esicup/datasets, is composed of large instances that were generated by combining zero-waste and non-zero waste instances.

Figs. 6–10 show the gaps for selected instances belonging to each of the benchmarks. The points identify the gaps between the heuristic solution (the highest point) and the different solu- tions obtained using the randomized version of the heuristic – without any local search added– as the number of parallel agents (executions) and the computing time increases. The lower the point the better the solution and the larger the gap with respect to the original solution provided by the heuristic. For each prob- lem, the graphs clearly show that the quality of the results increases (i.e., the negative gap increases in size) as the number of parallel executions and the total time spent for each of them increase. Therefore, in general, the most promising area of the sur- faces corresponds to the corner with the highest number of agents and highest time spent for each agent. Note that the maximum time spent is just a few seconds and there are many cases in which the improvement regarding the original heuristic is between 5% and 10%, so that a big leap in quality is obtained really fast. Being probabilistic algorithms driven by pseudo-random numbers, these biased-randomized algorithms could be easily run in parallel using multiple threads or computers, each of these employing a different seed for the generation of the random variates associated with the different skewed probability distributions. Consequently, investing the same time that the original constructive heuristic (i.e., real- time) it is possible to obtain much better solutions by simply com- bining BRPs with parallelization and multi-agent strategies, as shown in Juan, Faulin, Jorba, Caceres, and Marques (2013) and Martin et al. (2016), respectively. Of course, if more time (in the range of seconds) is permitted the quality of the solution improves even further.

6.3. A second experiment comparing BRPs with GRASP-like approaches

The previous subsection illustrates how biased-randomized and parallelized versions of different heuristics clearly outperform the heuristics themselves in a real-time optimization environment. Here, we compare the performance, also in a real-time optimiza- tion environment, of BRPs versus the traditional GRASP randomiza- tion process described in Section 2. Thus, for each of the five

Fig. 6. Gap evolution for a VRP instance (Kelly03).

Fig. 7. Gap evolution for an ARP instance (egl-s2-B).

224 A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228

instances selected (one for each optimization problem considered), both the BRP and the GRASP strategies have been executed using four different parameters (i.e., four different values of the geometric-distribution beta in the case of BRP and four different sizes of the RCL in the case of GRASP). Each of these executions consisted of ten runs (using a different seed for the random num-

ber generator at each run), allowing a maximum time of two sec- onds per run. For each problem, the best value obtained using each approach is shown in Table 1. This table also contains the original value provided by the heuristic as well as the best value obtained after running forty times a uniformly-randomized pro- cess (i.e., similar to a GRASP but without restricting the candidate

Fig. 8. Gap evolution for a PFSP instance (tai110_200_20).

Fig. 9. Gap evolution for an UFLP instance (20Fpp17).

A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228 225

list). The associated gaps of BRP, GRASP, and pure-uniform with respect to the heuristic value are also included (the more negative the gap, the higher the improvement).

The first thing to be noticed is that, even in real-time, both BRP and GRASP strategies are able to clearly improve the value pro- vided by the original heuristic approach, with negative gaps rang- ing from �1.11% in the PFSP to the �7.51% in the 2DSPP. The relatively low improvement in the case of the PFSP is probably due to the fact that the NEH heuristic used in this problem employs a simple but efficient local search mechanism. This local search might compensate from ‘bad’ decisions in the order in which jobs

are selected – from the list of potential candidates – during the constructive process. Another interesting observation is related to the extremely poor performance of the uniformly- randomization process. As expected, applying a pure-uniform selection process will completely destroy the logic behind the heuristic, thus leading to suboptimal solutions –frequently of lower quality than the one provided by the original heuristic itself – even for large computational times. Finally, observe that BRP seems to have a superior performance than GRASP, both in average values (�5.64% vs. �4.95%) as well as in the number of significant differences. As shown in Fig. 11, BRP clearly outperforms GRASP in

Fig. 10. Gap evolution for a 2DSPP instance (zdf3).

Fig. 11. Visual comparison between BRP and GRASP.

Table 1 Comparison of BRPs vs. GRASP randomization in real-time optimization.

Problem Heuristic (a) BRP (b) GRASP (c) Uniform (d) Gap (a) - (b) Gap (a) - (c) Gap (a) - (d)

VRP (kelly03) 12,594 11,718 11,860 91,181 �6.96% �5.83% 624.00% ARP (egl-S2-B) 14,124 13,476 13,692 33,118 �4.59% �3.06% 134.48% PFSP (tai110_200_20) 11,869 11,644 11,737 11,784 �1.90% �1.11% �0.72% UFLP (20FPP17) 123,245 114,327 114,330 123,349 �7.24% �7.23% 0.08% 2DSPP (zdf3) 213 197 197 219 �7.51% �7.51% 2.82% Averages – – – – �5.64% �4.95% 152.13%

226 A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228

three out-of five experiments, while showing a similar behavior in the remaining two (the more external the curve the lower the

improvement gap with respect to the solution provided by the associated heuristic).

A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228 227

7. Conclusions

This work reviews biased randomized procedures (BRPs), their different implementations, and some of their main applications in logistics, transportation, and production. The paper focuses on an emergent family of BRPs that rely on the use of skewed theoret- ical probability distributions, such as the geometric and the decreasing triangular distributions. These BRPs have two main advantages over more traditional BRPs based on empirical bias functions: (i) they are computationally faster, since they benefit from analytical expressions to generate random variates from the- oretical probability distributions; and (ii) they use at most one parameter that does not require complex and time-consuming set- ting processes.

By combining skewed probability distributions with random sampling, the logic behind the heuristic can be slightly randomized without losing its good properties. This strategy allows transform- ing deterministic heuristic procedures into probabilistic algo- rithms that can be run several times (either sequentially or in parallel) to obtain different promising solutions to the original problem, thus increasing the probability of obtaining better and diversified solutions. As the computational experiments show, the use of BRPs based on skewed probability distributions can easily and noticeably improve the performance of already existing or new heuristics.

Due to their relative simplicity, their fast execution times, and their ability to be parallelized, BRPs constitute an excellent alterna- tive to the use of simple heuristics without incurring in the compu- tational, implementation, and fine-tuning efforts required by most metaheuristics. This is especially the case in online optimization or whenever decisions must be made in real-time even for large-size instances, something that is becoming more frequent due to the growing dynamism, complexity, and responsiveness requirements of most real-life systems in areas such as logistics, transportation, production, telecommunication, finance, Internet computing, and health care.

Acknowledgements

This work has been partially supported by the Spanish Ministry of Economy and Competitiveness and FEDER (TRA2013-48180-C3- P, TRA2015-71883-REDT), the Ibero-American Program for Science and Technology for Development (CYTED2014-515RT0489), and the Erasmus+ Program (2016-1-ES01-KA108-023465).

References

Aarts, E., & Lenstra, J. K. (1997). Local search in combinatorial optimization. New York: John Wiley & Sons.

Arcus, A. L. (1965). A computer method of sequencing operations for assembly lines. International Journal of Production Research, 4(4), 259–277.

Berry, W. L., Marcus, M., & Williams, J. G. (1977). Inventory investment analysis using biased sampling techniques. Management Science, 23(12), 1295–1306.

Bresina, J. L. (1996). Heuristic-biased stochastic sampling. Proceedings of the thirteenth national conference on artificial intelligence (Vol. 1, pp. 271–278). Portland, OR: AAAI Press.

Bresina, J., Drummond, M., Swanson, K., & Edgington, W. (1994). Automated management and scheduling of remote automatic telescopes. In D. M. Pyper, & R. J. Angione (Eds). Optical Astronomy from the Earth and Moon, ASP Conference Series (Vol. 55, pp 216–233). Astronomical Society of the Pacific: San Francisco.

Buxey, G. M. (1979). The vehicle scheduling problem and Monte Carlo simulation. Journal of the Operational Research Society, 30(6), 563–573.

Cabrera, G., Gonzalez-Martin, S., Juan, A. A., Marquès, J. M., & Grasman, S.E. (2014). Combining biased random sampling with metaheuristics for the facility location problem in distributed computer systems. In A. Tolk, S. Y. Diallo, I. O. Ryzhov, L. Yilmaz, S. Buckley, & J. A. Miller (Eds). Proceedings of the 2014 Winter Simulation Conference (pp 3000–3011). IEEE Press: Savannah, GA.

Cáceres-Cruz, J., Grasas, A., Ramalhinho, H., & Juan, A. A. (2014). A savings-based randomized heuristic for the heterogeneous fixed fleet vehicle routing problem with multi-trips. Journal of Applied Operational Research, 6(2), 69–81.

Chan, Y. (2011). Location theory and decision analysis: Analytics of spatial information technology (2nd ed.). Berlin, Heidelberg: Springer.

Cicirello, V. A., & Smith, S. F. (2005). Enhancing stochastic search performance by value-biased randomization of heuristics. Journal of Heuristics, 11(1), 5–34.

Clarke, G., & Wright, J. (1964). Scheduling of vehicles from a central depot to a number of delivering points. Operations Research, 12(4), 568–581.

Clerc, M. (2015). Guided randomness in optimization (Vol. 1) London: Wiley-ISTE. Coelho, J., & Tavares, L. (2003). Comparative analysis of metaheuristics for the resource

constrained project scheduling problem. Technical report. Portugal: Department of Civil Engineering, Instituto Superior Tecnico.

Coffman, E. G., Garey, D. S., & Tarjan, R. E. (1980). Performance bounds for level oriented two-dimensional packing algorithms. SIAM Journal on Computing, 9(4), 808–826.

Collet, P., Rennard, J. P. (2006). Stochastic optimization algorithms. In J. P. Rennard (Ed). Handbook of research on nature inspired computing for economics and management (pp 28–44). Idea Group Inc.: Hershey, PA.

Cooper, D. F. (1976). Heuristics for scheduling resource-constrained projects: An experimental investigation. Management Science, 22(11), 1186–1194.

Corberán, A., & Laporte, G. (2014). Arc routing: Problems, methods, and applications. Philadelphia, PA: SIAM.

De Armas, J., Cadarso, L., Juan, A. A., & Faulin, J. (2016). A multi-start randomized heuristic for real-life crew rostering problems in airlines with work-balancing goals. Annals of Operations Research. http://dx.doi.org/10.1007/s10479-016- 2260-y.

De Armas, J., Juan, A. A., Marques, J. M., & Pedroso, J. (2017). Solving the deterministic and stochastic uncapacitated facility location problem: From a heuristic to a simheuristic. Journal of the Operational Research Society. http://dx. doi.org/10.1057/s41274-016-0155-6.

Dominguez, O., Guimarans, D., Juan, A. A., & Nuez, I. (2016). A biased-randomised large neighbourhood search for the two-dimensional vehicle routing problem with backhauls. European Journal of Operational Research. http://dx.doi.org/ 10.1016/j.ejor.2016.05.002.

Dominguez, O., Juan, A. A., Barrios, B., Faulin, J., & Agustin, A. (2016). Using biased randomization for solving the two-dimensional loading vehicle routing problem with heterogeneous fleet. Annals of Operations Research, 236(2), 383–404.

Dominguez, O., Juan, A. A., De, Nuez I., & Ouelhadj, D. (2016). An ILS-biased randomization algorithm for the two-dimensional loading HFVRP with sequential loading and items rotation. Journal of the Operational Research Society, 67(1), 37–53.

Dominguez, O., Juan, A. A., & Faulin, J. (2014). A biased-randomized algorithm for the two-dimensional vehicle routing problem with and without item rotations. International Transactions in Operational Research, 21(3), 375–398.

Dorigo, M., & Gambardella, L. M. (1997). Ant colony system: A cooperative learning approach to the traveling salesman problem. IEEE Transactions on Evolutionary Computation, 1(1), 53–66.

Dorigo, M., Maniezzo, V., & Colorni, A. (1996). The ant system: Optimization by a colony of cooperating agents. IEEE Transactions on Systems, Man, and Cybernetics-Part B, 26(1), 29–41.

Drexl, A. (1991). Scheduling of project networks by job assignment. Management Science, 37(12), 1590–1602.

Faulin, J., & Juan, A. A. (2008). The ALGACEA-1 method for the capacitated vehicle routing problem. International Transactions in Operational Research, 15(5), 599–621.

Faulin, J., Juan, A. A., Grasman, S. E., & Fry, M. J. (2012). Decision making in service industries: A practical approach. Boca Raton, FL: CRC Press.

Feo, T., & Resende, M. (1995). Greedy randomized adaptive search procedures. Journal of Global Optimization, 6(2), 109–133.

Fernández de Córdoba, P., Garcı́a Raffi, L. M., & Sanchis, J. M. (1998). A heuristic algorithm based on Monte Carlo methods for the Rural Postman Problem. Computers & Operations Research, 25(12), 1097–1106.

Fernández de Córdoba, P., García-Raffi, L. M., Mayado, A., & Sanchis, J. M. (2000). A real delivery problem dealt with Monte Carlo Techniques. Top, 8(1), 57–71.

Giffler, B., Thompson, G. L., & Van Ness, V. (1963). Numerical experience with the linear and Monte Carlo algorithm for solving production scheduling problems. In J. F. Muth & G. L. Thompson (Eds.), Industrial scheduling. Prentice-Hall Inc: Englewood Cliffs, NJ.

Glover, F. (1989). Tabu Search - Part I. ORSA Journal on Computing, 1(3), 190–206. Glover, F. (1990). Tabu search: A Tutorial. Interfaces, 20(4), 74–94. Golden, B., Wasil, E., Kelly, J., & Chao, I. (1998). The impact of metaheuristics on

solving the vehicle routing problem: Algorithms, problem sets, and computational results. In T. G. Crainic & G. Laporte (Eds.), Fleet management and logistics (pp. 33–56). New York: Springer.

Gonçalves, J. F., & Resende, M. G. C. (2011). Biased random-key genetic algorithms for combinatorial optimization. Journal of Heuristics, 17(5), 487–525.

González, S., Juan, A. A., Riera, D., Castellà, Q., Muñoz, R., & Pérez, A. (2012). Development and assessment of the SHARP and RandSHARP algorithms for the Arc Routing Problem. AI Communications, 25(2), 173–189.

Grasas, A., Cáceres-Cruz, J., Lourenço, H. R., Juan, A. A., & Roca, M. (2013). Vehicle routing in a Spanish distribution company: Saving using a savings-based heuristic. OR Insight, 26(3), 191–202.

Heller, J., & Logemann, G. (1962). An algorithm for the construction and evaluation of feasible schedules. Management Science, 8(2), 168–183.

Herrero, R., Rodríguez, A., Cáceres-Cruz, J., & Juan, A. A. (2014). Solving vehicle routing problems with asymmetric costs and heterogeneous fleets. International Journal of Advanced Operations Management, 6(1), 58–80.

228 A. Grasas et al. / Computers & Industrial Engineering 110 (2017) 216–228

Juan, A. A., Barrios, B., Vallada, E., Riera, D., & Jorba, J. (2014). SIM-ESP: A simheuristic algorithm for solving the permutation flow-shop problem with stochastic processing times. Simulation Modelling Practice and Theory, 46, 101–117.

Juan, A. A., Cáceres-Cruz, J., González-Martín, S., Riera, D., & Barrios, B. B. (2014). Biased randomization of classical heuristics. In J. Wang (Ed), Encyclopedia of Business Analytics and Optimization (Vol. 1, pp. 314–324). IGI Global Books: Hershey, PA.

Juan, A. A., Faulin, J., Cáceres-Cruz, J., Barrios, B. B., & Martinez, E. (2014). A successive approximations method for the heterogeneous vehicle routing problem: Analysing different fleet configurations. European J. Industrial Engineering, 8(6), 762–788.

Juan, A. A., Faulin, J., Ferrer, A., Lourenço, H. R., & Barrios, B. (2013). MIRHA: Multi- start biased randomization of heuristics with adaptive local search for solving non-smooth routing problems. Top, 21(1), 109–132.

Juan, A. A., Faulin, J., Grasman, S., Rabe, M., & Figueira, G. (2015). A review of Simheuristics: Extending metaheuristics to deal with stochastic optimization problems. Operations Research Perspectives, 2, 62–72.

Juan, A. A., Faulin, J., Jorba, J., Caceres, J., & Marques, J. (2013). Using parallel & distributed computing for solving real-time vehicle routing problems with stochastic demands. Annals of Operations Research, 207, 43–65.

Juan, A. A., Faulin, J., Jorba, J., Riera, D., Masip, D., & Barrios, B. (2011). On the use of Monte Carlo simulation, cache and splitting techniques to improve the Clarke and Wright savings heuristics. Journal of the Operational Research Society, 62(6), 1085–1097.

Juan, A. A., Faulin, J., Ruiz, R., Barrios, B., & Caballé, S. (2010). The SR-GCWS hybrid algorithm for solving the capacitated vehicle routing problem. Applied Soft Computing, 10(1), 215–224.

Juan, A. A., Goentzel, J., & Bektas�, T. (2014). Routing fleets with multiple driving ranges: Is it possible to use greener fleet configurations? Applied Soft Computing, 21, 84–94.

Juan, A. A., Lourenço, H. R., Mateo, M., Luo, R., & Castellà, Q. (2014). Using iterated local search for solving the flow-shop problem: Parallelization, parametrization, and randomization issues. International Transactions in Operational Research, 21 (1), 103–126.

Juan, A. A., Pascual, I., Guimarans, D., & Barrios, B. B. (2015). Combining biased randomization with iterated local search for solving the multidepot vehicle routing problem. International Transactions in Operational Research, 22(4), 647–667.

King, G. W. (1953). The Monte Carlo Method as a natural mode of expression in operations research. Operations Research, 1(2), 46–51.

Kochetov, Y. & Ivanenko, D. (2003). Computationally difficult instances for the uncapacitated facility location problem. In Proceedings of the 5th Metaheuristics International Conference (MIC) (pp. 41:1–41:6).

Kolisch, R. & Hartmann, S. (1999). Heuristic algorithms for the resource-constrained project scheduling problem: Classification and computational analysis. In J, Weglarz (Ed). Project scheduling: Recent models, algorithms and applications (pp 147–178). Kluwer Academic Publishers: Dordrecht, The Netherlands.

Lawler, E. L., Lenstra, J. K., Rinnooy Kan, A. H. G., & Shmoys, D. B. (1985). The traveling salesman problem. Chichester: Wiley & Sons.

Leung, S. C. H., & Zhang, D. (2011). A fast layer-based heuristic for non-guillotine strip packing. Expert Systems with Applications, 38(10), 13032–13042.

Li, L. Y. O., & Eglese, R. W. (1996). An interactive algorithm for vehicle routeing for winter-gritting. Journal of the Operational Research Society, 47(2), 217–228.

Løkketangen, A. & Glover, F. (1996). Probabilistic move selection in Tabu search for zero-one mixed integer programming problems. In I. H. Osman & J. P. Kelly

(Eds). Meta-heuristics: Theory & applications (pp 467–487). Kluwer Academic Publishers: Boston, MA.

Mabert, V. A., & Whybark, D. C. (1977). Sampling as a solution methodology. Decision Sciences, 8(1), 167–179.

Martí, R., Resende, M.G.C., & Ribeiro,C. C. (2013). Multi-startmethods for combinatorial optimization. European Journal of Operational Research, 226(1), 1–8.

Martin, S., Ouelhadj, D., Beullens, P., Ozcan, E., Juan, A. A., & Burke, E. (2016). A multi-agent based cooperative approach to scheduling and routing. European Journal of Operational Research, 254(1), 169–178.

Mazza, D., Pagès-Bernaus, A., Tarchi, D., Juan, A. A., & Corazza, G. E. (2017). Supporting mobile cloud computing in smart cities via randomized algorithms. IEEE Systems Journal. http://dx.doi.org/10.1109/JSYST.2016.2578358.

Mumford-Valenzuela, C., Wang, P. Y., & Vick, J. (2001). Heuristic for large strip packing problems with guillotine patterns: An empirical study. In Proc. 4th Metaheuristics Int. Conference (pp. 417–421).

Nawaz, M., Enscore, E. E., & Ham, I. (1983). A heuristic algorithm for the m-machine, n-job flow-shop sequencing problem. Omega, 11(1), 91–95.

Ntene, N., & van Vuuren, J. H. (2009). A survey and comparison of guillotine heuristics for the 2D oriented offline strip packing problem. Discrete Optimization, 6(2), 174–188.

Nugent, C. E., Vollmann, T. E., & Ruml, J. (1968). An experimental comparison of techniques for the assignment of facilities to locations. Operations Research, 16 (1), 150–173.

Pinedo, M. L. (2012). Scheduling: Theory, algorithms, and systems (4th ed.). Springer Science & Business Media: New Jersey.

Pinedo, M., & Chao, X. (1999). Operations scheduling with applications in manufacturing and services. Boston, MA: Irwin/McGraw-Hill.

Prais, M., & Ribeiro, C. C. (2000). Reactive GRASP: An application to a matrix decomposition problem in TDMA traffic assignment. INFORMS Journal on Computing, 12(3), 164–176.

Reese, J. (2006). Solution methods for the p-median problem: An annotated bibliography. Networks, 48(3), 125–142.

Resende, M. G. C., & Ribeiro, C. C. (2010). Greedy randomized adaptive search procedures: Advances, hybridizations, and applications. In M. Gendreau & J. Y. Potvin (Eds.), Handbook of metaheuristics (pp. 283–319). New York: Springer, US.

Schirmer, A. & Riesenberg, S. (1997). Parameterized Heuristics for Project Scheduling - Biased Random Sampling Methods. Technical Report 456, Manuskripte aus den Instituten für Betriebswirtschaftslehre der Universität Kiel.

Taillard, E. (1993). Benchmarks for basic scheduling problems. European Journal of Operational Research, 64(2), 278–285.

Talbi, E.-G. (2009). Metaheuristics: From design to implementation. New Jersey: Wiley Publishing.

Tonge, F. M. (1965). Assembly line balancing using probabilistic combinations of heuristics. Management Science, 11(7), 727–735.

Toth, P., & Vigo, D. (2014). Vehicle routing: Problems, methods, and applications (2nd ed.). Philadelphia, PA: SIAM.

Valls, V., Quintanilla, S., & Ballestı́n, F. (2003). Resource-constrained project scheduling: A critical activity reordering heuristic. European Journal of Operational Research, 149(2), 282–301.

Wäscher, G., Haußner, H., & Schumann, H. (2007). An improved typology of cutting and packing problems. European Journal of Operational Research, 183(3), 1109–1130.

Zhang, D., Wei, L., Leung, S. C. H., & Chen, Q. (2013). A binary search algorithm base on randomized local search for the rectangular strip packing problem. INFORMS Journal on Computing, 25(2), 332–345.