A Hybrid Genetic, Differential Evolution Optimization Algorithm

Abstract

This chapter presents a heuristic evolutionary optimization algorithm that is loosely based on the principles of evolution and natural genetics. In particular, this chapter describes an evolutionary algorithm that is a hybrid of a genetic algorithm and a differential evolution algorithm. This algorithm uses an elitist, ranking, random selection method, several mutation methods and both two level and three level Taguchi crossover. This algorithm is applied to 13 commonly used global numerical optimization test functions, including a spherical, three hyper-ellipsoid, the sum of different powers, Rastrigin’s, Schwefel’s, Griewank’s, Rosenbrock’s valley, Styblinski-Tang, Ackley’s Path, Price-Rosenbrock, and Eggholder’s functions. This algorithm is applied 1000 times to each of the 13 test functions, and the results shows that this algorithm always converges to each of the 13 test function’s global minimum.

1. Introduction

Optimization algorithms are systems that determine an optimal set of parameters that minimize or maximize a cost, or objective, function subject to constraints. Optimization applications are common in engineering and other scientific and mathematical fields. For a typical engineering optimization application, a cost, or objective, function mathematically describes a metric of the error between a desired performance and actual performance over a constrained solution space. For such applications, optimization algorithms would determine an optimal set of parameters that minimize the cost function subject to physical constraints, such as the optimal parameters result in a stable system. As computing power has increased, many multi-modal optimization problems have been solved using heuristic evolutionary optimization algorithms. An evolutionary algorithm is an optimization search algorithm that is loosely based on the principles of evolution and natural genetics and uses operators such as reproduction, selection, recombination and mutation [1]. Popular evolutionary algorithms include genetic algorithms [2], differential evolution [3], particle swarm optimization [4], simulated annealing [5] and colony optimization [6, 7]. Although no algorithm can solve all types of optimization problems [8, 9], genetic algorithms and differential evolution algorithms have become popular in engineering optimization applications because these method are simple, effective and flexible.

Because no algorithm can solve all types of optimization problems [8, 9], hybrid algorithms that combine the elements of an evolutionary algorithm with one or more evolutionary algorithms or search algorithms have been developed and have been shown to be effective search algorithms [10]. Because genetic algorithms and differential evolution algorithms have become popular in engineering optimization applications, this chapter presents a hybrid genetic, differential evolution algorithm. The algorithm uses an elitist, ranking, random selection method. Elitist selection methods assure the survival of the fittest individual, which is the candidate solution with the best optimization criterion cost, during the selection process. The fittest individual is also assured selection in all recombination and mutation operations. Except for the fittest individual which is guaranteed selection, the candidate solutions that survive the selection process are randomly selected for a differential evolution operator to improve convergence, a differential evolution mutation operator to improve diversity and a recombination operator that improves both convergence and diversity. The selection probabilities for the mutation and recombination operators are dynamic and change each generation, or algorithm iteration, to maintain a constant population size. After generating new candidate solutions using these operators, the new candidate solutions are added to the set of candidate solutions that survived the selection process. Except for the fittest individual (which is guaranteed selection), candidate solutions are randomly selected for Taguchi crossover [11] which is an effective recombination operator that creates near optimal new candidate solutions from two or more parent candidate solutions. Section 2 of this chapter describes the basic elements of genetic and differential evolution algorithms. Section 3 describes this chapter’s algorithm in detail. In Section 4, this algorithm is applied to 13 commonly used global numerical optimization test functions, including a spherical, three hyper-ellipsoid, the sum of different powers, Rastrigin’s, Schwefel’s, Griewank’s, Rosenbrock’s valley, Styblinski-Tang, Ackley’s Path, Price-Rosenbrock, and Eggholder’s functions.

2. Elements of genetic and differential evolution algorithms

Genetic algorithms and differential evolution algorithms are evolutionary algorithms that typically define objective functions so that the set of parameters being optimized are represented in a vector [3]. Parameter constraints are implemented by restricting the available solution spaces for each parameter in the vector. The basic design strategy for such genetic and differential evolution algorithms is to determine evolutionary operators that balance the algorithm’s ability to both effectively search the solution space and converge to an optimal solution.

Genetic algorithms and differential evolution algorithms typically begin by randomly selecting $K\left(0\right)$ vectors, called candidate solutions or individuals, in the solution space. This initial set of candidate solutions is called the initial population. A subset of $M\left(1\right)$ of these initial candidate solutions, or individuals, are selected as a function of their fitness, or cost, when evaluated with respect to the optimization criterion. Although many selection operators exist [12, 13], selection operators are typically designed to select a subset of $M\left(1\right)$ candidate solutions from the population in a manner that should improve the overall fitness of the population. The candidate solutions that survive the selection operator form the initial mating pool of $M\left(1\right)$ individuals where $1\le M\left(1\right)\le K\left(0\right)$. Candidate solutions from the mating pool are then selected for recombination and mutation. Recombination operators create new candidate solutions by combining vector elements from two or more of the candidate solutions from the mating pool. Mutation operators in genetic algorithms create new candidate solutions by altering elements of candidate solution vectors. In differential evolution, mutation operators create new candidate solutions by using vector operations, such as addition and scaling, on two or more of the candidate solutions from the mating pool. Regardless of the algorithm, the new candidate solutions created by recombination and mutation are added to the mating pool to create a new population of $K\left(1\right)$ candidate solutions. This process of selection, recombination and mutation forms one generation, or iteration, of a genetic or differential evolution algorithm. This process iterates until a convergence criteria is met and an optimal solution is determined. A generic genetic algorithm or differential evolution algorithm can be summarized as follows:

Initialize population {$K\left(0\right)$ candidate solutions)}

