Research Article  Open Access
Estimation of Distribution Algorithm Using Correlation between Binary Elements: A New BinaryCode Metaheuristic
Abstract
A new metaheuristic called estimation of distribution algorithm using correlation between binary elements (EDACE) is proposed. The method searches for optima using a binary string to represent a design solution. A matrix for correlation between binary elements of a design solution is used to represent a binary population. Optimisation search is achieved by iteratively updating such a matrix. The performance assessment is conducted by comparing the new algorithm with existing binarycode metaheuristics including a genetic algorithm, a univariate marginal distribution algorithm, populationbased incremental learning, binary particle swarm optimisation, and binary simulated annealing by using the test problems of CEC2015 competition and one realworld application which is an optimal flight control problem. The comparative results show that the new algorithm is competitive with other established binarycode metaheuristics.
1. Introduction
Nowadays in the economiccompetitive world, optimisation has become increasingly popular for real applications as it is a powerful mathematical tool for solving a wide range of engineering design types. Once an optimisation problem is posed, one of the most important elements in the optimisation process is an optimisation method or an optimiser used to find the optimum solution. Optimisers can be categorised as the methods with and without using function derivatives. The former are traditionally called mathematical programming or gradientbased optimisers whereas the latter have various subcategories. One of them is a metaheuristic (MH). The term metaheuristics can cover natureinspired optimisers [1–10], swarm intelligent algorithms [11–20], and evolutionary algorithms [21–24]. Most of them are based on using a set of design solutions, often called a population, for searching an optimum. The main operator usually consists of the reproduction and selection stages. The advantages of such an optimiser are simplicity to use, global optimisation capability, and flexibility to apply as it is derivativefree. However, it still has a slow convergence rate and search consistency. These issues have made researchers and engineers around the globe investigate how to improve the search performance of MHs.
A genetic algorithm (GA) [21] is probably the best known MH while other popular methods are differential evolution (DE) [22] and particle swarm optimisation (PSO) [17]. Among MH algorithms, they can be categorised as the methods using real, binary, or integer codes. The mix of those types of design variables and some other types can also be made. This makes MHs considerably appealing for use with realworld applications particularly for those design problems that function derivatives are not available or impossible to calculate. Most MHs are based on continuous design variables or real codes. For single objective optimisation, there have been numerous realcode MHs being developed. At the early stage, methods like evolutionary programming [25] and evolution strategies [26] were proposed. Then, DE and PSO were introduced. Until recently, there have been probably over a hundred new realcode MHs in the literature. Some recent algorithms include, for example, a sinecosine algorithm [27], a grey wolf optimiser [20], teachinglearningbased optimisation [2], and Jaya algorithm [28]. Meanwhile, powerful existing algorithms such as PSO and DE have been upgraded by integrating into them some types of selfadaptive schemes, for example, adaptive differential evolution with optional external archive (JADE) [29], SuccessHistory Based Parameter Adaptation for Differential Evolution (SHADE) [30], SHADE Using Linear Population Size Reduction (LSHADE) [31], and adaptive PSO [32–34]. MHs are even more popular when they can be used to find a Pareto front of a multiobjective optimisation problem within one optimisation run. Such a type of algorithm is usually called multiobjective evolutionary algorithms (MOEAs) where some of the best known algorithms are nondominated sorting genetic algorithm (NSGAI, NSGAII, and NSGAIII) [35–37], multiobjective particle swarm optimisation [38], strength Pareto evolutionary algorithm [39], multiobjective grey wolf optimisation [40], multiobjective teachinglearningbased optimisation [41], multiobjective evolutionary algorithm based on decomposition [42], multiobjective ant colony optimisation [43], multiobjective differential evolution [44], and so forth. One of the most challenging issues in MHs is to improve their ability for tackling manyobjective optimisation (a problem with more than three objectives). Some recently proposed algorithms are knee pointdriven evolutionary algorithm [45], an improved twoarchive algorithm [46], preferenceinspired coevolutionary algorithms [47], and so forth.
In practice, GA a metaheuristic using binary strings is arguably the most used method as it is included in engineering software such as MATLAB. Apart from GA, other MHs using a binary string representing a design solution include a univariate marginal distribution algorithm (UMDA) [48], populationbased incremental learning (PBIL) [24], binary particle swarm optimisation (BPSO) [49], binary simulated annealing (BSA) [50], binary artificial bee colony algorithm based on genetic operator (GBABC) [51], binary quantuminspired gravitational search algorithm (BQIGSA) [52], and selfadaptive binary variant of a differential evolution algorithm (SabDE) [53]. With the popularity of GA, a binarycode MH has been rarely developed and proposed while its realcode counterparts have over a hundred different search concepts reported in the literature. That means there are possible more than a thousand realcode MH algorithms being published. It should be noted that realcode MHs can be modified to solve binarycode optimisation by means of binarisation [54].
This paper is therefore devoted to the further development of a binarycode metaheuristic. The method is called estimation of distribution algorithm using correlation between binary elements (EDACE). Performance assessment is made by comparing the proposed optimiser with GA, UMDA, BPSO, BSA, and PBIL by using the CEC2015 test problems. Also, the realworld optimal flight control is used for the assessment. The comparative results are obtained and discussed. It is shown that EDACE is among the top performers.
2. Proposed Method
The simplest but efficient estimation of distribution algorithm is probably populationbased incremental learning (PBIL). Another MH that uses a similar concept is UMDA. Unlike GA which uses a matrix containing the whole binary solutions during the search, PBIL uses the socalled probability vector to represent a binary population. During an optimisation process, the probability vector is updated iteratively until approaching an optimum. In EDACE, a matrix called a correlation between binary elements (CBE) matrix is used to represent a binary population. The matrix can be denoted as , where the value of the element indicates the correlation between element and element of a binary design solution. The higher value of means the higher probability that binary elements and will have the same value. The algorithm is developed to deal with a boxconstrained optimisation problem:where is an objective function and is a vector containing design variables (a design vector). and are the lower and upper bounds of , respectively. Assuming that a design vector can be represented by a row vector of binary bits size , the CBE matrix thus has the size of . It should be noted that the details of converting a binary string to be a design vector can be found in [55]. In generating a binary string from the CBE matrix, a reference binary solution (RBS) is needed. It can be a randomly generated solution or the best solution found so far depending on a user preference. Then, a row of the matrix is randomly selected (say the th row). The th element of a generated binary solution is set to be the th element of the reference binary solution. The rest of the created binary elements are based on the value of ; . The procedure for creating a binary solution sized from the CBE matrix is detailed in Algorithm 1 where is a binary design solution, is the reference binary solution, is a population size, and is a uniform random number. The algorithm spends loops for creating binary solutions. The process for generating a binary solution from the CBE matrix is in steps –. For one binary solution, only one randomly selected row of CBE (say row ) is used (step ). Then, the th element of a generated binary solution is set equal to the th element of the reference binary solution, . The rest of the elements of the generated binary solution are created in such a way that their values depend on corresponding elements on the th row of CBE. From the computation steps –, the value of determines the probability of to be the same as . The higher value of means the higher correlation between elements and and consequently the higher probability that will be set equal to .

