## Abstract

This article presents an Evolution Strategy (ES)--based algorithm, designed to self-adapt its mutation operators, guiding the search into the solution space using a Self-Adaptive Reduced Variable Neighborhood Search procedure. In view of the specific local search operators for each individual, the proposed population-based approach also fits into the context of the Memetic Algorithms. The proposed variant uses the Greedy Randomized Adaptive Search Procedure with different greedy parameters for generating its initial population, providing an interesting exploration–exploitation balance. To validate the proposal, this framework is applied to solve three different -Hard combinatorial optimization problems: an Open-Pit-Mining Operational Planning Problem with dynamic allocation of trucks, an Unrelated Parallel Machine Scheduling Problem with Setup Times, and the calibration of a hybrid fuzzy model for Short-Term Load Forecasting. Computational results point out the convergence of the proposed model and highlight its ability in combining the application of move operations from distinct neighborhood structures along the optimization. The results gathered and reported in this article represent a collective evidence of the performance of the method in challenging combinatorial optimization problems from different application domains. The proposed evolution strategy demonstrates an ability of adapting the strength of the mutation disturbance during the generations of its evolution process. The effectiveness of the proposal motivates the application of this novel evolutionary framework for solving other combinatorial optimization problems.

## 1 Introduction

In this article, the class of Evolutionary Algorithms (EA) known as Evolution Strategies (ES) (Beyer and Schwefel, 2002) is investigated and a new ES variant is proposed, being able to solve challenging combinatorial optimization problems. We combine the diversity of the population-based algorithms with the power of Reduced Variable Neighborhood Search (RVNS) (Hansen and Mladenović, 2001), a well-known trajectory search algorithm. Moreover, we supply RVNS with adaptive rules, producing a new variant that we call Adaptive RVNS (ARVNS). EA are search and optimization methods inspired by well-known evolutionary principles, such as random mutations and selective pressure for evolution and adaptation of its population. On the other hand, VNS and other trajectory-based approaches such as Tabu Search (TS) (Glover, 1996) and Iterated Local Search (ILS) (Lourenço et al., 2003) take advantage of the flexibility in designing and exploring different Neighborhood Structures (NS) of the problem, using the simple fact that the local minimum with respect to one neighborhood structure is not necessarily so with respect to another.

Biological evolution in nature is an inspiration for the ES, which is mainly guided by operators: mutation and selection. This class of methods has often been applied for solving continuous optimization problems (Kashan et al., 2015; Chaquet and Carmona, 2012; Andersen and Santos, 2012; Aler et al., 2012; Costa and Oliveira, 2001). While EA have already been applied to solve several combinatorial optimization problems (Prado et al., 2014; Qaurooni and Akbarzadeh-T, 2013; Freitas and Guimarães, 2011), only a few articles in the literature address combinatorial optimization problems using ES (Cai and Thierauf, 1996; Rajasekaran, 2006; Kashan et al., 2015). Our current article focuses on a mechanism for guiding the search in the solutions space, mainly based on a weighting system for applying move operations from *k _{max}* distinct neighborhoods. The motivation to develop an ES for combinatorial optimization, combined with the use of different NS, comes from the successful applications in both fields, numerical and combinatorial optimization.

Adaptive local search techniques have been exploited by researchers (Dong et al., 2015; Li et al., 2015; Schneider et al., 2014; Hosny and Mumford, 2010). This family of methods has the ability of exploring attraction basins with iterative moves, combining it with smart strategies that, in general, check the success of previous steps done by the methods. In this context, an ILS with self-adaptive shaking procedure was applied for tackling a flow shop problem (Dong et al., 2015) and multi-depot Vehicle Routing Problem (VRP) (Li et al., 2015). Following the same idea, Schneider et al. (2014) proposed an adaptive mechanism for guiding the shaking step of a VNS applied on a VRP. Their approach selected and favored route and vertex according to their success within the search.

By Reduced VNS method a random point from the *k*-th neighborhood () of the current incumbent solution *x* is taken, and no descent from there is made. RVNS has been shown to be useful in solving large instances, for which local search is costly (Xiao et al., 2011; Hansen et al., 2009).

The adaptive variant of the RVNS designed in this article, the ARVNS, explores specific parts of each *N _{k}*, playing with probabilities evolved through the ES evolutionary process. Our proposal implicitly considers the problem-specific characteristics and the success of a given

*N*within the search. There is no need of mechanisms for analyzing previous success of the

_{k}*N*, since it is inherited through the genes, ES mutation operators, from the parents to the offspring. The latter are generated after mutations over the ARVNS application probabilities and are expected to survive if interesting changes have occurred.

_{k}On the other hand, the proposed self-adaptive evolution strategy also fits the context of Memetic Algorithms (MA) (Moscato and Cotta, 2003), comprising a population-based method in which individuals possess specific, unique, or special searching operators. By incorporating problem domain knowledge, such as heuristics, local search techniques, individual learning strategies, or other stochastic search operators, MA evolve their population through the generations. These techniques have shown a great potential for solving -Hard combinatorial and continuous optimization problems (Guimaraes et al., 2007; Wanner et al., 2008; de Freitas et al., 2014).

We introduce an evolution strategy--based algorithm, abbreviated as GES, which generates its initial population through a diversified and greedy procedure. Thus, we suggest the use of the Greedy Randomized Adaptive Search Procedures (GRASP) (Resende and Ribeiro, 2010). GRASP construction phase is responsible for filling the initial population with individuals generated with different random greedy parameters. Other versions with simpler solution generation procedures could also be used, such as generating solutions at random or initializing a homogeneous initial population. We kept the ES mutation operator vectors as the main search operator in the solution space, however, being guided by new data structures incorporated with each individual representation. The proposed adaptive operators regulate the rate of moves application of different NS, in a special case of an RVNS search. An optional intensification phase is also suggested for some problems, in order to accelerate the convergence of the algorithm.

Problems already tackled by the authors of this article were used as case studies. We took profit from the optimization framework OptFrame (see Coelho, Munhoz et al., 2011), a computational framework for the development of efficient metaheuristic algorithms for combinatorial optimization problems. Three different (and very challenging) combinatorial optimization problems are considered in this work as case studies:

An Open-Pit-Mining Operational Planning (OPMOP) problem (Souza et al., 2010);