$n\leftarrow 0$

   repeat

      $n\leftarrow n+1$

      Select $M\left(n\right)$ candidate solutions using a selection operator

      Generate new candidate solutions using a mutation operator

      Generate new candidate solutions using a recombination operator

      Add the new $K\left(n\right)-M\left(n\right)$ candidate solutions to the mating pool

   until convergence condition is met

where $n$ is the algorithm’s iteration number.

Genetic operators such as recombination and mutation generate a combination of new candidate solutions that can be either similar or diverse from the candidate solutions in the mating pool. Controlling the ratio of the diversity and similarity of new candidate solutions added to a population each generation is a fundamental design parameter of any search algorithm [14]. Creating diversity, or exploration, is the process of generating candidate solutions that lie in previously unevaluated regions of the search space. Creating similarity, or exploitation, on the other hand, is the process of generating candidate solutions within a neighborhood of previously visited points so as to converge to an optimal point in the neighborhood. Exploration and exploitation are typically conflicting processes of a search algorithm because a lack of diversity can result in a population converging to a local minima or maxima and a lack of similarity can impede convergence. Therefore, every search algorithm needs to design an effective ratio of exploration and exploitation of a search space. In general, an optimal ratio of diversity and similarity is not only dependent on the search algorithm but also the cost, or objective, function. For example, determining the optimal solution of a unimodal objective function typically requires less exploration than determining an optimal solution of a multi-modal objective function that typically requires more exploration. Also, different generations, or iterations, of a search algorithm typically have a different optimal ratios of exploration and exploitation. For example, a search algorithm’s early generations require more exploration than exploitation until the neighborhood of the optimal solution is found. After the neighborhood of the optimal solution is found, a search algorithm’s generations require more exploitation and less exploration. Therefore, the goal of any search algorithm is to design a ratio of adding diverse new candidate solutions and similar new candidate solutions to each generation so that the algorithm can effectively determine optimal solutions for different types of objective functions.

Genetic algorithms and differential evolution algorithms typically use three operators, selection, mutation and recombination, for controlling the ratio of adding diverse and similar new candidate solutions to their populations. Selection operators control the ratio of exploration and exploitation by varying the selection process. A selection operator that is designed to select the most fit candidate solutions, the candidate solutions with the best costs when evaluated with respect to the optimization criterion, biases the selection process away from exploration and towards exploitation. A selection operator that is designed to select the least fit candidate solutions biases the selection process away from exploitation and towards exploration.

Mutation operators for a typical genetic algorithm or differential evolution algorithm randomly modify individuals from the mating pool to increase the diversity of a population. As a result, a typical mutation operator increases the exploration of unevaluated regions of the search space. However, some mutation operators only slightly alter individuals from the mating pool. In such cases, these types of mutation operators can be classified as an exploitation operator because most of the mutated individual is preserved and the mutated individual still remains in the neighborhood of the parent candidate solution. A recombination operator for a typical genetic algorithm or differential evolution algorithm combines two or more parent individuals, or candidate solutions, from the mating pool to generate a new and possibly more fit candidate solution. As a result, a typical recombination operator generates new candidate solutions within a neighborhood of the parent candidate solutions. From this perspective, a recombination operator increases exploitation and improves the convergence of the algorithm. However, when recombination uses a mutated candidate solution, recombination can create a new candidate solution in a previously unevaluated region of the search space. In such cases, recombination operators can improve exploration of the search space. In most cases, mutation operators improve a population’s diversity and recombination operators improve an algorithm’s convergence rate; however, in practice, the combination of selection, mutation and recombination determines the ratio of exploration and exploitation in both genetic algorithms and differential evolution algorithms.

3. A hybrid genetic, differential evolution algorithm

This hybrid genetic, differential evolution algorithm determines a best solution with respect to an optimization criterion that has a solution space which is the subset of an $N$ dimensional hyper-rectangular solution space although the algorithm can be adapted for other types of $N$ dimensional spaces. The algorithm generates an initial population of $K\left(0\right)$ candidate solutions by randomly selecting $K\left(0\right)$ candidate solutions within the solution space. Each candidate solution in the initial population is evaluated by an optimization criterion and ranked. Except for the top ranked candidate solution that is assured selection, a subset of these candidate solutions are randomly selected as a function of their rank. The surviving $M\left(1\right)$ candidate solutions solutions from the initial mating pool, where $1\le M\left(1\right)\le K\left(0\right)$, are randomly selected and employed in creating new candidate solutions using a differential evolution operator to improve convergence, a differential evolution operator to improve diversity and a crossover operator that improves both convergence and diversity. To ensure that the new candidate solutions lie in the solution space, one of two methods is used to move any unfeasible candidate solutions into the solution space. All the new candidate solutions are then added to the mating pool and candidate solutions from this set are randomly selected for Taguchi crossover, a type of recombination operator. After Taguchi crossover, the population should have an average of approximately $K\left(0\right)$ candidate solutions. After the Taguchi crossover operation, each candidate in the population is evaluated, ranked, and randomly selected for the next iteration. The algorithm can be summarized as follows:

Initialize population {$K\left(0\right)$ candidate solutions)}

$n\leftarrow 0$

   repeat

      $n\leftarrow n+1$

      Ranking and stochastic selection {$M\left(n\right)$ candidate solutions}

      Differential evolution operator to improve convergence

      Differential evolution mutation operator to improve diversity

      Recombination operator to improve convergence and diversity

      Ensure new solutions are in the solution space

      Taguchi crossover {$K\left(n\right)$ candidate solutions}

   until convergence condition is met

where $n$ is the algorithm’s iteration number.

3.1 Population initialization

