Estimation of Distribution Algorithm Using Correlation between Binary Elements : A New Binary-Code Metaheuristic

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 binary-code metaheuristics including a genetic algorithm, a univariatemarginal distribution algorithm, population-based incremental learning, binary particle swarm optimisation, and binary simulated annealing by using the test problems of CEC2015 competition and one real-world application which is an optimal flight control problem. The comparative results show that the new algorithm is competitive with other established binary-code metaheuristics.


Introduction
Nowadays in the economic-competitive 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 gradient-based optimisers whereas the latter have various subcategories.One of them is a metaheuristic (MH).The term metaheuristics can cover nature-inspired optimisers [1][2][3][4][5][6][7][8][9][10], swarm intelligent algorithms [11][12][13][14][15][16][17][18][19][20], and evolutionary algorithms [21][22][23][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 derivative-free.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 real-world 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 real-code 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 real-code MHs in the literature.Some recent algorithms include, for example, a sine-cosine algorithm [27], a grey wolf optimiser [20], teaching-learning-based 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 self-adaptive schemes, for example, adaptive differential evolution with optional external archive (JADE) [29], Success-History Based Parameter Adaptation for Differential Evolution (SHADE) [30], SHADE Using Linear Population Size Reduction (LSHADE) [31], and adaptive PSO [32][33][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 (NSGA-I, NSGA-II, and NSGA-III) [35][36][37], multiobjective particle swarm optimisation [38], strength Pareto evolutionary algorithm [39], multiobjective grey wolf optimisation [40], multiobjective teaching-learning-based 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 many-objective optimisation (a problem with more than three objectives).Some recently proposed algorithms are knee point-driven evolutionary algorithm [45], an improved two-archive algorithm [46], preference-inspired 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], population-based 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 self-adaptive 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 real-code counterparts have over a hundred different search concepts reported in the literature.That means there are possible more than a thousand real-code MH algorithms being published.It should be noted that real-code MHs can be modified to solve binary-code optimisation by means of binarisation [54].
This paper is therefore devoted to the further development of a binary-code 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 real-world 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.

Proposed Method
The simplest but efficient estimation of distribution algorithm is probably population-based 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 so-called 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   ∈ [0, 1], 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 box-constrained optimisation problem: where  is an objective function and x is a vector containing design variables (a design vector).x  and x  are the lower and upper bounds of x, respectively.Assuming that a design vector can be represented by a row vector of binary bits size  × 1, 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  × 1 from the  ×  CBE matrix is detailed in Algorithm 1 where b is a binary design solution, b REF is the reference binary solution,   is a population size, and rand ∈ [0, 1] 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 (3)- (12).For one binary solution, only one randomly selected row of CBE (say row ) is used (step (4)).Then, the th element of a generated binary solution is set equal to the th element of the reference binary solution, b REF .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 (5)- (11), 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 (b best ).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 ( − 1)/2 Set a = { } a vector used to contain elements of a generated binary string.(4) Randomly select a position (th row) of P. (5) Set   =  REF, .% Set the th element of a as the th element of b REF .(6) For  = {1, 2, . . ., } − {} (7) If rand <   (8)   =   %   and   values are equal, which are either "0" or "1".
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   = 1 and   = 0.5.This implies that whengenerating 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 b REF for generating a binary solution and b best 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 b REF and b best are found.Then, new binary solutions are generated using Algorithm 1 while the greedy selection (steps (6)-( 8)) is activated with b REF and b best being determined.The CBE matrix is updated by using b best 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 quantum-inspired 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 b REF .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 b REF and b best , if both solutions are the same which is b best , 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 b REF to maintain intensification.For updating the CBE matrix, we use the new updating scheme as The solutions b best1 and b best2 are two types of best solutions.Firstly,   best solutions are selected from {b  } ∪ {b  new } (see Algorithm 2 for both solution sets), sorted according to their functions, and then saved to a set Best sol.For  = 1 to   (5) Calculate objective function values   new = fun(x  new ).( 6) balance between exploration and exploitation is maintained throughout the search process.Algorithm 3 shows the new CBE updating strategy.

Experimental Set-Up
To investigate the search performance of the proposed algorithm, fifteen learning-based 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 box-constrained optimisation while the latter is the real-world application.

CEC2015 Learning-Based Test
Problems.The CEC2015 learning-based test problems are box-constrained 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 learning-based 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.

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][58][59][60]:  The control vector u can be expressed as where u  is a pilot's control input vector while C and K are the gain matrices expressed as follows [59]: where parameters  1 - 7 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 More details about this aircraft dynamic model can be found in [58][59][60].To handle the constraints, the penalty function which was presented in [61] is used.
The proposed EDACE and several well-established binary-code metaheuristics are used to solve the fifteen CEC2015 learning-based 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.Population-based 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 V-shaped transfer function while the transfer function used is the V-shaped 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   = 0.1,   = 0.9,  , = 0.4,  , = 0.6, and  best = 10.
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 learning-based 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.