An Unrelated Parallel Machine Scheduling Problem with Setup Times (UPMSP-ST) (Al-Salem, 2004), implemented in a Java platform with OptFrame ideas;

Calibration of a hybrid fuzzy model calibration for the Short-Term Load Forecasting Problem (STLFP) (Coelho et al., 2014a).

Since we are dealing with many -Hard problems, exact solution methods have restricted applicability. This fact motivates us to search for solutions by means of metaheuristic procedures.

The abstraction of the concepts involving the proposed discrete ES and its application over the aforementioned problems was not found in the literature and shows up as a novel evolutionary framework for combinatorial optimization problems. The major contributions of the current work are to:

Use and adapt a population-based algorithm for combinatorial optimization problems, combining it with trajectory search techniques. ---Walk through the search space by using discrete moves, following atrajectory provided by random moves in a Reduced Variable NeighborhoodSearch;---Combine neighborhood structures in order to guide the search for newsolutions using self-adaptive mutation operators;

Introduce a flexible self-adaptive search framework implemented in the core of the open-source optimization framework OptFrame;

Apply the proposed methodology for solving different -Hard combinatorial optimization problems, including two large-scale real case problems (OPMOP and HFVRPMT).

The remainder of this article is organized as follows. Section 2 details the proposed self-adaptive ES. Section 3 describes an application for the OPMOP. Analogously, Sections 4 and 5 describe mutation operators’ behavior of the proposed algorithm over the STLFP and UPMSP-ST, respectively. Section 6 draws the final considerations and future works.

## 2 Self-Adaptive Evolution Strategy

Evolution Strategies (ES) were developed in the 1960s and 1970s by P. Bienert, I. Rechenberg, and H.-P Schwefell at the Technical University of Berlin. The initial version operated with single individuals, subjected to mutation and selection among their descendants. Beyer and Schwefel (2002) provide a detailed description about ES in their comprehensive introduction.

As mentioned by Meyer-Nieberg and Beyer (2007), self-adaptation can be seen as state-of-the-art methods for adjusting the setting of control parameters. A study about self-adaptation in EA applied for Combinatorial Optimization can be found in Smith (2008).

ES uses natural problem-dependent representations according to each problem that is being tackled. One advantage is its searching ability over the evolution process that is guided, primarily, by mutation and selection. Here, we take advantage of the basic principle, introducing the ARVNS, an RVNS guided by probabilities. A wide range of genetic operators can be used in order to generate the offspring population. Novel mechanisms are still being explored and developed in the literature (Chuang et al., 2015). Mutation operators from ES usually can change all components of a parent vector at the same time, but with minor changes since it is assumed that in the real biological evolution small mutations occur frequently but large ones only rarely.

In the 1990s, Cai and Thierauf (1996) proposed a general ES for solving discrete optimization problems, suggesting that not all components of a parent vector should be mutated, but only a few should be randomly changed every time. This strategy was quite smart, since, in discrete sets, differences between any two adjacent values are usually not small. This approach has been used/followed in different applications (Hasancebi, 2007; Chen and Chen, 2009; Yao et al., 2011). Li et al. (2013) introduced a Mixed Integer ES capable of handling parameters consisting of discrete and integer variables.

Motivated by successful applications in the literature, it was felt that the ES could deal with combinatorial optimization problems that Coelho, Souza et al. (2011) in partnership with a self-adaptive strategy based on moves generated from simple NSs, where neighborhood structures guide the individuals through the solution space. Thus, a special care in the design of the algorithm was given, exploring the RVNS ability of handling random moves in *N _{k}*, . The details will be discussed later in this article. The use of NS is well known in the literature and has been widely applied for solving -Hard problems (Johnson et al., 1988; Kirkpatrick, 1984; Glover, 1989; Mladenović and Hansen, 1997; Lourenço et al., 2003; Lust and Teghem, 2010; Pisinger and Ropke, 2010).

In this sense, we provide a compact and efficient encoding for adapting NS use and strength in connection with individual mutation operators, described in the next section.

### 2.1 Mutation Operators

According to the problem that is being solved, a desired data structure, or *representation*, is selected and a solution to the problem is defined as *s*. Here, each individual is comprised of two additional mutation vectors, defined as *P* and *A*, presented in Equations (1) and (2) respectively, embedded within the solution representation *s*.

*s*,

*P*, and

*A*.

The first mutation vector, *P*, represents the likelihood associated with the choice of each NS in a set composed of neighborhood structures. This vector guides the probability of applying each NS used to walk on the search space. Each position stores the probability of the application of a given move .

The second mutation vector, *A*, stores integer values for controlling the strength of the disturbance, once an NS is selected to be applied and shake the solution *s*. Each position of this vector indicates the number of random moves to be applied from neighborhood *k*, with *nap _{k}* representing the maximum

*number of applications*of different moves in each mutation event.

Both vectors are adapted across the generations of the evolutionary process, according to well-known probability distribution functions. The evolution of these two mutation vectors will be discussed for each application described in this current study. Other sets of parameters could also be included for adapting the distributions along the generations.

### 2.2 Generic Evolution Strategy Pseudocode Guideline

The proposed self-adaptive evolution strategy algorithm pseudocode is outlined in Algorithm 1. As emphasized by Lust and Teghem (2010), generating an initial population diversified and with good potential is a very important feature for the convergence of population-based algorithms. Thus, we suggest the GRASP procedure in partnership with the proposed algorithm.

The initial population (lines 1 to 7 of Algorithm 1) consists of generating a set of individuals. Line 3 calls the GRASP procedure and generates each solution of the initial population with different random greedy parameter . Achieving a diversified initial population is an important stage for the algorithm convergence, as can be verified in Lust et al. (2011); therefore, different parameters are used in order to control the size of the candidates' list. Thus, a random GRASP is designed here. In the second step (line 4 of Algorithm 1), the self-adaptive mutation vectors *P* and *A* are built for each individual. The procedure BuildMutationVectors (outlined in Algorithm 2) describes a generic and simple idea for generating initial values for the mutation vectors. Line 5 merges the triple formed by a GRASP solution *s* and the mutation operators *P* and *A*. In this sense, the following nomenclature is defined: let *ind ^{S}* be the solution

*s*of the individual ;

*ind*be the probability parameter vector;

^{P}*ind*be the application parameter vector.

^{A}Algorithm 2 fills the probabilities vector, *P*, random numbers generated between the interval and the same idea is applied for the vector of applications, *A*, respecting the range .