To generate the initial population of $K\left(0\right)$ candidate solutions, candidate solutions are randomly selected so that the population is uniformly distributed over the solution space. For a hyper-rectangular solution space, an initial population of candidate solutions can be generated in Matlab using

$$\begin{array}{l}G\left(0\right)=\mathit{\text{kron}}\left({\mathbf{X}}_c,\mathit{\text{ones}}\left(1,K\left(0\right)\right)\right)\\ +\mathit{diag}\left({\mathbf{X}}_s\right)\ast \left(\mathit{rand}\left(N,K\left(0\right)\right)-0.5\right);\end{array}$$E1

where $G\left(0\right)$ is a matrix containing the $K\left(0\right)$ initial candidate solutions, ${\mathbf{X}}_c$ is a vector containing the solution space’s center in rectangular coordinates, ${\mathbf{X}}_s$ is a vector containing the solution space’s size for each dimension in rectangular coordinates, and $N$ is the solution space’s dimension. For example, for a two dimensional ($N=2$) hyper-rectangular solution space of ${\left[0-1\right]}^T\le \mathbf{x}\le {\left[\begin{array}{cc}6& 1\end{array}\right]}^T$, ${\mathbf{X}}_c={\left[\begin{array}{cc}3& 0\end{array}\right]}^T$ and ${\mathbf{X}}_s={\left[\begin{array}{cc}6& 2\end{array}\right]}^T$.

This selection of initial candidate solutions can be adapted for other types of $N$ dimensional spaces. For example, the initial population of $K\left(0\right)$ candidate solutions for a hyper-ellipsoid solution space can be generated in Matlab using

$$\begin{array}{l}G\left(1:2:N,:\right)=\mathit{\text{kron}}\left({\mathbf{X}}_r,\mathit{\text{ones}}\left(1,K\right)\right).{}^{\ast }\mathit{rand}\left(\mathit{\text{ceil}}\left(N/2\right),K\right);\\ G\left(2:2:N,:\right)=2^{\ast }{\mathit{pi}}^{\ast }\mathit{rand}\left(\mathit{\text{floor}}\left(N/2\right),K\right);\end{array}$$E2

where ${\mathbf{X}}_r$ is a vector containing the the hyper-elliptical solution space’s radii, the terms, $G\left(1:2:N,:\right)$, represent the magnitudes of each candidate solution and the terms, $G\left(2:2:N,:\right)$, represent their respective phases. If the solution space is not centered at the origin, then candidate solutions of the form ${\left[\begin{array}{cc}r& \theta \end{array}\right]}^T$ and centered at ${\left[\begin{array}{cc}0& 0\end{array}\right]}^T$ can be moved to ${\left[\begin{array}{cc}R& \mathrm{\Theta }\end{array}\right]}^T$ and centered at ${\left[\begin{array}{cc}{r}_0& \gamma \end{array}\right]}^T$ using the transformations,

$$\begin{array}{l}R=\sqrt{r^2+r_0^2+2{\mathit{rr}}_{0}cos\left(\theta -\gamma \right)}\\ \mathrm{\Theta }=arctan\left(\frac{rsin\left(\theta \right)+r_{0}sin\left(\gamma \right)}{rcos\left(\theta \right)+r_{0}cos\left(\gamma \right)}\right).\end{array}$$E3

If appropriate, these initial candidate solutions could be converted to rectangular coordinates using

$$x_k=Rcos\mathrm{\Theta }\text{and}{x}_{k+1}=Rsin\mathrm{\Theta }.$$E4

3.2 Ranking and stochastic selection

This algorithm uses an elitist, linear ranking, random selection method. Because the selection operator is elitist, the fittest individual, or the candidate solution vector with the best optimization criterion cost, is guaranteed to survive the selection process. Elitist selection algorithms can increase an algorithm’s exploitation and therefore increase the algorithm’s ability to converge, especially when steady-state misadjustment is significant [15]. Linear ranking selection methods evaluate each candidate solution by the cost function and rank the candidate solutions according to their costs [16, 17]. Starting with the candidate solution with the best cost, each candidate solution is assigned a selection probability in linearly decreasing increments so that all candidate solutions have a nonzero probability of selection. This method of selection allows diverse candidate solutions that might contain useful vector elements but have a poor cost to survive the selection process. This can improve an algorithm’s exploration and prevent the algorithm from converging in a local minima or maxima.

The selection operator is the first operation performed for each generation, or iteration of the algorithm. At the start of the algorithm’s nth iteration, the selection operator evaluates the $K\left(n\right)$ candidate solutions, ${\mathbf{x}}_k\left(n\right)$, with respect to the cost function, $J$, and ranks the candidate solutions according to their cost. For a minimization problem, the ranked candidate solutions are sorted from highest cost to lowest cost and are assigned consecutive integers from 1 to $K\left(n\right)$ so that ${\mathbf{x}}_1\left(n\right)$ is candidate solution with the highest, or worst, cost and ${\mathbf{x}}_{K\left(n\right)}\left(n\right)$ is assigned to the candidate solution with the lowest, or best, cost. After ranking, each candidate solution is assigned a selection probability, $P\left({\mathbf{x}}_k\left(n\right)\right)$, so that

$$P\left({\mathbf{x}}_k\left(n\right)\right)=\sum _{m=1}^k\mathrm{\Delta }{p}_m$$E5

where

$$\mathrm{\Delta }{p}_m=\frac{1}{K\left(n\right)}\left[{\eta }^{-}+\left({\eta }^{+}-{\eta }^{-}\right)\frac{m-1}{K\left(n\right)-1}\right],$$E6