Optimum Results
4.1.CEC2015.After applying the proposed EDACE and several well-established binary MHs for solving the CEC2015 learning-based 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 quantum-inspired gravitational search algorithm (BQIGSA), and self-adaptive 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.
For the measure of search convergence based on the mean objective function values, the best performer for the unimodal test functions, 1 and 2, is EDACE while the second best is BPSO.For the simple multimodal functions, the best performer for 4 and 5 is SabDE while the best performer for 3 is BPSO.The second best performers for 3, 4, and 5 are SabDE, BEDACE, and UMDE, respectively.For the hybrid functions, the best performers for the functions 6, 7, and 8, are SabDE, EDACE, and BPSO, respectively, while the second best performer for 6 and 7 is BPSO and the second best for 8 is EDACE.For the final group of CEC2015 test problems, composition functions, the best performer for the 11, 12, and 14 is SabDE while the best performers for the 10 and 15 are BPSO and EDACE, respectively.For 9, the best performers are UMDA, BPSO, GA, PBIL, and EDACE, which obtain the same mean values while, for 13, the best performers are UMDA, BPSO, GA, PBIL, BSA, and EDACE, which obtained the same mean values.It should be noted that the results from [53] were obtained from using the total number of function evaluations as 1,000,000 with the binary length of 50 for each design variable whereas this work uses 10,000 function evaluations with the binary length of 5 for each design variable.This indirect comparison with GBABC, BQIGSA, and SabDE can only be used to show that the proposed EDACE also has good performance and cannot be used to claim which method is superior.
For the measure of search consistency based on the STD values, the most consistent methods for unimodal functions, 1 and 2, are BPSO and EDACE while the second most consistent methods are EDACE and BPSO, respectively.For the simple multimodal functions, the best for 3 and 5 is SabDE while the best for 4 is the proposed EDACE.EDACE is the best for the hybrid function of 7 while BPSO is the best for the hybrid functions 6 and 8.For the composition functions, EDACE is the best for the problems 9 and 12 while BPSO is the best for 10.For the composition functions, 11, 14, and 15, the best is SabDE while the best for 13 is BSA.

Results
reported in [53] with 1,000,000 function evaluations and 50 binary lengths for each design variable.
The value Min in Table 2 is the objective function value of the best run from a particular method.Note that only the UMDA, BPSO, GA, PBIL, BSA, and EDACE were compared.For the unimodal function, the minimum objective function values of 1 and 2 were obtained by BPSO and EDACE, respectively.For the simple multimodal functions, the minimum objective function values for 3 and 5 are obtained from BPSO and EDACE, respectively, while for 4, the minimum is obtained from UMDA, BSA, and EDACE.The EDACE obtained minimum objective function values for all test functions in the hybrid function group.However, for the hybrid function 8, three algorithms including BPSO, GA, and EDACE obtained the minimum values.For the composition functions, EDACE obtained the minimum function values for all test functions.However, for the functions 9 and 13, all algorithms obtained the same minimum values while for the 11, BPSO and EDACE obtained the same minimum function values.Similarly, for 12, UMDA, BPSO, BSA, and EDACE obtained the same minimum values.
Table 3 shows the summary of ranking based on the mean objective function values from 30 optimisation runs.It was found that the proposed EDACE is mostly ranked in top three best from solving fifteen CEC2015 learning-based test problems.After summing up the ranking score, it is found that EDACE and BPSO are equal best performer while the third best is UMDA.
In order to further investigate the performance comparison of the binary-code MHs, the statistical -test is employed.Table 4 shows a 9 × 9 comparison matrix of the 9 optimisers.If method  is significantly better than method  based on the -test at 5% significant level, the column  and row  of the matrix are set to be 1; otherwise, they are set to be 0. When  summing up along the columns, the highest score indicates the best optimiser based on this type of comparison.In the table, it means EDACE is the best.Table 5 shows the ranking of the 9 optimisers when solving all CEC2015 learning-based test problems based on the -test.After summing up the ranking numbers of all test problems, it is found that EDACE is the overall best optimiser while BPSO and UMDA are the second and the third best, respectively.
Figures 2-5 show the search history of the top three optimisers EDACE, BPSO, and UMDA on solving all CEC2015 learning-based test problems where the vertical axis is the average objective function from 30 runs of each method.For all test functions, it was found that EDACE and UMDA converged to the optimal values at higher speed while BPSO seems to converge slowly and consistently.However, for all functions, BPSO finally moves to the minimum or near  Table 6 shows performance of EDACE on solving unimodal function, 1, when the binary lengths for each design variable are 5, 10, 25, and 50 for 10 optimisation runs.It was found that when the number of binary bits increases, the computational time increases and the resulting mean objective function values decrease for the binary lengths less than 25.However, for the binary length of 50, the mean objective function value increases meaning EDACE performance deteriorates.Without considering computational time, the best number of binary length is 25.