In line 11 of Algorithm 1, individual parameters are updated; the pseudocode of the procedure “UpdateParameters” is described in Algorithm 3. Vectors of mutation parameters *A* and *P* are updated and adapted according to a Normal or Binomial Distribution, both centered at mean zero and standard deviation and , respectively. For the binomial distribution, an additional parameter indicating the number of trials is required. Parameter regulates the probability of successes in a sequence of independent yes/no experiments. Updates of the vectors *ind ^{P}* and

*ind*can be viewed in lines 2 and 3, respectively, in Algorithm 3. Line 5 checks if the limits of both mutation operators are respected after the mutation. As expected, the maximum value assigned for the application probability, of each cell from vector

^{A}*ind*, should be between 0% to 100%. Following a similar reasoning, the maximum number of applications, for each cell

^{P}*k*from vector

*ind*, should be

^{A}*nap*and no less than 1.

_{k}Line 12 of Algorithm 1 calls the mutation procedure, a special case of the classical RVNS, called ARVNS, illustrated in Algorithm 4. In line 2 of Algorithm 4, a random number is generated and then line 3 checks if this number satisfies the probability . In the positive case, the neighborhood structure *N _{k}*, allocated in the current index

*k*, is applied times. The neighborhood order of this parameter vector is chosen at random. An optional mutation rate parameter can be added in each individual representation in order to regulate the sequence in which the NS are applied.

Line 17 of Algorithm 1 opens the possibility of calling an optional intensification phase, usually done by local search techniques, such as Variable Neighborhood Descent (VND). When triggered, an intensification phase refines random solutions in the offspring population. We emphasize this phase as optional since several variations could be implemented here according to the combinatorial optimization problem that is being tackled.

Finally, the selection procedure, line 19 of Algorithm 1, can be any desired selection strategy, as long as it returns a population with cardinality . We used two basic forms of competition, both with the same notation of Beyer and Schwefel (2002). In the first one, denoted by (), there is competition between parents and offspring. In this strategy, best individuals are selected among the union of parents and offspring. In the second selection strategy, denoted by (), only the offspring compete for survival. It is clear that using the strategy as a way of selection, the population of the next generation suffers a considerably higher selective pressure than using the strategy ().

## 3 Open-Pit-Mining Operational Planning Problem

The OPMOP involves the allocation of mining equipment to the pits, which may be of ore or waste rocks, as well as determining the number of trips for each truck so that both the production goals and the desired mineral composition of the ore are fulfilled. The goal is to find a mining rate on every pit that minimizes deviations from production goals, quality, and also the number of trucks required for the process. Dynamic truck allocation is considered, which is the possibility to allocate the trips from a certain truck to a different pit. This allocation system contributes to an increased fleet productivity and, therefore, to reduce the number of trucks needed for the production process. Figure 1 shows a graphical example of the OPMOP, composed of pits with different mineral compositions, two shovels, and different trucks.

A brief literature review is described next.

Costa (2005) developed a heuristic algorithm based on GRASP and VNS using six different types of movements to explore the solution space. A comparison was made between the results obtained by this heuristic algorithm and those found by the solver LINGO, version 7, applied to a mathematical programming model developed in Costa et al. (2004). Results showed that the heuristic algorithm was able to find better solutions faster. Guimarães et al. (2007) presented a computer simulation model to validate results obtained by applying a mathematical programming model to determine the mining rate in open-pit mines (e.g., occurrence of queues).

Souza et al. (2010) proposed an algorithm, called GGVNS, which combines the metaheuristics General Variable Neighborhood Search (GVNS) (Hansen et al., 2008) and GRASP procedure. The GVNS was chosen due to its simplicity, efficiency, and capacity of its natural local search to deal with different neighborhoods. The authors compared the results generated by GGVNS with those achieved by CPLEX optimizer 11.0.1, using eight test problems. Computational experiments showed that the algorithm was competitive and, in most instances, capable of finding new optimal solutions—with a gap —requiring a short computational time.

Coelho et al. (2012) developed the first multi-objective application to the OPMOP. Three multi-objective heuristic algorithms were validated based on Two-phase Pareto Local Search with VNS (2PPLS-VNS), proposed by Lust and Teghem (2010), Multi-objective Variable Neighborhood Search (MOVNS), presented by Geiger (2004), and Non-dominated Sorting Genetic Algorithm II (NSGA-II) developed by Deb et al. (2002). Approximations of Pareto sets generated by the developed algorithms were compared considering the hypervolume and spacing metrics. Computational experiments have shown the superiority of the algorithms based on VNS methods, which were able to find better sets of nondominated solutions, more diversified and with an improved convergence.

### 3.1 Representation and Evaluation of a Solution

Given a set of mining pits *F*, a set of trucks *T*, and a set of shovels *K*, a solution for OPMOP is represented by a matrix , where *Y* is a matrix and *N* a matrix . Each cell *y _{i}* of the matrix represents shovel allocated to the pit . If there aren’t trips made to the pit

*i*, the shovel

*k*associated with this pit is considered inactive and it is not penalized for a production below the minimum limit.

In the matrix , each cell *n _{il}* represents the number of trips performed by the truck to the pit . The value 0 (zero) means no trip to that truck. The value

*X*means that the truck is incompatible with the shovel allocated to the pit.

In Table 1, there is an example of a possible solution to the OPMOP. At the column , line *F*_{1}, the pair , indicates that the loading equipment is allocated to the pit *F*_{1} and the number 1 means that it is operating. At the column , line *F*_{3}, the pair indicates that the loading equipment is allocated to the pit *F*_{3}, but it is not operative. Finally, in line *F*_{2}, the value means that there is no loading equipment allocated to the pit *F*_{2} and, therefore, this pit is available. The other columns represent the number of trips from the truck to the corresponding pit, considering the compatibility between the truck and loading equipment allocated to the front. Cells with values *X* indicate incompatibility between a truck and its loading equipment.

### 3.2 Mathematical Model and Solution Evaluation

Equation (4) seeks to minimize the positive and negative deviations from the goals of each control parameter *j* of the mixture, and , respectively, as well as the positive and negative deviations from the production goals of ore and waste rocks, represented, in this order, by decision variables , , , and . This function also considers the minimization of the number of used trucks, represented by the binary variable *U _{l}*, which is 1 if the truck

