Unconstrained numerical optimization using real-coded genetic algorithms: a study case using benchmark functions in R from Scratch

Unconstrained numerical problems are common in solving practical applications that, due to its nature, are usually devised by several design variables, narrowing the kind of technique or algorithm that can deal with them. An interesting way of tackling this kind of issue is to use an evolutionary algorithm named Genetic Algorithm. In this context, this work is a tutorial on using real-coded genetic algorithms for solving unconstrained numerical optimization problems. We present the theory and the implementation in R language. Five benchmarks functions (Rosenbrock, Griewank, Ackley, Schwefel, and Alpine) are used as a study case. Further, four different crossover operators (simple, arithmetical, non-uniform arithmetical, and Linear), two selection mechanisms (roulette wheel and tournament), and two mutation operators (uniform and non-uniform) are shown. Results indicate that non-uniform mutation and tournament selection tend to present better outcomes.


Introduction
Numerical Optimization problems exist widely in di erent areas of science research and engineering practice (Zang et al., 2018), i.e., it is an essential tool in decision science and the analysis of physical systems (Nocedal and Wright, 2006). In other words, it is a tool for solving practical problems devised by many variables and no constraints, also known as unconstrained problems. Its primary purpose is to discover the best values for design variables and/or objective functions that are not known precisely (Muc and Sanetra, 2017). In general, unconstrained problems can be classi ed into two categories: test functions and real-world problems. Test functions are arti cial problems and can be used to evaluate the behavior of an algorithm in sometimes diverse and challenging situations (Jamil and Yang, 2013). On the other hand, real-world problems originate from di erent elds such as physics, chemistry, engineering, mathematics, etc. In this work, we will focus on ve test functions, also known as benchmark functions. A set of realworld problems can be seen in Averick et al. (1992). Those functions have been used to test di erent kind of algorithms as we can see in Borges et al. (2018), Maucec and Brest (2018), Karaboğa and Kaya (2018), Cavalca and Fernandes (2018), etc.
There are several optimization techniques for solving unconstrained numerical problems.
The traditional ones aim to discover the optimum solutions of continuous and di erentiable functions, i.e., they use analytical methods and calculus to locate the best solutions. In fact, the classical methods are fast; however, they are limited because they can only deal with unconstrained function and a small number of variables. Moreover, practical applications usually deal with non-di erentiable functions. Thus, evolutionary algorithms (Eiben and Smith, 2007) appear as a viable solution for optimizing constrained and nondi erentiable functions.
In this context, the idea of this work is to present a tutorial on Genetic Algorithms in Numerical Optimization using benchmark function with a study case in R language (RStudio Team, 2018). Then a question can appear: why not using an R package? The main reason is that when we use a package, such as genalg (Willighagen and Ballings, 2015), we are limited to those features o ered by the package. Particularly in the genalg package, the user has no control on genetic operators whatsoever. The only control provided by the referred package is mostly concerning parameters, i.e., the user cannot choose di erent crossover or mutation operators. Regarding the GA package (Scrucca, 2017), which is a more generic and exible package includes other evolutionary algorithms such as Di erential Evolution, present an advantage of using di erent genetic operators and parallel algorithms. However, having the code made from scratch, the user can quickly implement di erent operators or even create hybrid algorithms that are tter to the problem being solved.
Thus, this tutorial is divided as follows: Section 2 presents some basics on numerical optimization and benchmark functions; Section 3 shows the theory of real-coded genetic algorithms and their operators; Section 4 explains important concepts in R that are essential to understand the code; Section 5 implements all GA concepts in R; Section 6 illustrates how the GA works in ve benchmark functions; nally, Section 7 draws the conclusions of this tutorial.