Flight Dynamic Control System
Design.After applying the six binary-code MHs to solve the real engineering application of flight dynamic and control system for 30 optimisation runs, the comparison results are shown as box-plots of the objective and constraint violation values (Figure 6).The upper and lower horizontal lines of each box represent the maximum and minimum of objective function values, respectively, while the internal line shows the median of objective function values.From this figure, based on median values of objective function, it is found that the best performer is EDACE while the second best and the third best are BPSO and UMDA, respectively.The most consistent method having the smallest gap between the maximum and minimum for all of optimisation runs is UMDA.However, the worst function value found by EDACE is almost as good as the best found by UMDA.Thus, the proposed EDACE is superior.Based on the figure, it was found that GA failed to solve the problem as it cannot obtain a feasible optimum point.The minimum objective function value is obtained from using the proposed EDACE.
Figure 7 shows the best run search history of all optimisers (selection based on the minimum objective function values of feasible solutions).From the figure, UMDA and PBIL seem to be the fastest convergent methods initially.However, after the process goes on for about 4,000 function evaluations, the proposed EDACE converged to the minimum objective function value with a faster rate than the others.It has better exploration rate as the best function value is still decreased at the late iteration numbers.BPSO, on the other hand, seems to be slower than UMDA, PBIL, and BSA in the beginning.It however can converge to the better results after around 8,000 function evaluations.studies.The use of EDACE for hyperheuristic development is also possible.The extension to multiobjective optimisation and many-objective optimisation is also under investigation.Appling EDACE for the more complex problems such as large scale problems, mixed-variable problems, and reliability optimisation is for future work.The fight control optimisation problem, one of our recent research focuses, has more than three objective functions to be optimised; thus, it should be formulated as many-objective optimisation.This along with aircraft path planning dynamic optimisation still needs considerably more investigation while EDACE will be one of optimisers to be used for solving such design problems.

1 ) 3 )
Four  × 1 vectors are created as b 1 the so far best solution, b 2 a solution whose elements are averaged from the elements of the first  best (default = 10) best solutions found so far, b 3 a solution whose elements are averaged from the elements of the members of Best sol, and b 4 a solution whose elements are averaged from the elements of the current binary population.b best1 is randomly chosen from the aforementioned solutions (b 1 , b 2 , b 3 , and b 4 ) with equal probability while b best2 is randomly chosen from the members of Best sol.With this idea, the Input: number of generation ( iter ), population size (  ), binary length () Output: b best ,  best Initialisation: (0.1) Assign   = 0.5 and   = 1, sized  × .(0.2) Randomly generate   binary solutions b  and decode them to be x  .(0.3) Calculate objective function values   = fun(x  ) where fun is an objective function evaluation.(0.4) Find  best , b REF , b best Main iterations (For iter = 1 to  iter (2)Update P using Equation (2) (Generate b  new from P using Algorithm 1, and decode them to be x  new .(4) 7 ( = 10) 1 5 0 0 where x = {, , , }  ,  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, A is the kinetic energy matrix, B is Coriolis matrix, u = {    } is the control vector,   is the aileron deflection, and   is the rudder deflection.

Figure 1 :
Figure 1: Stability axes of an aircraft.

Figure 2 :
Figure 2: Search history of the top three best optimisers based on the -test for the unimodal function.

Figure 3 :
Figure 3: Search history of the top three best optimisers based on the -test for the simple multimodal functions.

Figure 4 :
Figure 4: Search history of the top three best optimisers based on the -test for the hybrid functions.
In this work, a new concept of a binary-code optimiser is proposed.Fifteen CEC2015 learning-based test problems and a real engineering design problem of flight dynamic and control system are used to investigate the search performance of the proposed algorithm.Several well-established binarycode MHs are used in comparison.The results obtained show that the proposed EDACE is the best performer on solving the 15 CEC2015 learning-based test problems and real engineering design problem of flight dynamic and control.Further improvement of EDACE by means of self-adaptation will be investigated in the future.The choice for b REF needs further

6 Figure 5 :
Figure 5: Search history of the top three best optimisers based on the -test for the composition functions.
, ,  , ].  best, and  best, are the th and th elements of b best , 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

Table 1 :
Summary of CEC2015 learning-based functions.

Table 3 :
Ranking of all optimisers based on the Mean values.

Table 4 :
Comparison based on the statistical -test of the test problem.

Table 5 :
Ranking of all optimisers for all CEC2015 learning-based test problem based on statistical -test.

Table 6 :
The table shows performance of EDACE for various number of binary bits.