*l*is used and 0, otherwise.

The constants , , , , , , and are weights that reflect the importance of each component of the objective function.

*f*to be minimized, composed of two parts. The first one is the actual objective function,

*f*, given by Equation (4), and the second one consists of functions that penalize the occurrence of infeasibility in the solution. Thus, the function

^{PM}*f*, given by Equation (5), measures the deviation of the goals and penalizes any violation of the constraints in the problem. where:

*f ^{PM}*(

*s*) is a function that evaluates

*s*with regard to the production goals and the quality of the final mixture;

*f ^{p}*(

*s*) penalizes

*s*if the limits for the production of ore and waste rocks are not respected;

(*s*) penalizes *s* if the limits for the *j*-th control parameter of the mixture are not respected;

(*s*) penalizes *s* if the maximum utilization rate for the *l*-th truck is exceeded;

(*s*), penalizes *s* if the productivity limits for the shovel *k* are not respected.

### 3.3 Neighborhood Structures

In order to explore the solution space, as described in Section 2, eight neighborhood structures, introduced by Souza et al. (2010), are used to analyze the convergence of the proposed GES. As detailed in Algorithm 3, the NS are polarized during the GES algorithm execution and may, by chance, “generate/create” new structures. With this new proposed mechanism it is possible to create new ways to shake a given solution, since a solution , generated from *s*, is a combination of random moves from all , as described in Algorithm 4.

A short description of the movements that will guide the GES walk through the solutions space are described next:

**Movement number of trips**—: This move increases or decreases in one unit the number of trips a truck *l* performs to a pit *i*, in which there is a compatible load equipment. Thus, in this movement a cell *n _{il}* of matrix

*N*has its value increased or decreased by one.

**Movement load**—: It consists of swapping two distinct cells *y _{i}* and

*y*of matrix

_{j}*Y*, that is, swapping the load equipments allocated to pits

*i*and

*j*, if both pits have an allocated loading equipment. When there is an allocated load equipment on only one of the pits, this movement will relocate the load equipment to the available pit. To maintain compatibility between shovels and trucks, the trips made to the pits are relocated along with the load equipments.

**Movement relocate trip from a truck**—: In this movement, two cells *n _{il}* and

*n*of matrix

_{kl}*N*are selected and one unit of

*n*is transferred to

_{il}*n*. Thus, the truck

_{kl}*l*does one trip less to pit

*i*and it does one trip more to pit

*k*. Compatibility between equipments are observed, with relocation of the trip only if there is a match between them.

**Movement relocate trip from a pit**—: Two cells *n _{il}* and

*n*of matrix

_{ik}*N*are selected and one unit of

*n*is relocated to

_{il}*n*. This move relocates one trip that the truck

_{ik}*l*performs to pit

*i*for the truck

*k*. Compatibility between equipment restrictions are respected and there is a relocation only when there is a match between them.

**Movement pit operation**—: This operation consists of removing the load equipment that is operating in pit *i*. The procedure removes all trips made to the pit *i*, leaving the shovel equipment allocated to it inactive. The equipment returns to operation only once a new route is associated to the pit.

**Movement truck operation**—: Consists of removing from operation one truck *l* that is related to a pit *i*. Thus, the cell *n _{il}* of matrix

*N*has its value set to zero.

**Movement swap trips**—: Two cells of the matrix *N* are selected and one unit of a cell is transferred to another one, which means the journey of a truck made to a pit is forwarded to another truck on another pit.

**Movement swap shovels**—: Two distinct cells *y _{i}* and

*y*matrix

_{k}*Y*have their values exchanged, that is, the load equipments operating on pits

*i*and

*k*are exchanged. Similar to the neighborhood structure

*NS*, the load equipments are exchanged, but the trips made to the pits are not. To maintain compatibility between shovels and trucks, the incompatible trips made to the pit are removed.

^{LD}### 3.4 Computational Experiments and Analysis

Computational experiments were carried out on a Pentium Core 2 Quad (Q6600) with 8 GB of RAM, and operating system Ubuntu 14.04.

A first batch of experiments, discussed in Section 3.4.1, seeks to find good parameters of population size and selection strategy rates for the OPMOP. Time-to-Target plots (TTTplots) were used in order to find a set of parameters able to find good targeted solutions. This approach was used since we are only trying to demonstrate the convergence of the GES and, if possible, achieve comparable solution with the GGVNS algorithm of Souza et al. (2010). Thus, the main focus of the experiments section is to demonstrate the ability of the proposed algorithm to self-adapt its parameters during the evolution process. This feature is discussed in Section 3.4.2. In addition, the proposed hybrid self-adaptive GES is tested on a set of standard benchmark problem instances from the literature.^{1} These test problems were the same used in Souza et al. (2010) to validate the GGVNS algorithm.

#### 3.4.1 GES Calibration Using Time-to-Target Plots

Six variants with different sets of parameters are analyzed here; Table 2 shows the parameters of these GES variants. The variants GES1-VND and GES2-VND include VND as the local search procedure for refining some solutions. In these two variants, a portion of individuals from the offspring population are chosen for the local search procedure (as suggested in Line 17 of Algorithm 1). Only a small group of those NS described in Section 3.3 is used, namely: *NS ^{LD}*,

*NS*,

^{NT}*NS*, and

^{TT}*NS*. The expensive computational cost of the local search justifies this restriction. The maximum number of application was set as for each .

^{TP}Acronym . | . | . | Selection . | VND . |
---|---|---|---|---|

GES1 | 30 | 160 | () | |

GES2 | 30 | 160 | () | |

GES3 | 100 | 600 | () | |

GES4 | 100 | 600 | () | |

GES1-VND | 30 | 160 | () | |

GES2-VND | 30 | 160 | () |

Acronym . | . | . | Selection . | VND . |
---|---|---|---|---|

GES1 | 30 | 160 | () | |

GES2 | 30 | 160 | () | |

GES3 | 100 | 600 | () | |

GES4 | 100 | 600 | () | |

GES1-VND | 30 | 160 | () | |

GES2-VND | 30 | 160 | () |

Two TTTplots experiments were performed for checking the efficiency of the proposed variants in achieving targeted solutions. Runtime distributions or TTTplots display, on the ordinate axis, the probability that an algorithm will find a solution at least as good as a given target value within a given running time, shown on the abscissa axis. These plots were first used in Feo et al. (1994). Runtime distributions have been advocated also in Ribeiro and Resende (2011) as a way to characterize the running times of stochastic algorithms for combinatorial optimization.