The CBE matrix is a square symmetric matrix with equal size to the length of a binary solution whose all diagonal elements are equal to one. For an iteration, the matrix will be updated according to the so far best solution (). The learning rate () with be used to control the changes in updating as with PBIL. Once is updated, the value of is set to be which means the process requires updates since is always set to be 1. The updated denoted by can be calculated fromwhere is the learning rate randomly generated in the interval . and are the th and th elements of , respectively. From the updating equation, if the th and th elements are similar, it means they are correlated; consequently, the value of (and ) is increased. If they are dissimilar or uncorrelated, is then decreased. Nevertheless, the value of must be limited to the predefined interval where and are the predefined lower and upper limits of . Equation (3) is used to maintain diversity in optimisation search. In the original PBIL, a mutation operator is used with the same purpose. Therefore, the procedure of EDACE starts with an initial matrix for correlation between binary elements where and . This implies that when generating a binary solution, its elements have equal probability to be 1 or 0 where its th element can be 1 or 0, created at random. The procedure for general purpose of EDACE is given in Algorithm 2. The decision on selecting for generating a binary solution and for updating the CBE matrix is dependent on a preference of a user. This means other versions of EDACE can be developed in the future.

An initial binary population is randomly created. The binary solutions are then decoded to be real design variables where function evaluations are performed and and are found. Then, new binary solutions are generated using Algorithm 1 while the greedy selection (steps –) is activated with and being determined. The CBE matrix is updated by using as detailed in (2)(3). The search process is repeated until termination criterion is reached. The generation of a binary design solution of EDACE is, to some extent, similar to those used in binary PSO [49] and binary quantuminspired gravitational search algorithm (BQIGSA) [52] in the sense that the binary solution is controlled by the probability of being “1” or “0”. However, in EDACE, a generated solution relies not only on such probability but also on the reference binary solution . Apart from that, the update of CBE tends to be similar to the concept employed in PBIL with a learning rate and this is totally different from binary PSO and BQIGSA.
In selecting and , if both solutions are the same which is , it could lead to a premature convergence. If both are set to be a solution randomly selected solution from the current binary population, the diversification increases but the convergence rate will be slower. Therefore, the balance between intensification and diversification must be made. In this work, the so far best binary solution is set to be to maintain intensification. For updating the CBE matrix, we use the new updating scheme asThe solutions and are two types of best solutions. Firstly, best solutions are selected from (see Algorithm 2 for both solution sets), sorted according to their functions, and then saved to a set Best_sol. Four vectors are created as the so far best solution, a solution whose elements are averaged from the elements of the first (default = 10) best solutions found so far, a solution whose elements are averaged from the elements of the members of Best_sol, and a solution whose elements are averaged from the elements of the current binary population. is randomly chosen from the aforementioned solutions (, , , and ) with equal probability while is randomly chosen from the members of Best_sol. With this idea, the balance between exploration and exploitation is maintained throughout the search process. Algorithm 3 shows the new CBE updating strategy.