${\eta }^{+}$ is a constant, and ${\eta }^{-}$ is a constant that is selected so $P\left({\mathbf{x}}_1\left(n\right)\right)={\eta }^{-}/K\left(n\right)$ which is the selection probability of the worst candidate solution [16, 17].

Because this algorithm uses an elitist selection method, the best candidate solution is assured survival during the selection process which implies that

$$P\left({\mathbf{x}}_{K\left(n\right)}\left(n\right)\right)=1.$$E7

Substituting Eq. (6) into Eq. (5), and the resulting equation into Eq. (7),

$$\sum _{m=1}^{K\left(n\right)}\frac{1}{K\left(n\right)}\left[{\eta }^{-}+\left({\eta }^{+}-{\eta }^{-}\right)\frac{m-1}{K\left(n\right)-1}\right]=1.$$E8

Solving Eq. (8), an elitist selection method requires that

$${\eta }^{+}=2-{\eta }^{-}$$E9

where $0<{\eta }^{-}<{\eta }^{+}$.

The set of surviving candidate solutions are referred to as the mating pool. After this selection method, the mating pool’s mean size is

$$E\left[M\left(n\right)\right]=\sum _{k=1}^{K\left(n\right)}P\left(x_k\left(n\right)\right)=\frac{2{\eta }^{-}+{\eta }^{+}}{6}\left(K\left(n\right)+1\right)$$E10

where $E$ is the expectation operator and $M\left(n\right)$ is the number candidate solutions in the mating pool after selection during the nth iteration. Because this algorithm uses an elitist selection method, Eq. (10) can be simplified by substituting Eq. (9) into Eq. (10) which results in

$$E\left[M\left(n\right)\right]=\frac{2+{\eta }^{-}}{6}\left(K\left(n\right)+1\right)$$E11

which is the expected number of candidate solutions that survive this elitist linear ranking selection process during the nth iteration. Ref. [17] shows that setting ${\eta }^{-}\approx 0.9$ often provides an adequate balance between selective pressure which allows for exploitation of the objective function and population diversity which allows for exploration of the objective function.

3.3 Differential evolution operator to improve convergence

Differential evolution algorithms generate new candidate solutions by adding a weighted difference between two randomly selected candidate solutions to a third randomly selected candidate solution. For this algorithm, the differential evolution operator to improve convergence generates a new candidate solution, $\mathbf{v}$, using

$$\mathbf{v}={\mathbf{x}}_k\left(n\right)+R\left[{\mathbf{x}}_m\left(n\right)-{\mathbf{x}}_j\left(n\right)\right]$$E12

where ${\mathbf{x}}_k\left(n\right)$ is the candidate solution randomly selected for differential evolution, ${\mathbf{x}}_m\left(n\right)$ and ${\mathbf{x}}_j\left(n\right)$ are two randomly selected candidate solutions from the mating pool, and $R$ is a uniformly distributed random number from the interval [0,1]. The two candidate solutions, ${\mathbf{x}}_m\left(n\right)$ and ${\mathbf{x}}_j\left(n\right)$, should be distinct and chosen so that ${\mathbf{x}}_m\left(n\right)\ne {\mathbf{x}}_j\left(n\right)$; however, this can become difficult when the algorithm is converging.

Because this algorithm is an elitist algorithm, the best candidate solution, ${\mathbf{x}}_{K\left(n\right)}\left(n\right)$, is always selected for this differential evolution operator. The other candidate solutions are selected randomly for differential evolution with a probability of $P_{\mathit{DE}1}\left(n\right)$. Because $R$ in Eq. (12) is a uniformly distributed random number from the interval [0,1], the value $R$ attenuates the difference between the two randomly selected candidate solutions. When this attenuated difference is added to the candidate solution selected for differential evolution, it creates a new candidate solution within a neighborhood of the candidate solution selected for differential evolution. As a result, this differential evolution operator improves the algorithm’s ability converge to an optimal point in the neighborhood. Figure 1 shows a plot of the contours of a two dimensional cost function, $J$, and three candidate solutions selected for differential evolution. The figure illustrates how this differential evolution operator creates candidate solutions within a neighborhood of the candidate solution, ${\mathbf{x}}_k\left(n\right)$, selected for differential evolution.

Figure 1.

A plot showing an example of the differential evolution operator that improves convergence for a two dimensional cost function, $J$. The plot shows the contour lines of the cost function and the candidate solution vectors involved in the differential evolution operation.

On average, this operator creates

$$\left(M\left(n\right)-1\right)P_{\mathit{DE}1}\left(n\right)+1$$E13

new candidate solutions.

3.4 Differential evolution mutation operator to improve diversity

Because differential evolution algorithms generate new candidate solutions by adding a weighted difference between two randomly selected candidate solutions to a third randomly selected candidate solution, the differential evolution operator creates new candidate solution vectors that contain elements that are different from the candidate solutions that formed the new solution vector. As a result, the differential evolution operator is often referred to as a mutation operator whether the operator creates similarity or diversity. In this chapter, the differential evolution operator that increases diversity is referred to as the differential evolution mutation operator.

For this algorithm, the differential evolution mutation operator generates a new candidate solution, $\mathbf{v}$, using

$$\mathbf{v}={\mathbf{x}}_k\left(n\right)+\frac{1}{4}\mathit{diag}\left(\mathbf{R}\right)\mathit{diag}\left({\mathbf{X}}_s\right)\left[{\mathbf{x}}_m\left(n\right)-{\mathbf{x}}_j\left(n\right)\right]$$E14