Aiex et al. (2007) described a Perl program to create TTTplots for measuring times that are assumed to fit a shifted exponential distribution, closely following Aiex et al. (2002). Such plots are very useful to compare different algorithms or strategies for solving a given problem and have been widely used as a tool for algorithm design and comparison.

For a better comparison among the variants, their empirical probability curves were superimposed. On the first experiment, the algorithms were applied to the test problem opm1; the target was set at 230.00 (2% of optimal value), as shown in Figure 2a. In the second one, the algorithms were applied to the instance opm8; the target was set at 164,024.00 (0.0033% of optimal value), as shown in Figure 2b. A battery of 100 executions was made and the performance ended only when the algorithm had found the target value. These times were then sorted in ascending order, and for each algorithm, were associated with the *i*-th largest running time *t _{i}*, a probability and the points plotted , for each . The results of the experiments are shown in Figures 2a and 2b.

Analyzing the empirical probability curves, it is possible to see that the variants that used the selection strategy () prevailed over the versions that used the selection (). This fact shows that the competition between parents and offspring made those individuals with a good optimality potential persist for more generations. This result is consistent with the report of Beyer and Schwefel (2002) and Herdy (1992), recommending the use of the () selection in discrete finite size search spaces.

In Figure 2a there is a total supremacy of the GES4, achieved by its selection strategy () combined with population size and . Since the initial instants of the search, the variant GES4 was able to generate better solutions than the other algorithms proposed. However, analyzing the curves of Figure 2b, one can notice that from 60 seconds and on, GES4 lost its performance, being surpassed by the variant GES2-VND, which continues to progress systematically, being the first to reach the desired target with a probability of approximately 100%.

In order to deal with the situation shown in Figure 2b, a probability experiment according to Ribeiro and Rosseti (2009) is presented. Let A1 and A2 be two stochastic search algorithms applied to the same problem and let X1 and X2 be the continuous random variables represeting the time required for algorithms A1 and A2, respectively, to find a solution as good as the given target. Ribeiro and Rosseti (2009) developed a numerical tool to calculate the probability of the runtime of the algorithm A1 being less than or equal to the runtime of the algorithm A2, that is, . This tool approximates the absolute error in the integration, by selecting the appropriate value of . The latter optimizes the resulting approximation errors, called , in order to make it sufficiently small. Using this tool to validate the analysis of the empirical experiment in Figure 2b, Table 3 was generated.

. | GES1 . | GES2 . | GES3 . | GES4 . | GES1-VND . | GES2-VND . |
---|---|---|---|---|---|---|

GES1 | 48.23% | 42.74% | 32.17% | 66.34% | 63.13% | |

GES2 | 50.16% | 43.77% | 34.19% | 65.58% | 62.65% | |

GES3 | 56.05% | 53.94% | 37.68% | 72.52% | 69.91% | |

GES4 | 65.68% | 61.72% | 59.25% | 80.28% | 78.68% | |

GES1-VND | 33.52% | 34.15% | 27.28% | 19.36% | 45.85% | |

GES2-VND | 36.83% | 37.27% | 30.02% | 21.20% | 54.15% |

. | GES1 . | GES2 . | GES3 . | GES4 . | GES1-VND . | GES2-VND . |
---|---|---|---|---|---|---|

GES1 | 48.23% | 42.74% | 32.17% | 66.34% | 63.13% | |

GES2 | 50.16% | 43.77% | 34.19% | 65.58% | 62.65% | |

GES3 | 56.05% | 53.94% | 37.68% | 72.52% | 69.91% | |

GES4 | 65.68% | 61.72% | 59.25% | 80.28% | 78.68% | |

GES1-VND | 33.52% | 34.15% | 27.28% | 19.36% | 45.85% | |

GES2-VND | 36.83% | 37.27% | 30.02% | 21.20% | 54.15% |

Analyzing Table 3, it appears that even though the variant GES2-VND has a greater probability of finding the target starting from 60 seconds of execution, the variant GES4 has a probability of 78.68% to have the runtime less than or equal to the variant GES2-VND. In addition, it’s clear that the variant GES4 outperforms all other variants.

This first set of experiments allows us to define, at least, a first set of parameters for the hybrid self-adaptive GES able to converge in a similar computational time than the one used by the literature. From now on, the variant GES4 will simply be named GES.

#### 3.4.2 Evolution Strategy Self-Adaptive Mechanism

In order to verify the effect of the maximum number of application for each NS of the OPMOP, a first batch of experiments composed of 1,800 executions, 225 for each of the eight instances, was performed with different limits *nap _{k}* (as presented in Section 2.1). Objective functions were normalized for each instance and an Analysis of Variance (ANOVA) test (Shapiro and Wilk, 1964) was done for analyzing the differences between the limits

*nap*. Figure 3 shows an effect plot with limits , , , , for . As can be noticed, the only significant difference detected, with 95% confidence level, was the worse performance of the GES with strict limits . We believe that the model could be free to adapt the NS application naturally; thus, we might have left it with a large number . However, in order to keep the disturbances slighter, we decided to keep the maximum number of application for each

_{k}*nap*as 15.

_{k}Instance opm1 was now solved and the mean values of the mutation vectors of typically 120 seconds execution are reported. The mutation operators were able to adapt its application probability during the evolutionary process, as can be seen in Figure 4.

A special relationship was detected between some NS; thus, two other plots (see Figures 5a and 5b), were generated focusing only on the specific ability of combining NS. As described in Section 3.3, the neighborhood is able to swap shovels from different active pits. Due to the different ore composition among pits, it is interesting to reallocate trips for improving the final quality of the mixture after a swap done by . The GES was able to increase the number of moves from the neighborhood at the same time that the probability of applying swaps from increased. This was an interesting fact, since we verified that improvements on the trip allocation had to be performed and was able to attend to these changes respecting the other pits, shovels, and trucks.

Another interesting behavior of the mutation operators was detected between the neighborhood , able to remove a pit from operation, and the neighborhoods that deal with trucks' trips and . Again, GES was able to take profit of free trips from the trucks now free due to the higher probability of using moves.

#### 3.4.3 Benchmark Results