3. Experimental SetUp
To investigate the search performance of the proposed algorithm, fifteen learningbased test problems from CEC2015 and one flight dynamic control optimisation problem are used. The former are used for testing the performance of EDACE for general types of boxconstrained optimisation while the latter is the realworld application.
3.1. CEC2015 LearningBased Test Problems
The CEC2015 learningbased test problems are boxconstrained single objective benchmark functions proposed in [56]. The problems consist of 2 Unimodal Functions, 3 Simple Multimodal Functions, 3 Hybrid Functions, and 7 Composition Functions. The summary of CEC2015 learningbased test problems is shown in Table 1. It should be noted that the details and the codes for the test problems can be downloaded from the website of CEC2015 competition.

3.2. Flight Dynamic Control Optimisation Problem
Flight dynamic control system design is a classical important application for real engineering problems. The motion of an aircraft can be described using the body axes which is herein the stability axes consisting of roll axis (), pitch axis (), and yaw axis () as shown in Figure 1. The motion of the aircraft is described by Newton’s 2nd law or equations of motion for both translational and rotational motions. The dynamical model is nonlinear but can be linearised by applying aerodynamic derivatives. Due to aircraft symmetry with respect to the plane, the linearised dynamical model can be decoupled into two groups as longitudinal motion and the lateral/directional motion. For more details of deriving the equations of motion, see [57]. In this work, only the lateral/directional motion control is considered. A state equation representing the dynamic motion of an aircraft is expressed as follows [57–60]:where , is the sideslip, a velocity in direction, is the yaw rate, rate of change of rotation about the axis, is the roll rate, rate of change of rotation about the axis, is the bank angle, rotation about the axis, is the kinetic energy matrix, is Coriolis matrix, is the control vector, is the aileron deflection, and is the rudder deflection.
The control vector can be expressed aswhere is a pilot’s control input vector while and are the gain matrices expressed as follows [59]:where parameters are control gain coefficients which need to be found.
From (5)(6), the state equation for lateral/directional motion of an aircraft can be expressed as
Design optimisation of the control system of an aircraft is found to have many objectives as there are several criteria that need to be satisfied such as control stability, accuracy, sensitivity, and control effort, while the control gains coefficients are set to be design variables for an optimisation problem. In this work, the optimal flight control of an aircraft focuses on only the stability aspect. The objective function is posed to minimise spiral root subjected to stability performance constraints. The optimisation problem can then be written aswhere , , , and are spiral root, roll damping, damping ratio of Dutchroll complex pair, and Dutchroll frequency, respectively. These parameters can be calculated based on the eigenvalues associated with the matrix . The design variables are control gain coefficients in the matrix . The kinetic energy matrix () and the Coriolis matrix () are defined as
More details about this aircraft dynamic model can be found in [58–60]. To handle the constraints, the penalty function which was presented in [61] is used.
The proposed EDACE and several wellestablished binarycode metaheuristics are used to solve the fifteen CEC2015 learningbased test problems and the flight dynamic control test problem. The metaheuristic optimisers are as follows: Genetic algorithm (GA) [21] used binary codes with crossover and mutation rates are 1 and 0.1, respectively. Binary simulated annealing (BSA) [50] used binary codes with exponentially decreasing temperature. The starting and ending temperature are set to be 10 and 0.001, respectively. The cooling step is set as 10. Populationbased incremental learning (PBIL) [24] used binary codes with the learning rate, mutation shift, and mutation rate as 0.5, 0.7, and 0.2, respectively. Binary particle swarm optimisation (BPSO) [49] used binary codes with Vshaped transfer function while the transfer function used is the Vshaped version 4 (V4) as reported in [49]. It is noted that this version is said to be the most efficient version based on the results obtained in [49]. Univariate marginal distribution algorithm (UMDA) [48] used binary codes. The first 20 best binary solutions are used to update the probability matrix. Estimation of distribution algorithm with correlation of binary elements (EDACE) (Algorithm 2) used binary codes with , , , , and .
Each algorithm is used to solve the problems for 30 optimisation runs. The population sizes are set to be 100 and 20 while number of generation is set to be 100 and 500 for the CEC2015 learningbased test problems and the flight dynamic control test problem, respectively. For an algorithm using different population size and number of generations such as BSA, it will be terminated at the same number function evaluations, which is 10,000 for all test problems. The binary length is set to be 5 for each design variable for all optimisers.
4. Optimum Results
4.1. CEC2015
After applying the proposed EDACE and several wellestablished binary MHs for solving the CEC2015 learningbased benchmark functions, the results are shown in Tables 2–4. Note that, apart from the algorithms used in this study, the results of solving CEC2015 test suit obtained from efficient binary artificial bee colony algorithm based on genetic operator (GBABC), binary quantuminspired gravitational search algorithm (BQIGSA), and selfadaptive binary variant of a differential evolution algorithm (SabDE) as reported in [53] are also included in the comparison. From Table 2, the mean (Mean) and standard deviation (STD) values of the objective functions are used to measure the search convergence and consistency of the algorithms. The lower Mean is the better convergence while the lower STD is the better consistency. The value of Mean is more important; thus, for method A with lower Mean but higher STD than method B, method A is considered to be superior.