where ${\mathbf{x}}_k\left(n\right)$ is the candidate solution randomly selected for differential evolution mutation, ${\mathbf{x}}_m\left(n\right)$ and ${\mathbf{x}}_j\left(n\right)$ are two randomly selected candidate solutions from the mating pool, $\mathbf{R}$ is a vector whose elements are uniformly distributed random number from the interval [0,1], and ${\mathbf{X}}_s$ is a vector containing the solution space’s size for each dimension in rectangular coordinates or the diameters of an elliptic solution space. Again, the randomly selected candidate solutions, ${\mathbf{x}}_m\left(n\right)$ and ${\mathbf{x}}_j\left(n\right)$, should be distinct and chosen so that ${\mathbf{x}}_m\left(n\right)\ne {\mathbf{x}}_j\left(n\right)$.

Because this algorithm is an elitist algorithm, the best candidate solution, ${\mathbf{x}}_{K\left(n\right)}\left(n\right)$, is always selected for this differential evolution mutation operator. The other candidate solutions are selected randomly for differential evolution with a probability of $P_{\mathit{DE}2}\left(n\right)$. Because the term, $\frac{1}{4}\mathit{diag}\left(\mathbf{R}\right)\mathit{diag}\left({\mathbf{X}}_s\right)$, in Eq. (14) is a diagonal matrix with uniformly distributed random numbers from the interval zero to $\frac{1}{4}$ the size of each dimension of the solution space, the term, $\frac{1}{4}\mathit{diag}\left(\mathbf{R}\right)\mathit{diag}\left({\mathbf{X}}_s\right)$, typically increases each dimension of the difference between the two randomly selected candidate solutions randomly. As the entire population begins to converge and the differences between any two randomly selected candidate solutions begins to decrease, the term, $\frac{1}{4}\mathit{diag}\left(\mathbf{R}\right)\mathit{diag}\left({\mathbf{X}}_s\right)$, increases these small differences and these increased differences are added to the candidate solutions selected for differential Evolution. Therefore, the new candidate solutions typically lie outside the neighborhood of the candidate solutions selected for differential evolution. As a result, this differential evolution operator improves the algorithm’s diversity until the entire population begins to converge within very small differences.

The mean number of mutant solutions created by this process is

$$\left(M\left(n\right)-1\right)P_{\mathit{DE}2}\left(n\right)+1.$$E15

3.5 A recombination operator to improve convergence and diversity

Taguchi crossover can greatly increase convergence rates [11, 18]. As a result, when the differential evolution operators discussed earlier are combined with Taguchi crossover, this algorithm can converge too quickly. To prevent this algorithm from converging too quickly into a local minima or maxima, a recombination operator that creates a pair of new candidates is added to this algorithm. To improve convergence, this recombination operator generates a new candidate solution, $\mathbf{v}$, by averaging the selected candidate solution, ${\mathbf{x}}_k\left(n\right)$, with another randomly selected candidate solution, ${\mathbf{x}}_m\left(n\right)$, from the mating pool so that

$$\mathbf{v}=\left[{\mathbf{x}}_k\left(n\right)+{\mathbf{x}}_m\left(n\right)\right]/2.$$E16

To improve diversity, this recombination operator generates the another candidate solution by circularly shifting the elements of the newly formed candidate solution, v, by a uniformly distributed integer and then randomly changing the signs of the elements. In Matlab, this new vector, $\mathbf{w}$, can be created by

$$w=\mathit{sign}\left(\mathit{\text{randi}}\left(2,N,1\right)-1.5\right).{}^{\ast }\mathit{\text{circshift}}\left(v,\mathit{\text{randi}}\left(N\right)\right).$$E17

Because this algorithm is an elitist algorithm, the best candidate solution, ${\mathbf{x}}_{K\left(n\right)}\left(n\right)$, is always selected for this recombination operator. The other candidate solutions are selected randomly with a probability of $P_{\mathit{cr}}\left(n\right)$. On average, this operator creates

$$\left(M\left(n\right)-1\right)2P_{\mathit{cr}}\left(n\right)+2$$E18

new candidate solutions.

3.6 Solution space

A candidate solution is considered infeasible if it does not lie within the solution space. If a new candidate solution is infeasible, that solution is made feasible by one of two methods. If a convergence operator, such as the differential evolution for convergence or the recombination operator for convergence, creates an infeasible candidate solution, the infeasible solution vector is moved to the nearest edge of the solution space by changing the vector’s elements that lie outside solution space to the nearest edge of the solution space. This method attempts to generate feasible solutions within the neighborhood of the original infeasible solution so that the intent of the convergence operator that created the infeasible solution is maintained.

If a diversity operator, such as the differential evolution mutation for diversity operator or the recombination operator for diversity, creates an infeasible solution, the infeasible solution vector is moved into the solution space by performing a spatially circular shift of the infeasible solution vector’s elements. For example, if an infeasible solution, v, is created by a diversity operator in a hyper-rectangular solution space, the infeasible solution vector is moved into the solution space using

$$\mathbf{v}=\mathit{mod}\left[\mathbf{v}-\left({\mathbf{X}}_c-0.5{\mathbf{X}}_s\right),{\mathbf{X}}_s)\right]+\left({\mathbf{X}}_c-0.5{\mathbf{X}}_s\right)$$E19

where ${\mathbf{X}}_c$ is a vector containing the center of the solution space in rectangular coordinates, and ${\mathbf{X}}_s$ is a vector containing the size of each dimension of the solution space in rectangular coordinates. Similarly, if an infeasible solution, v, is created by a diversity operator in an elliptical solution space centered at the origin, the infeasible solution vector, v, expressed in polar coordinates, ${\mathit{re}}^{\mathit{j\theta }}$, is moved into the solution space using

$$\begin{array}{l}{r}_k=\mathit{rem}\left(r_k,r_{kmax}\right)\\ {\theta }_k={\theta }_k+\pi \end{array}$$E20