Finally, in this section, the GES performance is verified against two algorithms already applied for the OPMOP solution. A batch of 30 executions was done and objective functions values were calculated for comparing the performance of the proposed model over each of the eight test problems. Both algorithms from the literature, used for comparison, are based on VNS methods. The GGVNS algorithm is a trajectory search--based algorithm and its best results and average objective function in 30 executions were also previously reported. The second one is a multi-objective algorithm, so-called G2PPLS-VNS, for which only the best results were reported.

*imp*refers to the relative average percentage values gain of the GES algorithm against the GGVNS; that is:

^{d}According to Table 4, we could verify that the proposed GES algorithm was able to generate good quality solutions with low variability, being competitive with the GGVNS algorithm. It was noticed that GES was able to slightly improve the quality of final solutions up to 0.87% and reduce the variability of these solutions up to 0.70%. Furthermore, the GES algorithm was able to find a better solution than the best one known by the literature, for instance, opm5. A similar behavior happened in the batch of experiments performed in Section 5.3.1. Thus, even though the proposal uses the same NS of both algorithms in analysis, individuals reached, guided by GES self-adaptive strategy, different solutions result in the search space.

. | GGVNS . | G2PPLS . | GES . | . | . | ||
---|---|---|---|---|---|---|---|

Instance . | Average . | Best . | Best . | Average . | Best . | imp
. ^{b} | imp
. ^{d} |

opm1 | 230.12 | 230.12 | 228.12 | 228.50 | 228.12 | 0.87 | 0.70 |

opm2 | 256.56 | 256.37 | 256.37 | 256.43 | 256.37 | 0.00 | 0.05 |

opm3 | 164064.68 | 164039.12 | 164046.32 | 164044.68 | 164031.28 | 0.00 | 0.01 |

opm4 | 164153.92 | 164099.66 | 164074.32 | 164097.61 | 164057.04 | 0.03 | 0.03 |

opm5 | 228.09 | 228.09 | 227.04 | 227.21 | 226.66 | 0.63 | 0.39 |

opm6 | 237.97 | 236.58 | 236.35 | 237.07 | 236.58 | 0.00 | 0.38 |

opm7 | 164021.89 | 164021.38 | 164018.81 | 164020.24 | 164018.22 | 0.00 | 0.00 |

opm8 | 164027.29 | 164023.73 | 164022.63 | 164022.38 | 164020.26 | 0.00 | 0.00 |

. | GGVNS . | G2PPLS . | GES . | . | . | ||
---|---|---|---|---|---|---|---|

Instance . | Average . | Best . | Best . | Average . | Best . | imp
. ^{b} | imp
. ^{d} |

opm1 | 230.12 | 230.12 | 228.12 | 228.50 | 228.12 | 0.87 | 0.70 |

opm2 | 256.56 | 256.37 | 256.37 | 256.43 | 256.37 | 0.00 | 0.05 |

opm3 | 164064.68 | 164039.12 | 164046.32 | 164044.68 | 164031.28 | 0.00 | 0.01 |

opm4 | 164153.92 | 164099.66 | 164074.32 | 164097.61 | 164057.04 | 0.03 | 0.03 |

opm5 | 228.09 | 228.09 | 227.04 | 227.21 | 226.66 | 0.63 | 0.39 |

opm6 | 237.97 | 236.58 | 236.35 | 237.07 | 236.58 | 0.00 | 0.38 |

opm7 | 164021.89 | 164021.38 | 164018.81 | 164020.24 | 164018.22 | 0.00 | 0.00 |

opm8 | 164027.29 | 164023.73 | 164022.63 | 164022.38 | 164020.26 | 0.00 | 0.00 |

## 4 Short-Term Load Forecasting Problem

The importance of load forecasting has been increasing lately and improving the use of energy resources remains a great challenge to the emerging Smart Grid (SG) systems. Energy consumption forecasting in the context of economic development of a country was highlighted by Lee and Tong (2011). SG are considered the future of power grids, able to manage production, transmission, and electricity distribution. The task of optimizing the SGs has been done mainly by using an Artificial Intelligence (AI) technique (Raza and Khosravi, 2015; Olivares et al., 2014; Rigo-Mariani et al., 2014).

In this sense, improving the calibration of a forecasting model through the aid of an evolutionary algorithm seems reasonable. Here, we verify the ability of the GES in calibrating the Hybrid Forecasting Model (HFM) proposed by Coelho et al. (2014a; 2014b; 2016). In its previous version, the model’s fuzzy rules were being calibrated by a trajectory search--based algorithm and a classic evolution strategy using a mutation matrix fulfilled with standard deviations.

We will use this problem as a didactic case of study that combines the use of nine different NS. An acceptable convergence of the model is checked and verified in Section 4.3.2.

### 4.1 Representation and Evaluation of a Solution

A solution to the HFM is represented as a matrix of continuous values indicating fuzzy rules intervals and respective weights. An example of a solution representation *s* with three columns can be seen in Figure 6. For the batch of experiments analyzing the mutation operators, we are going to use a solution with several different inputs and approximately 1,000 columns to be calibrated.

The evaluating process is simple; the solution depicted in Figure 6 is evaluated by its ability in forecasting a given validation historical load time series. The rules of the matrix are applied considering previously measured data. Results of the combination of all rules and its weights give the next point forecast. These are very similar to artificial neural network--based models (Drezga and Rahman, 1999) and fuzzy time series (Song and Chissom, 1993). Errors between each forecast and the real measured point from the validation set are calculated using quality indicators (Goodwin and Lawton, 1999).

### 4.2 Neighborhood Structures

We designed a new didactic NS able to change each cell of the fuzzy matrix *s* with a disturbance *X*. The NS will be called and is described next:

**Movement add X**——This move increases or decreases the value of a random cell of the rules and weights matrix of a solution

*s*.

From this proposed NS, we generated nine different structures with different disturbances parameters *X* : , , , , , , , , and . Some of them use the average values of the historical load time series, namely *M*, as disturbance value. The special character *B* in indicates a big value multiplied by the average *M*.

### 4.3 Computational Experiments and Analysis

Computational experiments were carried out on an Intel Core i7-3537U CPU (2.00 GHz), with 4 GB of RAM, on the operating system Ubuntu 14.04. For simplicity, the configuration achieved after the OPMOP TTTplots calibration (Section 3.4.1) and confirmed in literature instances is used here.