Numerical Optimization and Benchmarks
The unconstrained optimization aims to minimize or maximize an objective function that depends on real variables, with no restrictions at all on the values of these variables (Nocedal and Wright, 2006). Mathematically, it is min or maxf(x), where x ∈ R n and n ≥ 1. Thus, a solution x * is a global solution of a Regardless of the kind of optimization, if we want to use a GA for this kind of problem, it is mandatory n > 1. Actually, n regards to the dimensionality of the search space, which is an important factor in the problem complexity, since the higher the dimension, the higher the probability of getting trapped in a local optima (Cortes et al., 2012). A study of the dimensionality problem and its features was carried out by Friedman (1994).
Two other properties are essential in numerical optimization: separability and multi-modality. The separability concerns the possibility of dividing f(x) into two or more functions. Consequently, non-separable functions are harder to optimize then separable ones. Multi-modality regards to the existence of many local optima. In this context, non-separable and multimodal functions are harder to solve than the other ones.
We will test our code using four unconstrained continuous numerical benchmarks functions: Rosenbrock (1960), Griewank (1981), Ackley (1987), Schwefel (1981), and Rahnamyan et al. (2007) as presented in Table 1. Table 2 presents the benchmarks properties (Separability, Modality, and Di erentiability), the domain, and the global optima. The domain is a constraint for each gene, i.e., the lower and upper bounds. The optimal solution is the minimum value that the benchmark can reach. The separability represents if the function is separable, i.e., if the function can be split into two or more functions. In other words, a function of p variables is called separable, if it can be written as a sum of p functions of just one variable (Boyer et al., 2005). Finally, the modality regards to the existence of many local optima. In this context, non-separable and multi-modal functions are harder to solve than the other ones. Fig. 1 shows the pseudo code for a GA. Firstly, the GA creates a random population, then evaluates it to select individuals undergoing genetic operators. Usually, methods such as the roulette wheel or tournament, for example, select the new population. Afterward, the crossover exchanges information (genes) between two parents based on the probability of mutation (pc), creating one or more o spring. Then, the mutation operator can change zero or more genes on  I f ( best ( f i t ) > best ( f i t ' ) ) 12

Real-Coded Genetic Algorithms
Swap( ) 13 End-I f 14 End-I f 15 End-While

Representation
There are two main representations of genetic algorithms: binary-coded and real-coded. In the binary-code representation, an individual or chromosome, which is a possible solution to the problem being solved, is represented by a vector of {0, 1}. On the other hand, as expected, real numbers devise a real-coded chromosome as presented in Fig. 2. Because we are dealing whit real numbers, each gene requires a domain constraint represented by a lower and an upper bound, respectively (LB and UB). In other words, assuming that c i is a gene within a chromosome i, we have LB i ≤ c i ≤ UB i .

Crossover
The main purpose of the crossover operator is to exchange information (genes) between parents, creating one or more o spring.
In real-coded representation, there are several ways of doing that. A thorough list of operators can be seen in Herrera et al. (1998). In this tutorial, we describe four of them: simple, arithmetical, non-uniform arithmetical, and Linear. In the simple crossover, a cut point is randomly chosen, then the o spring are formed, making a combination of parts. Considering that p1 = c 1 1 , c 1 2 , . . . , c 1 n and p2 = c 2 1 , c 2 2 , . . . , c 2 n are two parents, and j is the cutting point, the rst o spring is (1) and (2) illustrates how to perform the arithmetical crossover for two descendants, in which r is a random number in the range [0, 1].
The di erence between arithmetical and nonuniform arithmetical crossover relies on the fact that the variable r is not random anymore but computed by dividing the current generation t by the maximum number of generations (Tmax) as shown in Eq. (3).
Finally, the Linear Crossover creates three o spring as presented in Eqs. (4) to (6). The rst one is similar to the arithmetical crossover with r = 0.5. The other ones explore the outer limits of c 1 i and c 2 i , respectively.
It is essential to observe that not every individual from the selected population undergoes mutation. Only those whose probability is less than the probability of crossover, pc, participate in the operation.

Mutation
The mutation operator changes genes from a chromosome. Concerning real-coded individuals, the uniform mutation, also known as a random mutation, randomly replaces a gene inside its domain based on a parameter known as probability of mutation pm.
Another well-known mutation in the literature is the non-uniform mutation operator, in which the selected genes are mutated according to Eqs. (7) and (8), in which s i is the gene being mutated, t represents the current generation, LB and UB are the lower and the upper bound of the variable i, respectively, r is a random number between 0 and 1, T is the maximal number of generations, and b is the degree of dependency (usually b = 5).
The non-uniform mutation is one of the operators responsible for the ne-tuning capabilities of the system (Michalewicz, 1999).

Important Concepts in R
The purpose of this chapter is not to provide a thorough vision of R programming but gives concepts that are essential to understanding the code in the next sections. Details of how to program in R can be found in Lander (2015) and Crawley (2012) books.
The rst essential concept is the notion of indirect indexing. Programming languages, such as C and Java, access elements in a matrix using a speci c index devised by row and column. If the programmer wants to access a set of elements, it is necessary to use a for loop and work on them element-by-element. In R, the indexes of a matrix can also be a matrix. For example, suppose that we have to replace some elements obeys a condition by a random number in a matrix Mat, Fig. 3 illustrates how to perform this task. The which() function returns a matrix containing two columns (row and column) of the elements that satisfy the condition. Then, a vector operation (line 2) replaces the corresponding elements. We have to note that the number of elements created by the runif() function has to be the same number of rows in idx, therefore, we have to use the nrow() function. Fig. 4 illustrates how the operation works.

Figure 3: Indirect indexing example
The other concept we have to know how to deal with is the logical indexing. The concept is quite similar to the indirect indexing; however, in this case, all indexes are logical values. Let us suppose that we have a matrix Mat of integer numbers randomly created as presented in Fig. 5. Then we want to replace all values less than 10 with the lower bound 10. The instruction idx < -Mat < 10 returns a matrix in which all positions obeys the condition are true. Afterward, all values are replaced. This operation makes things easier when the domain is the same for all genes, which is common in numerical optimization. On the other hand, if the domain is di erent for each gene, then we have to use the instruction which(idx, arr.ind = TRUE) to locate the positions that ful ll the condition.
The third important concept we have to tackle is called group operations or group functions. This kind of instruction, as the name suggests, executes in a group of data. Moreover, we preferably execute group operations instead of for loops because the rst one is usually faster than loops. The rst set of group functions is: sum(), mean(), and sd(). These functions receive a vector or a matrix as a parameter and return the sum, the mean, and the standard deviation of the entrance data, respectively. If the parameter is a matrix and ones wants to perform row or column based operations, we have to use, for example, the functions rowSum() or colSum().
On the other hand, if we want to perform row or column-based operations using a pre-existed function, we have to use the function apply(), which the syntax is apply (obj, margin, function, parameters), where obj is a data structure, usually, a matrix, margin sets the kind of operation (1 -row-based, 2 -column based), function is an implemented function, and parameters are the parameters required by the implemented function. An indispensable extension of this function is the lapply() function in which obj has to be a list. Next, we present some general remarks about R that are important to the correct implementation of GAs.

Remarks
• Vectors and matrices start with 1 instead of 0; • Operations between matrices are element-wise. If a traditional multiplication matrix is required, we have to use the symbol %*%; • Avoid for loops;

Initializing
Let us start implementing each function separately in the following order: initialize population, selection, crossover, and mutation. Thus, Fig. 6 shows the code that initializes the population, in which each gene is within the domain [lb, ub]. Then, we evaluate the population using the apply() function that receives the objective function as a parameter. Finally, the initialization function returns a list containing the rst population and the tness for each chromosome, corresponding to lines 1 and 2 from

Selection
The next step is the selection method that chooses which chromosomes will try to participate in the crossover stage. Two selection mechanisms, roulette wheel and tournament, are presented in Figs for (   In the roulette function, the variables new.pop and new. t will contain the new population and its tness, respectively. Then q will contain the cumulative probability matrix that is the wheel. Afterward, the vector r will contain one random number in the domain [0, 1] for each chromosome, and the for loop will check which one will form the temporary population undergoes the crossover process. In the selection by tournament method, the variable t.size is responsible for setting the tournament size, i.e., how many individuals will participate in each round. Then the winner of each round is chosen to undergo mutation. It is a simple but e ective form of selection.

Crossover
The next step in the GA algorithm is to apply the crossover operation on the new population. Fig. 9 shows the code for this task. The main parameters of this function are the selected population and pc (probability of crossover). Why is pc so important? Because not every single chromosome will undergo crossover, only those ones which a random number r is less than pc, i.e., r < pc. Lines 2 to 5 perform this selection. Then the for loop produces two o spring using the cross point as a divisor, returning only the population because the evaluation is necessary only after mutation.  The arithmetical crossover, also known as uniform arithmetical crossover, is an interesting crossover method because it explores the search space between two genes and does not create invalid individuals, i.e., all genes are within their domain. Fig. 10 presents the implementation in which we can see that l controls how further from a gene the crossover goes, i.e., closer to the gene of the rst or second parent.
The arithmetical crossover, also known as a uniform arithmetical crossover, is an interesting crossover method because it explores the search space between two genes and does not create invalid individuals, i.e., all genes are within their domain. Fig. 10 presents the implementation in which we can see that l controls how further from a gene the crossover goes, i.e., closer to the gene of the rst or second parent.
Finally, linear crossover generates three individuals as previously explained; however, only the best two will form the next population. Also, because of the Eqs. (5) and (6), we have to control the boundaries of the new individuals. Thus, the instruction nrow(idx) > 0 means that a gene is out of its limits. After correcting the chromosome boundaries, the concern is how to select the best two chromosomes elegantly. The answer is 1 arith . crossover <-function ( lb , ub , pop   for ( Fig. 13 shows the listing for the mutation operator that receives two essential parameters: the new population and the probability of mutation pm. Then, line 4 identi es which genes, row, and column, will undergo mutation in the whole population. Finally, line 7 evaluates the new population. The next mutation operator is the non-uniform mutation, which is based on the number of iterations or generations since as iterations go on the mutation size goes down. Fig. 14 shows the mutation operator in R, in which we have to expand the boundaries using the variables tmp.lb and tmp.ub, because these vectors can be larger than the domain vectors lb and ub. Then we use the logical indexing to decide which equation to use (θ = 0 or θ = 1 from Eq. (7)) in order to perform the mutation. 1 n. uniform . mutation <-function ( func , lb , ub , pop , pop . size , dimension , 2 pm,max. i t , i t ){ 3 r <-matrix ( runif (pop . size*dimension ) ,nrow=pop . size ) 4 f i t n e s s <-rep (NA, pop . size ) 5 idx <-which( r < pm, arr . ind=TRUE) 6 values <-pop[ idx ] 7 r <-sample( c ( 0 , 1 ) , nrow( idx ) , replace = TRUE) 8 tmp. idx <-r == 0 9 tmp. lb <-lb [ idx [ , 2  f i t n e s s <-apply (pop , 1 , func ) 17 return ( l i s t (pop = pop , f i t = f i t n e s s ) ) 18 } 19 20 delta <-function ( i t , max. i t , y , b=5){ 21 r <-runif ( length (y ) ) 22 y <-y * ( 1 -r^((1-i t /max. i t )^b ) ) 23 return (y) 24 } Figure 14: Non-Uniform mutation

Main Function and Elitism
Now, it is time to get all functions together as presented in Fig. 15. As we can see, the code is pretty similar to the algorithm presented previously in Fig. 1. To test other operators, we only need to use the respective functions replacing the symbol "#". Moreover, if someone wants to test everything at once, it is trivial to add control structures to execute all operators.
1 GA <-function ( func , lb , ub , pop . size = 10 , 2 dimension = 10 , max. i t = 100 , 3 pc = 0.6 , pm = 0.005 , sel , t . size = 4 , 4 elitism = TRUE){ 5 tmp. pop <-matrix ( rep (NA, pop . size*dimension ) , 6 nrow = pop . size ) 7 i n i t <-i n i t . population ( func , lb , ub , pop . size , dimension) 8 pop <-i n i t$pop 9 f i t n e s s <-i n i t$ f i t 10 for (   Regarding the elitism, having the nal population, we have to identify whether the elitism is set or not. If so, we have to check in which population the best individual is. In case the best individual is in the new population, the current one replaces the old one. Otherwise, we swap the worst individual from the current population by the best one from the previous one; thus, we guarantee that the best individual will always be in the population.

Benchmarks
Those benchmarks previously presented in Section 2 are implemented in Fig. 16. It is important to remark that we did not use any loop, as suggested in Section 4. Even though some R packages implement these functions, the main idea is to program everything from scratch.

Experiments
All experiments have been conducted in a Windows 10, 64 bits, 16 GB of RAM, 500 GB of SSD, R version  (Project, 2018), and RStudio 1.0.136 (RStudio Team, 2018). The GA con guration is: pc = 0.6, pm = 0.01, populatiosize = 50, and dimension = 30. In the case of tournament selection, we use four individuals in the competition. Further, we set a seed for random numbers; thus, the result of the experiment will be the same in all computers in terms of the quality of solutions. Furthermore, all experiments have been performed using 30 trials; then, we are able to present the best result, the mean, and the standard deviation. Fig. 17 presents the main script used for testing the code. Firstly, we load all functions presented previously utilizing the instruction source(). The le 'crossover.R', for instance, contains all crossover operators. The other R les follow the same principle. Either, we automatized the test only for the benchmark functions, which are stored in a vector of lists. Changing the code for executing all operators or adding new operators will be a trivial task. Tables 3 to 6 show the results for 1000 iterations with uniform mutation and roulette wheel, uniform mutation and tournament, non-uniform mutation and roulette wheel, and non-uniform mutation and tournament, respectively.
As we can see, the best combinations for those parameters is the nonuniform mutation and tournament selection, possibly because of the ne-tuning capability of the nonuniform mutation. On the other hand, Schwefel function presented the best combination using uniform mutation and tournament. Other combinations using arithmetical crossover presented good results in Griewank, Ackley, and Alpine functions; however, they are not the best ones.

Conclusions
This tutorial presented how to implement Genetic Algorithms to solve unconstrained numerical optimization problems in R from scratch. As we could see, the code is concise and straightforward, for ( Figure 17: Testing allowing anyone to implement new GA variations or even hybrid algorithms. We explore the ability to perform vectorial operations and group functions in R. Either, we implemented the GA code to be easily extended to functions in which all genes have di erent domains in the same individual. Implementing the code using a single scalar to control the boundaries of each gene could simplify parts of the code by using logical indexing; however, extend it for multiples domains in the same chromosome would be much harder.
We expect that this tutorial is an incentive to those who want to explore the numerical optimization capability of GAs and/or to those who crave to enter in the R language world. In this context, we also incentive all interested people to implement di erent operators or apply the presented code to di erent kind of numerical problems, especially those involving the optimization of real-world problems.   Table 5: Results of the optimization of the ve benchmarks functions using the following parameters: it = 1000, non-uniform mutation, roulette wheel