where $\mathit{rem}$ is the remainder function, $r_k$ is a radius that places the candidate solution outside of the solution space, $r_{kmax}$ is the maximum value of $r_k$ that keeps the candidate solution inside the solution space and ${\theta }_k$ is the angle associated with the radius, $r_k$. This method attempts to generate feasible solutions away from the neighborhood of the original infeasible solution so that the intent of the diversity operator that created the infeasible solution is maintained.

3.7 Taguchi crossover

A crossover operator is a recombination operator that combines the elements from two or more parent candidate solutions to generate a new offspring candidate solution. Taguchi crossover generates new candidate solutions by intelligently selecting elements from the two or more parent solutions vectors [11]. Taguchi crossover is a simple design of experiments method that creates a near optimal candidate solution from the parent candidate solutions. Consequently, Taguchi crossover can greatly increase an algorithm’s rate of convergence [11, 18].

Before selecting candidate solutions for Taguchi crossover, all new candidate solutions created by the other operators are added to the mating pool. The mean number of candidate solutions in the mating pool at this stage can be obtained by summing Eq. (14), Eq. (15) and Eq. (18) which implies that the mean number of candidate solutions in the mating pool at this stage is

$$M\left(n\right)+\left(M\left(n\right)-1\right)P_3\left(n\right)+4$$E21

where $P_3\left(n\right)=P_{\mathit{DE}1}\left(n\right)+P_{\mathit{DE}2}\left(n\right)+2P_{\mathit{cr}}\left(n\right)$.

Because this is an elitist algorithm, the best candidate solution is always selected for Taguchi crossover. The other candidate solutions from the mating pool are selected randomly for Taguchi crossover with a probability of $P_{\mathit{Tc}}$. For two level Taguchi crossover, crossover involving two parent solutions, one other candidate solution is selected randomly from the mating pool. For three level Taguchi crossover, crossover involving three parent solutions, two other candidate solutions are selected randomly from the mating pool.

On average, the Taguchi crossover operator creates

$$\left[\left(M\left(n\right)-1\right)\left(1+P_3\left(n\right)\right)+4\right]P_{\mathit{Tc}}+1$$E22

new candidate solutions.

3.8 Managing population size

Because the selection operator, the differential evolution operators, the recombination operators and Taguchi crossover operator generate a random number of new candidate solutions, the population size and mating pool size vary each generation, or iteration of the algorithm. After the Taguchi crossover operator, the average number, $E\left[K\left(n\right)\right]$, of the candidate solutions can be calculated by adding Eq. (21) and Eq. (22) which results in

$$E\left[K\left(n\right)\right]=\left[\left(M\left(n\right)-1\right)\left(1+P_3\left(n\right)\right)+4\right]\left(1+P_{\mathit{Tc}}\right)+2.$$E23

To maintain the population’s size, $K\left(n\right)$, at the population’s initial size, $K\left(0\right)$, the probabilities of at least one of the operators must vary so that

$$E\left[K\left(n\right)\right]=K\left(0\right).$$E24

Substituting Eq. (23) into Eq. (24) and solving for $P_3\left(n\right)$,

$$P_3\left(n\right)=\frac{K\left(0\right)-6-4P_{\mathit{Tc}}}{\left(1+P_{\mathit{Tc}}\right)\left(M\left(n\right)-1\right)}-1$$E25

where $P_{\mathit{Tc}}$ is assumed to be fixed and $P_3\left(n\right)$ varies. In this algorithm, the requirement in Eq. (25) is met by fixing $P_{\mathit{DE}1}\left(n\right)$ and $P_{\mathit{DE}2}\left(n\right)$, and letting

$$P_{\mathit{cr}}\left(n\right)=\left[P_3\left(n\right)-P_{\mathit{DE}1}\left(n\right)-P_{\mathit{DE}2}\left(n\right)\right]/2.$$E26

4. Optimization test function results

To evaluate this algorithm’s ability to solve optimization problems, the algorithm was applied to 13 commonly used global numerical optimization test functions. Table 1 lists these 13 cost functions, $J_1\left(\mathbf{x}\right)$ through $J_{13}\left(\mathbf{x}\right)$, where $\mathbf{x}={\left[\begin{array}{cccc}{x}_1& x_2& \dots & x_N\end{array}\right]}^T$, and their solution spaces. These functions include a spherical, three hyper-ellipsoid, the sum of different powers, Rastrigin’s, Schwefel’s, Griewank’s, Rosenbrock’s valley, Styblinski-Tang, Ackley’s Path, Price-Rosenbrock, and Eggholder’s functions. The first 11 functions, $J_1\left(\mathbf{x}\right)$ through $J_{11}\left(\mathbf{x}\right)$, are multidimensional functions and are tested for two dimensions ($N=2$) and 35 dimensions ($N=35$). Functions $J_{12}\left(\mathbf{x}\right)$ and $J_{13}\left(\mathbf{x}\right)$ are two dimensional functions and were only tested for $N=2$. For all 13 functions, $J_{\mathit{min}}=0$ where $J_{\mathit{min}}$ is the global minimum value of the cost function. Figure 2 shows a plot of the two dimensional Schwefel function and Figure 3 shows a plot of the two dimensional Eggholder function.