The dataset used to check the mutation operators was kindly provided by Liu et al. (2014). It is composed of microgrid user data from a small residential area with maximum load of 273 KW. The dataset is composed of 1,368 hourly samples for training and 672 samples used as blind validation. Section 4.3.2 reports obtained Mean Absolute Percentage Error (MAPE) for the validation set, as a way of certifying the success of the GES in calibrating the metaheuristic fuzzy model.

#### 4.3.1 Evolution Strategy Self-Adaptive Mechanism

A first batch of two-minutes training, typical for online microgrids load forecasting, was performed and the behavior of the mutation operators are discussed in Figures 7, 8, and 8.

Figure 7 plots the vector of probability multiplied by the current number of applications of each NS; average values for the population are presented for each generation. The ability of the GES mutation mechanism to regulate the probability of the NS regarding its power of disturbance is highlighted. Moves that slightly change the solution matrix *s* are the most likely to be applied, as is the case of , , and . On the other hand, was adapted and adjusted, through the generations, to be rarely applied.

Figures 8b and 8a show groups of NS aggregated regarding the impact they have in mutating the solutions.

The stem plot (see Figure 8b) shows the higher number of moves that are applied for some NS and differences between application probabilities.

Another batch with two-hour training was performed. Figure 9a shows the average values with and Figure 9b presents the results for a larger limit of 1,000. We highlight that the operators seldom present slopes with higher probabilities and more moves applications. Specially, for the neighborhood it happened in both cases, but a few generations later it converged to a steady state, returning the average values of the population to low values of application probability. Even after generation 3000, Figure 9b, we believe that the number of applications would follow the same decrease after some generations.

#### 4.3.2 GES Convergence

A batch of 30 executions for the aforementioned dataset was executed and average MAPE errors of 9.5%, 8.5%, 9.8%, and 8.2% were obtained for the 1st, 2nd, 3rd, and 4th weeks, respectively. Thus, as described by Liu et al. (2014), it is known that when the MAPE is less than 10%, the applicability of the forecast model over microgrids becomes interesting and does not increase its cost sharply. The obtained results indicate, again, the ability of the GES in calibrating the matrix of weights and rules and achieving useful forecasting models.

## 5 Unrelated Parallel Machine Scheduling Problem with Setup Times

The UPMSP-ST tackled here has as the main goal the makespan minimization. This problem has great practical and theoretical importance. It belongs to the -hard class, since it can be seen as a generalization of Parallel Machine Scheduling Problem with Identical Machines and without Setup Times (Garey and Johnson, 1979). The UPMSP-ST is found in different industry sectors, such as textile, chemicals, painting, semiconductors, and paper production (Rabadi et al., 2006).

In the UPMSP-ST there is a set of *N* jobs and a set of *M* machines, with the following features: (i) each job must be allocated to only one machine; (ii) each job *j* has a different processing time *p _{ij}* for each machine ; (iii) there is also a setup time

*S*for calibrating the machine

_{ijk}*i*after processing job

*j*and before processing job

*k*; (iv) there is a calibration time for processing the first job of a given machine .

Different approaches were used for solving the UPMSP-ST; initially, in Al-Salem (2004), the authors developed a heuristic procedure called Partitioning Heuristic and introduced the UPMSP-ST. Rabadi et al. (2006) proposed a metaheuristic for random prioritized search and described a mathematical formulation for the problem. In Ying et al. (2012), a Simulated Annealing algorithm was implemented with a smart strategy for eliminating unpromising jobs. Vallada and Ruiz (2011) analyzed two genetic algorithms and published a benchmark set of instances at the SOA website: http://soa.iti.es/problem-instances. In the works of Haddad et al. (2014), Cota, Haddad, Souza, and Coelho (2014), Cota, Haddad, Souza, and Martins (2014), and Haddad et al. (2015), the authors proposed several trajectory search--based algorithms, named AIV, AIRP, AIA, and HIVP, respectively. These methods were designed based on the ILS and VND, which obtained better results than Vallada and Ruiz (2011), and among them, AIRP had the best performance. The AIRP combines a greedy constructive procedure with the ILS and RIV (described in Section 5.2.1) metaheuristics. Also, periodically, the search is intensified and diversified by a Path Relinking (Glover, 1996) procedure. In this sense, comparing the proposed evolutionary framework against the AIRP (see Section 5.3) is reasonable.

### 5.1 Representation and Evaluation of a Solution

A solution *s* for the UPMSP-ST is represented as a vector of integers with *m* positions, where each position represents one machine. Each active machine is associated with a list containing all jobs allocated to it.

Figure 10a shows an example of a possible scheduling for an instance with two machines and six jobs. In this example, machine *M*_{1} will process jobs , and 1, in this order; and machine *M*_{2}, in turn, will process , and 6, following this order. The conclusion time of machine *M*_{1} is calculated by the expression , while machine *M*_{2} conclusion time is given by .

Figure 10b illustrates the solution representation of the example given in Figure 10a. Jobs 3, 5, and 1 are scheduled to machine *M*_{1} while the rest are scheduled to *M*_{2}.

In this single objective optimization, solution *s* is evaluated by the makespan, in other words, by the processing time of the last machine to finish its jobs.

### 5.2 Neighborhood Structures

Three well-known NS are tested in this application and they are described as follows:

**Movement swap in the same machine**—: The movement of swapping jobs in the same machine originates an neighborhood. It consists of changing the positions of two jobs that belong to the same machine.

Figure 11a illustrates the swap of jobs 5 and 6 in machine *M*_{2}.

**Movement swap between different machines**—: Similar to the , but swapping one job from one machine with another job that belongs to a different machine.

**Multiple insertion—**: It consists of relocating one job from one machine to any position of all machines.

Figure 11b shows how works; job 4 from machine *M*_{2} is transferred to machine *M*_{1}, just before job 1.

#### 5.2.1 RIV as an Optional Local Search Intensification Phase

The optional local search phase of the GES, Line 17 of the Algorithm 1, activates the proposal of Cota, Haddad, Souza, and Coelho (2014), namely RIV. The latter is inspired by the ILS with a Random VND.

The RIV uses the neighborhoods and with first and best improvements strategies, respectively, and explores, sporadically, the neighborhood for perturbing the solutions.

### 5.3 Computational Experiments and Analysis