Function	Solution space $x_k\in \left[•,•\right]$
$J_1\left(\mathbf{x}\right)={\sum }_{k=1}^N{x}_k^2$	$\left[-\mathrm{5.12,5.12}\right]$
$J_2\left(\mathbf{x}\right)={\sum }_{k=1}^N{\mathit{kx}}_k^2$	$\left[-\mathrm{5.12,5.12}\right]$
$J_3\left(\mathbf{x}\right)={\sum }_{k=1}^N{\left[k\left(x_k-5k\right)\right]}^2$	$\left[-\mathrm{500,500}\right]$
$J_4\left(\mathbf{x}\right)={\sum }_{m=1}^N{\sum }_{k=1}^m{x}_k^2$	$\left[-\mathrm{65.536,65.536}\right]$
$J_5\left(\mathbf{x}\right)={\sum }_{k=1}^N{\left|x_k\right|}^{k+1}$	$\left[-1,1\right]$
$J_6\left(\mathbf{x}\right)=10N+{\sum }_{k=1}^N\left[x_k^2-10cos\left(2\pi x_k\right)\right]$	$\left[-\mathrm{5.12,5.12}\right]$
$J_7\left(\mathbf{x}\right)=418.9828872724338N-{\sum }_{k=1}^N{x}_{k}sin\left(\sqrt{\left|x_k\right|}\right)$	$\left[-\mathrm{500,500}\right]$
$J_8\left(\mathbf{x}\right)=1+{\sum }_{k=1}^N\frac{x_k^2}{4000}-{\prod }_{k=1}^{N}cos\left[\frac{x_k}{\sqrt{k}}\right]$	$\left[-\mathrm{600,600}\right]$
$J_9\left(\mathbf{x}\right)={\sum }_{k=1}^{N-1}\left[100{\left(x_k^2-x_{k+1}\right)}^2+{\left(x_k-1\right)}^2\right]$	$\left[-5,10\right]$
$J_{10}\left(\mathbf{x}\right)=39.16616570377142N+{\sum }_{k=1}^N\left(x_k^4-16x_k^2+5x_k\right)$	$\left[-5,5\right]$
$J_{11}\left(\mathbf{x}\right)=20+e-20e^{-0.2\sqrt{\frac{1}{N}{\sum }_{k=1}^n{x}_k^2}}-e^{{\sum }_{k=1}^N\frac{cos\left(2\pi x_k\right)}{N}}$	$\left[\mathrm{32.768,32.768}\right]$
$J_{12}\left(\mathbf{x}\right)=100{\left(x_2-x_1^2\right)}^2+{\left[6.4{\left(x_2-0.5\right)}^2-x_1-0.6\right]}^2$	$\left[-5,5\right]$
$\begin{array}{l}{J}_{13}\left(\mathbf{x}\right)=959.640662720851-x_{1}sin\left(\sqrt{\left|x_1-x_2-47\right|}\right)\\ -\left(x_2+47\right)sin\left(\sqrt{\left|\frac{x_1}{2}+x_2+47\right|}\right)\end{array}$	$\left[-\mathrm{512,512}\right]$

Table 1.

Optimization test functions and their solution spaces.

Figure 2.

A plot of the two dimensional Schwefel’s function.

Figure 3.

A plot of the two dimensional Eggholder’s function.

For all functions and dimensions, the initial population, $K\left(0\right)$, was set to 50, ${\eta }^{-}=0.9$, $P_{\mathit{DE}1}\left(n\right)=0.16$, $P_{\mathit{DE}2}\left(n\right)=0.2$, $P_{\mathit{Tc}}\left(n\right)=0.22$ and $P_{\mathit{cr}}\left(n\right)$ was set each iteration according to Eq. (18). The algorithm was applied 1000 times to each function, and the algorithm was assumed to converge when a solution, ${\mathbf{x}}_{\mathit{opt}}$, was determined so that $J\left({\mathbf{x}}_{\mathit{opt}}\right)\le J_{\mathit{tol}}$ where $J_{\mathit{tol}}$ for each function is listed in Table 2. Table 2 lists the mean and standard deviation of the number of iterations that the algorithm required to converge for two level ($L=2$) and three level ($L=3$) Taguchi crossover for each function. Table 2 also lists the means and standard deviations of the number of cost function evaluations that the algorithm required to converge for two level ($L=2$) and three level ($L=3$) Taguchi crossover for each function.

Cost Fn	$J_{\mathit{tol}}$	N = 2; Mean Iter; St Dev Iter		N = 2; Mean $J$ Evals; St Dev Evals		N = 35; Mean Iter; St Dev Iter		N = 35; Mean $J$ Evals; St Dev Evals.
		$L=2$	$L=3$	$L=2$	$L=3$	$L=2$	$L=3$	$L=2$	$L=3$
$J_1\left(\mathbf{x}\right)$	${10}^{-300}$	43	30	3950	4250	130	105	58,000	90,300
		22	17	2000	2400	4	3	2600	4200
$J_2\left(\mathbf{x}\right)$	${10}^{-300}$	45	31	4150	4300	137	110	61,000	94,400
		22	17	2050	2400	4	3	2500	4300
$J_3\left(\mathbf{x}\right)$	${10}^{-8}$	44	39	4000	5450	1170	975	$5.3e5$	$8.3e5$
		15	20	1350	2800	1540	1650	$7.0e5$	$1.4e6$
$J_4\left(\mathbf{x}\right)$	${10}^{-300}$	68	71	6200	9950	6917	2990	$3.1e6$	$2.5e6$
		34	33	3150	4700	1520	765	$6.9e5$	$6.6e5$
$J_5\left(\mathbf{x}\right)$	${10}^{-300}$	42	33	3800	4700	922	343	$4.1e5$	$2.9e5$
		17	15	1600	2200	108	47	48,550	40,200
$J_6\left(\mathbf{x}\right)$	${10}^{-300}$	31	24	2800	3300	84	65	37,500	55,400
		12	10	1150	1450	4	3	2420	3450