The proposed algorithm for the UPMSP-ST was implemented in JAVA using the Netbeans 8.0.2 IDE. Computational experiments were carried out on a Core i7 (1.9 GHz) with 6 GB of RAM and Windows 7.

The method was re-implemented, following OptFrame first version in C++, and the ideas introduced in Section 2. This seems to be a good opportunity for analyzing and ensuring the convergence of the proposal in a different environment, considering new codes programmed in a similar language.

#### 5.3.1 Benchmark Results

Eight different configurations of the GES algorithm were designed, as described in Table 5. Six, among these variants, were allowed to activate RIV local search in individuals from the offspring population, picked at random (see Section 3.4.1). The other two variants, with the same population size as the ones used for the previous case studies were analyzed without the RIV intensification phase.

Acronym . | . | . | Selection . | RVI . |
---|---|---|---|---|

GES_{1} | 50 | 150 | () | |

GES_{2} | 20 | 60 | () | |

GES_{3} | 50 | 150 | () | |

GES_{4} | 20 | 60 | () | |

GES_{5} | 100 | 600 | () | |

GES_{6} | 100 | 600 | () | |

GES_{7} | 100 | 600 | () | |

GES_{8} | 100 | 600 | () |

Acronym . | . | . | Selection . | RVI . |
---|---|---|---|---|

GES_{1} | 50 | 150 | () | |

GES_{2} | 20 | 60 | () | |

GES_{3} | 50 | 150 | () | |

GES_{4} | 20 | 60 | () | |

GES_{5} | 100 | 600 | () | |

GES_{6} | 100 | 600 | () | |

GES_{7} | 100 | 600 | () | |

GES_{8} | 100 | 600 | () |

The stopping criterion adopted was the processing time, given by milliseconds, with *n* the total number of jobs, *m* the total number of machines, and three different values for *t* (10, 30, and 50). These values adopted for *t* were the same ones by Cota, Haddad, Souza, and Coelho (2014).

A batch of 30 executions was performed using 36 SOA instances. These instances involve combinations of 50, 100, and 150 jobs with 10, 15, and 20 machines. The average values of each instance were calculated and the deviations between the best results for each instance were measured. Due to the stochastic character of the search, the different GES configurations were executed 30 times for each instance and for each *t* value, as has been done already for the AIRP.

Figure 12 shows a box plot graph with average objective function values.

It is noteworthy that the selection strategy () was again able to direct the evolutionary process to better performance. We were able to obtain better solutions (10% better) than the ones found by the AIRP. This fact shows the potential of the proposal to combine NS and to find new solutions in the search space.

It should be noticed that the two variants with the intensification procedure RIV achieved better average objective function values. This fact induces that the RIV procedure was well designed by Cota, Haddad, Souza, and Coelho (2014), being able to collaborate within the search as an intensification phase. Furthermore, it gives a brief motivation for using other smart strategies, based on metaheuristic procedures, in partnership with the proposed GES.

#### 5.3.2 Evolution Strategy Self-Adaptive Mechanism

Three different sizes of instances (large, medium, and small size) were used for analyzing the mutation operators. Parameters and were kept for this analysis. We aim now at two variants, verifying the convergence with and without the RIV intensification phase. For the sake of clarity, the best results among the two independent runs (using different seeds) are depicted in Figures 13, 14, and 15.

In these figures, the lower bounds (some are optimum values) were plotted for enhancing bounds comprehension of the mutation operators' behavior. The evaluation of the best known solution can be followed across the generations, exhibited with blue lines. Values were normalized for each of the three instances by dividing them by the maximum makespan found in the first generation.

The convergence of the variant using RVI was close to the lower bound in all executions. Figure 13 shows how operators keep trying to escape from local optimum, considerably improving the best known solution until generation 150.

Figure 14 depicts an execution applied for solving a medium size instance, in which the population was guided only by the ARVNS. Analyzing Figure 14b, it can be verified that the mutation operators suffered an intensive change after generation 50. This fact might have happened because of a local trap, which the population was able to escape near generation 600 (see Figure 14a). Four other improvements were also achieved from generation 600 to 1,000.

## 6 Conclusions and Extensions

A new hybrid self-adaptive algorithm based on the concepts of evolution strategies was introduced in this article. Three different combinatorial optimization problems were used as cases of study. For each of them, algorithm convergence and mutation operators' behavior were analyzed.

The results have shown that the proposed evolutionary method is able to achieve competitive solutions. For the OPMOP, the proposed GES algorithm was compared with two other algorithms from the literature. Computational experiments emphasized the competitiveness of the GES algorithm, since it was able to improve the solutions’ quality and reduce their variability.

Even though the average performance of the GES, applied to the UPMSP-ST, was not considerably better than the literature, it was able to enhance the quality of the solutions in 10% of the analyzed problems.

The ability of adapting the probabilities of application was also verified in each of the three analyzed problems. The operators handled different NS and its application across different phases of the evolution process. The flexibility of the proposal makes it suitable for a wide area of practical applications, such as in mining, scheduling, and electric load forecasting.

As could be verified in this article, the self-adaptive ES was able to adapt the mutation operators in such a way that there is a balance between exploration and exploitation throughout the generations of the evolutionary process, being able to escape from local optima attraction basins. This ability is related to the issue of performing small changes in the solution for finding local optima, and, on the other hand, enhancing the probability and strength of the mutation operators reflected in increasing the shaking and hence escaping from local traps. This fact promoted a self-adaptive balance between exploiting an attraction basin and jumping out of it.

Further experiments should focus on how close the solutions are from local optima. Thus, future work might investigate the fitness landscape of the problems discussed here. Improving and designing novel self-adaptive mechanisms is another possible extension, as well as parameter calibration. Finally, we suggest implementing a parallel version of the GES algorithm in order to take advantage of the multicore technology available in current devices.

## Acknowledgments

The authors are indebted to the anonymous reviewers for their constructive suggestions that have helped us improve the original manuscript. The authors acknowledge CNPq (grants 306694/2013-1 and 312276/2013-3) and FAPEMIG (grants CEX PPM 772/15 and APQ-02449-14) for supporting the development of this research. Vitor N. Coelho would like to thank his brother, Bruno, for all his motivation and friendship.

## References

*Computers and intractability: A guide to the theory of*

*NP-completeness*

^{nd}Conference of Parallel Problem Solving from Nature

## Note

^{1}

Available at http://www.decom.ufop.br/prof/marcone/projects/mining.html