$J_7\left(\mathbf{x}\right)$	${10}^{-10}$	40	33	3650	4500	155	100	$69350$	86,350
		6	5	550	650	12	7	6200	6600
$J_8\left(\mathbf{x}\right)$	${10}^{-300}$	39	35	3550	4900	90	77	40,300	65,700
		19	18	1750	2600	4	4	2500	4650
$J_9\left(\mathbf{x}\right)$	${10}^{-10}$	65	67	5950	9500	14,325	5776	$6.5e6$	$4.9e6$
		82	88	7500	12,500	6920	4804	$3.1e6$	$4.1e6$
$J_{10}\left(\mathbf{x}\right)$	${10}^{-11}$	31	26	2850	3650	105	70	47,000	60,000
		3	3	300	450	6	3	2650	3250
$J_{11}\left(\mathbf{x}\right)$	${10}^{-15}$	40	28	3700	3900	110	72	49,600	61,400
		22	16	2050	2250	17	16	7400	14,200
$J_{12}\left(\mathbf{x}\right)$	${10}^{-2}$	19	18	1750	2550	NA	NA	NA	NA
		24	19	2200	2700	NA	NA	NA	NA
$J_{13}\left(\mathbf{x}\right)$	${10}^{-11}$	869	881	$7.9e4$	$1.2e5$	NA	NA	NA	NA
		1068	1162	$9.7e4$	$1.6e5$	NA	NA	NA	NA

Table 2.

Results for the optimization of the 2-D and 35-D test functions where $K\left(0\right)=50$, ${\eta }^{-}=0.9$, $P_{\mathit{DE}1}=0.16$, $P_{\mathit{DE}2}=0.2$ and $P_{\mathit{Tc}}=0.22$ for two level ($L=2$) and three level ($L=3$) Taguchi crossover. Results are averages over 1000 runs.

The number of cost function evaluations that the algorithm required to converge can also be calculated as a function of the algorithm’s average population size, average mating pool size and operator probabilities. For example, 2-D functions require four cost function evaluations for two level Taguchi crossover and nine cost function evaluations for three level Taguchi crossover. Therefore, the average number of cost function evaluations per algorithm iteration is

$$\overline{K}+4+4P_{\mathit{Tc}}\left[\left(\overline{M}-1\right)\left(1+{\overline{P}}_3\right)+4\right]$$E27

for two level Taguchi crossover and

$$\overline{K}+9+9P_{\mathit{Tc}}\left[\left(\overline{M}-1\right)\left(1+{\overline{P}}_3\right)+4\right]$$E28

for three level Taguchi crossover where $\overline{K}$ is the average population size, $\overline{M}$ is the average mating pool size after selection and ${\overline{P}}_3$ is the average of $P_3\left(n\right)$.

Similarly, for the 35-D functions, two level Taguchi crossover requires 40 cost function evaluations, and three level Taguchi crossover requires 81 cost function evaluations. Therefore, the average number of cost function evaluations per algorithm iteration is

$$\overline{K}+40+40P_{\mathit{Tc}}\left[\left(\overline{M}-1\right)\left(1+{\overline{P}}_3\right)+4\right]$$E29

for two level Taguchi crossover and

$$\overline{K}+81+81P_{\mathit{Tc}}\left[\left(\overline{M}-1\right)\left(1+{\overline{P}}_3\right)+4\right]$$E30

for three level Taguchi crossover.

Although no algorithm can solve all types of optimization problems [8, 9], the data in Table 2 shows that the algorithm converged below the specified $J_{\mathit{min}}$ for 100% of the 1000 runs for all the test functions. The data in Table 2 also shows that this algorithm requires significantly more iterations to converge for Eggholder’s function, $J_{13}$, and for Rosenbrock’s valley, $J_9$ when $N=35$ which implies that the algorithm has not been optimized for all types of cost functions. The data also shows that although the three level Taguchi crossover algorithm typically converges using few iterations than the two level Taguchi crossover algorithm, the two level Taguchi crossover algorithm typically requires fewer cost function evaluations than the three level Taguchi crossover algorithm.

5. Summary and conclusions

This chapter presents a hybrid genetic, differential evolution algorithm that represents the set of parameters being optimized in a vector. The algorithm uses an elitist, ranking, random selection method to generate a mating pool. Candidate solutions from the mating pool are randomly selected for two differential evolution operators, and two recombination operators. The new candidate solutions generated by these operators are added to the mating pool. Candidate solutions from this expanded mating pool are selected randomly for Taguchi crossover.

To evaluate this algorithm’s ability to solve optimization problems, the algorithm was applied to 13 commonly used global numerical optimization test functions, including a spherical, three hyper-ellipsoid, the sum of different powers, Rastrigin’s, Schwefel’s, Griewank’s, Rosenbrock’s valley, Styblinski-Tang, Ackley’s Path, Price-Rosenbrock, and Eggholder’s functions. The algorithm was evaluated using two and three level Taguchi crossover. For both two and three level Taguchi crossover, the algorithm converged below the specified $J_{\mathit{min}}$ for 100% of the 1000 runs for all the test functions. Although the three level Taguchi crossover algorithm typically converged using fewer iterations than the two level Taguchi crossover algorithm, the two level Taguchi crossover algorithm typically required fewer cost function evaluations than the three level Taguchi crossover algorithm.

Although this algorithm required significantly more iterations to converge for Eggholder’s function and for 35-D Rosenbrock’s valley function, ref. [19] shows that this algorithm has been successfully used to design digital infinite impulse response (IIR) filters with arbitrary magnitude responses. As a result, it can be expected that the simple optimization algorithm described in this chapter can be used successfully for similar engineering optimization applications.