*Article* **An Improved Butterfly Optimization Algorithm for Engineering Design Problems Using the Cross-Entropy Method**

#### **Guocheng Li 1,***∗***, Fei Shuang 2, Pan Zhao <sup>1</sup> and Chengyi Le 3,***<sup>∗</sup>*


Received: 4 July 2019; Accepted: 9 August 2019; Published: 14 August 2019

**Abstract:** Engineering design optimization in real life is a challenging global optimization problem, and many meta-heuristic algorithms have been proposed to obtain the global best solutions. An excellent meta-heuristic algorithm has two symmetric search capabilities: local search and global search. In this paper, an improved Butterfly Optimization Algorithm (BOA) is developed by embedding the cross-entropy (CE) method into the original BOA. Based on a co-evolution technique, this new method achieves a proper balance between exploration and exploitation to enhance its global search capability, and effectively avoid it falling into a local optimum. The performance of the proposed approach was evaluated on 19 well-known benchmark test functions and three classical engineering design problems. The results of the test functions show that the proposed algorithm can provide very competitive results in terms of improved exploration, local optima avoidance, exploitation, and convergence rate. The results of the engineering problems prove that the new approach is applicable to challenging problems with constrained and unknown search spaces.

**Keywords:** global optimization; meta-heuristic; butterfly optimization algorithm; cross-entropy method; engineering design problems

#### **1. Introduction**

Real-world engineering design optimization problems are very challenging to find the global optimum of a highly complex and multiextremal objective function, involving many different decision variables under complex constraints [1,2]. A general engineering design optimization problem can be stated as follows:

$$\min f(\mathbf{x}), \mathbf{x} = (\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n)^T \in R^n,\tag{1}$$

subject to

$$\mathcal{S}\_j(\mathbf{x}) \le 0, j = 1, 2, \dots, p,\tag{2}$$

$$h\_k(\mathbf{x}) = 0, k = 1, 2, \cdots, q,\tag{3}$$

$$l!b\_i \le x\_i \le u b\_{i\prime} i = 1, 2, \cdots, n,\tag{4}$$

where *n* is the number of search space dimensions, *gj*(*x*) and *hk*(*x*)are the *j*th inequality constraint and *k*th equality constraint, respectively. *lbi* and *ubi* represent the lower and upper bound of the value of *xi*.

Most of the constraints of the global optimization problem are nonlinear. Such nonlinearity often results in a multimodal response landscape [3]. Due to their high complexity, nonlinearity, and multimodality, traditional optimization methods, such as gradient-based optimization algorithms, are no longer suitable for this problem [4]. Many researchers have developed various derivative-free global optimization methods for it. In general, these methods can be divided into two classes: deterministic methods and stochastic meta-heuristic algorithms [4,5]. The former, such as the Hill-Climbing [6], Newton–Raphson [7], and DIRECT algorithm [8], can obtain the same optimal results based on the same set of initial values. However, this behavior results in local optima entrapment and a loss of reliability in finding the global optimum, since most real engineering design problems have extremely high numbers of local solutions [9]. The latter, such as the Genetic Algorithm (GA) [10,11], Particle Swarm Optimization (PSO) [12], Cuckoo Search (CS) [13], Bat Algorithm (BA) [14], Grey Wolf Optimizer (GWO) [15], Forest Optimization Algorithm (FOA) [16,17], Ant Lion Optimizer (ALO) [9], Whale Optimization Algorithm (WOA) [18], Crow Search Algorithm (CSA) [19], Salp Swarm Algorithm (SSA) [20], and Butterfly Optimization Algorithm (BOA) [21,22], mostly benefit from their simplicity, flexibility, and stochastic operators, which makes them different from deterministic methods [4], and they have become very popular optimization techniques for such optimization problems in recent years.

Despite the merits of nature-inspired meta-heuristic algorithms, it has been pointed out that most of them are unable to guarantee global convergence [23]. At the same time, the No Free Lunch (NFL) theorem has proved that none of these methods are able to solve all global optimization problems [24]. They perform very well when dealing with certain optimization problems, but fail in most cases. An excellent meta-heuristic algorithm has two symmetric search capabilities: local search and global search, and achieves a proper balance between exploration and exploitation to effectively avoid it falling into a local, and find the global optimum. In order to enhance global search capabilities, many hybrid meta-heuristic algorithms have been developed by combining meta-heuristics with exact algorithms or other meta-heuristics for solving more complicated optimization problems [25–29].

Inspired by the idea of global random search algorithms [5,30], and based on co-evolution, this paper explores an improved butterfly optimization algorithm (BOA) using the cross-entropy method, which is a global stochastic optimization method based on a statistical model. The original BOA was proposed by Arora and Singh [21] in 2015, and was inspired by the food foraging behavior of butterflies, while the cross-entropy (CE) method was developed by Rubinstein [31] in 1997 to estimate the probability of rare events in complex random networks. Since the BOA has a tendency to prematurely converge to local optima [32], this paper embeds the CE method into the BOA to obtain a good balance between exploration and exploitation and improve the BOA's global search capability.

The rest of this paper is structured as follows. In Section 2, the BOA and CE method are briefly introduced, and a study of the improved BOA is presented in Section 3. The results and a discussion of benchmark functions are provided in Section 4. Section 5 solves three classical engineering design problems. In Section 6, conclusions and future research directions are presented.

#### **2. Methods**

#### *2.1. Butterfly Optimization Algorithm*

The BOA is a new nature-inspired meta-heuristic algorithm which was inspired by the food foraging behavior of butterflies [21,22]. In this method, the characteristics of butterflies are ideally summarized as follows:

1. All butterflies are assumed to emit some fragrance that makes the butterflies attractive to each other;

2. Each butterfly moves randomly or toward the best butterfly—i.e., the one which emits the most fragrance;

3. The stimulus intensity of a butterfly is determined by the landscape of the objective function.

In the BOA, the perceived magnitude of the fragrance (*f*) is defined as a function of the stimulus physical intensity as follows:

$$f = \mathfrak{cl}^a,\tag{5}$$

where *c* ∈ [0, ∞] is the sensory modality; *I* is the stimulus intensity, which is associated with the encoded objective function; and *a* ∈ [0, 1] is the power exponent dependent on modality, which represents the varying degree of fragrance absorption.

There are two key steps in the BOA: global search and local search. The former can make the butterflies move toward the best butterfly, which can be represented as

$$\mathbf{x}\_{i}^{t+1} = \mathbf{x}\_{i}^{t} + (r^{2} \times \mathbf{g}^{\*} - \mathbf{x}\_{i}^{t}) \times f\_{i},\tag{6}$$

where *x<sup>t</sup> <sup>i</sup>* is the position of the *i*th butterfly at time *t*. Here, *g*<sup>∗</sup> represents the current best position. *fi* represents the fragrance of the *i*th butterfly, and *r* is a random number in [0, 1].

The latter is implemented through a local random walk, which can be represented as

$$\mathbf{x}\_{i}^{t+1} = \mathbf{x}\_{i}^{t} + (r^{2} \times \mathbf{x}\_{j}^{t} - \mathbf{x}\_{k}^{t}) \times f\_{i\nu} \tag{7}$$

where *x<sup>t</sup> <sup>j</sup>* and *<sup>x</sup><sup>t</sup> <sup>k</sup>* are the positions of the *j*th and *k*th butterflies in the solution space, respectively, and *r* is a random number in [0, 1].

Additionally, a switch probability *p* is used in the BOA to switch between a common global search and an intensive local search. Based on the above, the original BOA can be implemented with pseudo-code as shown in Algorithm 1:



BOA shows better performance compared to some other optimization algorithms [21,22] and has attracted the attention of many researchers due to its simplicity and interesting nature-inspired interpretations as an optimization approach for global optimization problems and for various real-life applications [32,33].

#### *2.2. The Cross-Entropy (CE) Method*

Based on Monte Carlo technology, Rubinstein [31] developed the CE method in 1997 and uses the cross-entropy or Kullback–Leibler divergence to measure the distance between two sampling distributions, solve an optimization problem by minimizing this distance, and obtain the optimal parameters of probability distribution. The CE method has good global search capability, excellent adaptability, and strong robustness. It is frequently used for many complex optimization problems such as continuous multi-extremal optimization [34], multi-objective optimization [35], combination optimization [36,37], vehicle routing problems [38], energy-efficient radio resource management [39], and complex optimization problems from other fields [40–43].

A general optimization problem can be stated as follows:

$$\min\_{\mathbf{x}\in\mathcal{X}} S(\mathbf{x}),\tag{8}$$

where X is a finite set of states, and *S* is a real-valued performance function on X .

Next, we can construct a probability distribution estimation problem associated with the above problem, and the auxiliary problem can be defined as follows:

$$l(\gamma) = P\_{\mathbb{H}}(S(X) \le \gamma) = E\_{\mathbb{H}}[I\_{\{S(X) \le \gamma\}}]\_{\prime} \tag{9}$$

where *Pu* is the probability measure under which the random state *X* has some probability density function *f*(*x*; *u*) on X ; *Eu* is the corresponding expectation operator; *γ* represents a threshold parameter; and *I* represents the indicator function, whose value is 1, if *S*(*X*) ≤ *γ*, and 0, otherwise. The importance sampling approach is used to reduce the sample size in the CE method. Consequently, Equation (9) can be rewritten as follows:

$$I(\gamma) = \frac{1}{N} \sum\_{i=1}^{N} I\_{\mathbb{S}(X)} \le \gamma \frac{f(\mathbf{x}^i; \mathbf{v})}{\mathcal{g}(\mathbf{x}^i)},\tag{10}$$

where *x<sup>i</sup>* represents a random sample from *f*(*x*; *v*) with importance sampling density *g*(*x*). the Kullback–Leibler divergence, i.e., the cross-entropy is introduced to measure the distance between two sampling distributions for obtaining the optimal importance sampling density. Thus, the optimal density *g*∗(*x*) can be found by minimizing the Kullback–Leibler divergence, which is equivalent to solving the following optimization problem [34]:

$$\min\_{\upsilon} \frac{1}{N} \sum\_{i=1}^{N} I\_{\mathcal{S}(X) \le \gamma} \ln f(\mathbf{x}^{i}; \upsilon). \tag{11}$$

In order to improve the convergence speed of the CE method, adaptive smoothing *v*ˆ*<sup>k</sup>* is demoted by *<sup>v</sup>* , as follows:

$$
\vartheta\_{k+1} = \alpha \widetilde{\upsilon} + (1 - \alpha) \vartheta\_{k^\*} \tag{12}
$$

where 0 ≤ *α* ≤ 1 is a smoothing parameter.

The CE method for optimization problems is described in Algorithm 2.


#### **3. Hybrid BOA-CE Method**

This section presents the details of the improved algorithm, named BOA-CE. A meta-heuristic method should have two main functions—exploration and exploitation—and should try to find a proper balance between them for better performance [44]. The nature-inspired BOA has shown the advantages of excellent local search capability and fast convergence; however, it tends to fall into a local optimum rather than finding the global optimum [21,22,32,33]. In order to improve the global search capability of the BOA, and, based on a co-evolutionary technique, this paper proposes the BOA-CE algorithm by embedding the CE method into it. The new method contains BOA and CE optimization operators, and co-updates the BOA population *P* and the CE sample *S* using co-evolutionary technology in each iteration. The CE operator obtains the initial probability parameters using the BOA population *P* in order to improve its convergence rate while the BOA operator updates its population *P* and the best butterfly using the elite sample *Se* of the CE operator to enhance the population diversity. The pseudo-code of the BOA-CE hybrid algorithm is described in Algorithm 3. Figure 1 presents the flow chart of the BOA-CE algorithm, and the co-evolutionary process between BOA operator and CE operator can be clearly seen from Figure 1.

In addition, as an example, we use the BOA-CE algorithm to find the global optimum of the 2D Sphere function in order to more clearly describe the co-evolutionary process of the BOA operator and the CE operator:

$$f(\mathbf{x}) = \mathbf{x}\_1^2 + \mathbf{x}\_2^2, \quad -100 \le \mathbf{x}\_1, \mathbf{x}\_2 \le 100,\tag{13}$$

which has a global minimum *f* <sup>∗</sup> = 0 at *x*<sup>∗</sup> = (0, 0)*T*. Now, let us use the BOA-CE algorithm to find the optima and this process is summarized as follows: (1) the values of the parameters of the BOA operator are *c* = 0.01, *a* = 0.1, *p* = 0.8, the population size *n* = 40, and the maximum number of iteration *t*<sup>1</sup> = 50, while the values of the parameters of the CE operator are *α* = 0.8, the sample size *N* = 80, the elite sample size *Ne* = 30, and the maximum number of iteration *t*<sup>2</sup> = 2; (2) initializing the population of the BOA operator *P*, calculating *fi*(*x*) = *x*<sup>2</sup> *<sup>i</sup>*<sup>1</sup> + *<sup>x</sup>*<sup>2</sup> *<sup>i</sup>*2, *i* = 1, 2, ··· , 40, and finding the current best *<sup>f</sup>* <sup>∗</sup> = 13.5036 and *<sup>x</sup>*<sup>∗</sup> = (−3.5586, −0.9164)*T*; (3) generating a rand number *r*, and updating *xi* with *xnew <sup>i</sup>* = *<sup>x</sup>old <sup>i</sup>* + (*r*<sup>2</sup> × *<sup>f</sup>* <sup>∗</sup> − *<sup>x</sup>old <sup>i</sup>* ) × *fi* if *<sup>r</sup>* ≤ *<sup>p</sup>* and *<sup>x</sup>new <sup>i</sup>* = *<sup>x</sup>old <sup>i</sup>* + (*r*<sup>2</sup> × *<sup>x</sup>old j* − *xold <sup>k</sup>* ) × *fi* otherwise for each butterfly in the *P*; (4) calculating the mean *μ* = (−4.9571, −11.4528) and standard *σ* = (60.2592, 53.9977) deviation of P, and using them to initialize the probability parameter *v* = (−4.9571, −11.4528, 60.2592, 53.9977); (5) generating a sample *S* by the probability parameter *v*, and calculating *fi*(*x*) = *x*<sup>2</sup> *<sup>i</sup>*<sup>1</sup> + *<sup>x</sup>*<sup>2</sup> *<sup>i</sup>*2, *i* = 1, 2, ··· , 80, and finding the new current best *f* <sup>∗</sup> = 4.4290 and *x*<sup>∗</sup> = (1.8693, 0.9668)*T*; (6) updating the probability parameter *v* using the elite sample, and regenerating a sample *S* using the new parameter *v*, and finding the new current best *f* <sup>∗</sup> = 2.4542 and *x*<sup>∗</sup> = (0.8951, 1.2857)*T*; (7) co-updating the population P using the sample S and the population P; (8) repeating the above operations until the termination condition (*t*<sup>1</sup> > 50) is

met, and outputting the approximate optimal solution *<sup>x</sup>*<sup>∗</sup> = (1.1154 × <sup>10</sup>−31, −2.5727 × <sup>10</sup>−32)*<sup>T</sup>* and function value *<sup>f</sup>* <sup>∗</sup> = 1.3104 × <sup>10</sup><sup>−</sup>62. Figure <sup>2</sup> clearly describes the of co-updating the current best*<sup>f</sup>* <sup>∗</sup> by the BOA operator and CE operator.


The BOA-CE hybrid algorithm has two main operators: BOA and CE. The BOA operator has two inner loops for the butterflies population size *n*1, and one outer loop for iteration *t*1. While the CE operator has one inner loops for the sample size *n*2, and two outer loops for iterations *t*<sup>1</sup> and *t*2, respectively. Therefore, for our proposed hybrid algorithm (BOA-CE), the time complexity is *O*(2*n*1*t*<sup>1</sup> + *n*2*t*1*t*2). It is linear in terms of *t*1*t*2, which represents the total number of iterations in the BOA-CE.

**Figure 1.** The flow chart of the BOA-CE algorithm.

**Figure 2.** The BOA operator and CE operator co-update the current best *f* ∗ of the 2D Sphere function.

#### **4. Experiment and Results**

In this section, a total of 19 test functions [45,46] with different characteristics were employed to benchmark the performance of the proposed BOA-CE algorithm from different perspectives [9,15,16,20,28]. The test functions can be divided into three classes: unimodal (Table 1), multimodal (Table 2), and composite functions (Table 3). Generally, unimodal functions are employed to benchmark the exploitation and convergence of a method, while multimodal functions are selected to evaluate the performance of exploration and local optima avoidance [9,20]. In contrast, composite functions can be used to evaluate the combined performance of exploration and exploitation.


**Table 1.** Unimodal benchmark functions.



To verify the results, the proposed BOA-CE algorithm was compared against a number of well-known and recent algorithms: Genetic Algorithm (GA) [10], PSO [12], BA [14], GWO [15], CSA [19], SSA [20], and BOA [31]. The variants were coded in the MATLAB R2018b, running on a PC with an Intel Core i7-8700 processor (Gainesville, FL, USA), a 3.19 GHz CPU, and 16 GB of RAM. To provide a fair comparison, the following test experimental conditions and settings were used: (1) the population size of the BOA operator in the BOA-CE algorithm was set to 40, while the CE operator's sample size was 80. The maximum number of iterations of the BOA operator and CE operator were 500 and 2, respectively; (2) the other method's population sizes were 100 for fair comparison, and their maximum number of iterations was 1000; (3) all of the other parameters in each algorithm for comparison were as the same as those used in the original references [10,12,14,15,19,20,31]. Each of the algorithms was run 30 times on each of the test functions, and the results are shown in Tables 4–6. The average value and standard deviation of the best solution obtained by each method is given to compare their overall performance. The winner (best value) is highlighted in bold.


**Table 3.** Composite benchmark functions.

<sup>1</sup> Mathematical formulation of Sphere, Ackley, Rastrigin, Weierstrass, and Griewank can be found in Appendix A A1.

#### *4.1. Results of the Algorithms Using Unimodal Test Functions*

The results of the application of the algorithms to unimodal test functions are shown in Table 4. The best results are highlighted in bold. From Table 4, it can be found that: (1) the proposed method outperforms GA, PSO, BA, CSA, SSA, and BOA for the majority of the test functions; and (2) BOA-CE and GWO each provide the best results for three of these problems. This is due to the co-evolutionary technology that are adopted between the BOA and CE operators to enhance exploitation.

The progress of the average best value over 30 runs for some of the benchmark functions, such as F1, F2, F6, and F7, is illustrated in Figure 3. The figure clearly shows that the proposed BOA-CE algorithm rapidly converges towards the optimum and exploits it accurately. This benefit is mainly due to the algorithm's excellent exploitation.

#### *4.2. Results of the Algorithms Using Multimodal Test Functions*

The results of the application of the algorithms to multimodal test functions are shown in Table 5. Table 5 shows that the BOA-CE method outperforms the other approaches in the majority of the test cases. Multimodal test functions are employed to benchmark the performance in terms of exploration and local optima avoidance. Therefore, these results prove that the BOA-CE algorithm has an excellent exploration which helps it to explore the promising regions of the search space. Additionally, the local optima avoidance of the improved method is also excellent, since it is able to avoid all of the local optima and find the global optimum for the majority of the multimodal test functions.





**Table 6.** Results of composite benchmark functions.

**Figure 3.** Convergence of algorithms when applied to some unimodal benchmark functions.

Figure 4 shows the convergence curves of the algorithms when applied to some of the multimodal test functions, such as F10, F11, F12, and F13. The figure shows that the BOA-CE algorithm has the fastest convergence for the majority of the multimodal test functions.

**Figure 4.** Convergence of algorithms on some unimodal benchmark functions.

#### *4.3. Results of the Algorithms on Composite Test Functions*

The results of composite test functions are reported in Table 6. From Table 6, we can find that the BOA-CE algorithm outperforms others in most of the composite test functions. Considering the characteristics of composite test functions and these results, it may be stated that the BOA-CE algorithm appropriately balances exploration and exploitation to focus on the high-performance areas of the search space.

The convergence curves of the methods for some of the composite test functions are shown in Figure 5. Figure 5 clearly shows that the BOA-CE algorithm has the fastest convergence on the composite test functions.

In addition, Table 7 presents the time consumption of all the algorithms to solve the 19 benchmark functions. From Table 7, it can be seen that the least total time cost is PSO to solve all of the benchmark functions, followed by CSA, and the BOA-CE is ranked third, which is better than BOA.

The results of three test functions show that the local optima avoidance of the improved method was enhanced by embedding the cross-entropy method. This global stochastic optimization method enables the BOA-CE algorithm to achieve a proper balance between exploration and exploitation, and enhance its global search capability. Furthermore, the new method is particularly outstanding in solving multimodal function optimization problems.This provides a new method for solving complex engineering design optimization problems.

#### *4.4. Analysis of the Hybrid BOA-CE Algorithm*

The main reasons for the excellent performance of the proposed BOA-CE algorithm in solving global optimization problems may be stated as follows:


diversity to the BOA operator so that it can effectively avoid falling into a local optimum and improve its global search capability.

• The BOA-CE algorithm employs a co-evolutionary technique to co-update the BOA operator's population and the CE operator's probability parameters. This enables the improved method to obtain an appropriate balance between exploration and exploitation and have more superior performance in terms of exploitation, exploration, and local optima avoidance when solving complex optimization problems.

**Figure 5.** Convergence of algorithms when applied to some composite benchmark functions.


**Table 7.** Time consumption of BOA-CE solving 19 benchmark functions.

#### **5. Using the BOA-CE Algorithm for Classical Engineering Design Problems**

In this section, the BOA-CE algorithm was tested on three classical constrained engineering design problems: a tension/compression spring, a welded beam, and a pressure vessel [9,15,19,20,28,31,47–52]. Constraint handling is a challenging task for a method when solving these problems. Penalty functions are divided into different types: static, dynamic, annealing, adaptive, co-evolutionary, and death penalty [47]. For the sake of simplicity, we equipped the BOA-CE algorithm with a death penalty function to handle constraints. Table 8 shows the BOA-CE parameters chosen to solve these design problems, where *N* is the population sizes of the BOA and CE operator, and *Itermax* represents the maximum number of iterations.

**Table 8.** Parameters of the BOA-CE algorithm for solving design problems.


#### *5.1. Tension/Compression Spring Design*

The objective of this problem is to minimize the weight of the tension/compression spring shown in Figure 6. It contains three design variables: the wire diameter (*d*), mean coil diameter (*D*), and number of active coils (*N*), which are subject to one linear and three nonlinear inequality constraints on shear stress, surge frequency, and deflection. The optimization problem is formulated as follows:

$$\min f(\mathbf{x}) = (\mathbf{x}\mathbf{z} + \mathbf{2})\mathbf{x}\mathbf{z}\_1^2 \tag{14}$$

subject to

$$\log\_1(\mathbf{x}) = 1 - \frac{\mathbf{x}\_2^3 \mathbf{x}\_3}{71785 \mathbf{x}\_1^4} \le 0,\tag{15}$$

$$\log g\_2(\mathbf{x}) = \frac{4\mathbf{x}\_2^2 - \mathbf{x}\_1\mathbf{x}\_2}{12566(\mathbf{x}\_2\mathbf{x}\_1^3 - \mathbf{x}\_1^4)} + \frac{1}{5108\mathbf{x}\_1^2} - 1 \le 0,\tag{16}$$

$$\text{g3}(\mathbf{x}) = 1 - \frac{140.45\mathbf{x}\_1}{\mathbf{x}\_2^2 \mathbf{x}\_3} \le 0,\tag{17}$$

$$g\_4(x) = \frac{x\_1 + x\_2}{1.5} - 1 \le 0,\tag{18}$$

where 0.05 ≤ *x*<sup>1</sup> ≤ 2.00, 0.25 ≤ *x*<sup>2</sup> ≤ 1.30, and 2.00 ≤ *x*<sup>3</sup> ≤ 15.00.

**Figure 6.** The tension/compression spring design problem.

This problem has been popular among researchers and has been optimized using various meta-heuristic algorithms. To solve this problem, Coello employed GA [48], He and Wang employed PSO [49], Gandomi et al. used BA [50], Mirjalili employed GWO [15] and SSA [20], Lee and Geem [51] used HS, and Askarzadeh [19] used CSA. The best weight values are highlighted in bold. Table 9 shows a comparison of the results obtained using the BOA-CE algorithm and those obtained using the other algorithms. It can be seen that the BOA-CE and CSA algorithms find a design with the minimum weight for this problem and outperform all other algorithms.


**Table 9.** A comparison of results for the tension/compression spring design problem.

#### *5.2. Welded Beam Design*

The objective of this problem is to minimize the fabrication cost of the welded beam shown in Figure 7. It contains four design variables: the thickness of the weld (*h*), length of the clamped bar (*l*), height of the bar (*t*), and thickness of the bar (*b*), which are subject to two linear and five nonlinear inequality constraints on shear stress, bending stress in the beam, buckling load, and end deflection of the beam. The optimization problem can be stated as follows:

$$\min f(\mathbf{x}) = 1.1047 \mathbf{x}\_1^2 \mathbf{x}\_2 + 0.04811 \mathbf{x}\_3 \mathbf{x}\_4 (14.0 + \mathbf{x}\_2),\tag{19}$$

subject to

$$g\_1(\mathbf{x}) = \tau(\mathbf{x}) - \tau\_{\max} \le 0,\tag{20}$$

$$\log\_2(\mathbf{x}) = \sigma(\mathbf{x}) - \sigma\_{\text{max}} \le 0,\tag{21}$$

$$\mathcal{g}\_{\mathcal{B}}(\mathbf{x}) = \delta(\mathbf{x}) - \delta\_{\text{max}} \le 0,\tag{22}$$

$$\mathbb{S}\_4(\mathbf{x}) = \mathbf{x}\_1 - \mathbf{x}\_4 \le \mathbf{0},\tag{23}$$

$$\mathcal{g}\_5(\mathbf{x}) = P - P\_\varepsilon(\mathbf{x}) \le 0,\tag{24}$$

$$\mathbf{g}\_6(\mathbf{x}) = 0.125 - \mathbf{x}\_1 \le \mathbf{0},\tag{25}$$

$$\log\_7(\mathbf{x}) = 1.10471\mathbf{x}\_1^2 + 0.04811\mathbf{x}\_3\mathbf{x}\_4(14.0 + \mathbf{x}\_2) - 5.0 \le 0,\tag{26}$$

 $\text{where } \tau(\mathbf{x}) = \sqrt{(\tau')^2 + 2\tau'\tau''\frac{x\_2}{2\hbar} + (\tau'')^2}, \ \tau' = \frac{p}{\sqrt{2x\_1x\_2}}, \ \tau'' = \frac{\hbar\overline{R}}{\frac{1}{\tau}}, M = P(L + \frac{x\_2}{2}), \ \overline{R} = \frac{p}{\sqrt{2x\_1x\_2}}, \ \tau'' = \frac{\hbar\overline{R}}{\frac{1}{\tau}}, \ \tau'' = \frac{\overline{R}}{\frac{1}{\tau}}$  $\text{2x1}\_2(\frac{x\_2}{4} + (\frac{x\_1 + x\_3}{2})^2), \ \sigma(\mathbf{x}) = \frac{\overline{6PL}}{\frac{1}{\tau\mathbf{x}\_1\overline{x}\_3}}, \ \delta(\mathbf{x}) = \frac{\overline{6PL}^3}{\frac{1}{\tau\mathbf{x}\_1\overline{x}\_3}}, P\_\mathbf{c}(\mathbf{x}) = \frac{4.015\overline{\mathcal{E}}\frac{x\_2\overline{x}\_4}{36}}{L^2}$  $\text{3x1}\_2(\frac{x\_2}{4} + (\frac{x\_1}{2} - \tau)\overline{R}), \ \tau = 0.001\overline{\mathcal{E}}\text{b}, L = 14\text{ in.}$  $\sigma\_{\text{max}} = 0.25 \text{ in.}$  $\mathcal{E} = 30 \times 10^6 \text{ psi}, \ G = 12 \times 10^6 \text{ psi}, \ \tau\_{\text{max}} = 30,000 \text{ psi}, \ \tau\_{\text{max}} = 0.1 \le \mathbf{x}\_1, \mathbf{x}\_4 \le 2, \text{ and } 0.1 \le \mathbf{x}\_2, \mathbf{x}\_3 \le 10.$ 

The welded beam design problem has been tackled using many meta-heuristic algorithms, such as GA [48], PSO [49], BA [50], GWO [15], HS [51], WOA [18], and CSA [19]. The optimization results using the BOA-CE algorithm are compared with those from the literature in Table 10. The minimum cost is highlighted in bold. The table shows that BOA-CE outperforms all other algorithms except CSA, with only a slight difference in its result.

**Figure 7.** The welded beam design problem.



#### *5.3. Pressure Vessel Design*

The objective of this problem is to minimize the total cost of a pressure vessel considering the cost of material, forming and welding shown in Figure 8. There are two discrete and two continuous design variables: thickness of the shell (*Ts*), thickness of the head (*Th*), inner radius (*R*) and length of the cylindrical section of the vessel (*L*), not including the head, which are subjected to three linear and one nonlinear inequality constraints. The thicknesses of the variables are integer multiples of 0.0625 inches. This optimization problem can be mathematically formulated as follows:

$$\min f(\mathbf{x}) = 0.6224 \mathbf{x}\_1 \mathbf{x}\_3 \mathbf{x}\_4 + 1.7781 \mathbf{x}\_2 \mathbf{x}\_3^2 + 3.1661 \mathbf{x}\_1^2 \mathbf{x}\_4 + 19.84 \mathbf{x}\_1^2 \mathbf{x}\_3. \tag{27}$$

subject to

$$\mathbf{x}\_1 \mathbf{g}\_1(\mathbf{x}) = -\mathbf{x}\_1 + 0.0193\mathbf{x}\_3 \le \mathbf{0},\tag{28}$$

$$\mathcal{g}\_2(\mathbf{x}) = -\mathbf{x}\_3 + 0.00954\mathbf{x}\_4 \le 0,\tag{29}$$

$$\log\_3(\mathbf{x}) = -\pi \mathbf{x}\_3^2 \mathbf{x}\_4 - \frac{4}{3} \pi \mathbf{x}\_3^3 + 1296000 \le 0,\tag{30}$$

$$\mathbf{g\_4(x)} = \mathbf{x\_4} - \mathbf{240} \le \mathbf{0},\tag{31}$$

where 0 < *x*1, *x*<sup>2</sup> ≤ 99, 0 < *x*3, *x*<sup>4</sup> ≤ 200.

In the past, this problem is solved by GA [48], PSO [49], BA [50], GWO [15], HS [51], WOA [18], and DE [52]. Table 11 shows optimization results of BOA-CE are compared with those found by other methods. It can be seen that BOA-CE and BA outperform all other algorithms except GWO. However the value of *Th* obtained by GWO is not an integer multiple of 0.0625 inches and thus does not satisfy the constraint.

**Figure 8.** The tension/compression spring design problem.

**Table 11.** A comparison of results for the tension/compression spring design problem.


#### **6. Conclusions**

In order to improve the global search ability of the BOA, an improved BOA algorithm, named BOA-CE, was constructed by embedding the CE method into the BOA using a co-evolutionary technique. A total of 19 test functions were used to evaluate the performance of the new method in terms of its exploration, exploitation, local optima avoidance, and convergence rate. The results of the test functions demonstrated that the proposed algorithm can effectively avoid falling into a local optima and has excellent local and global search capacity. This is mainly due to the co-evolutionary technique, which enables the new method to obtain an appropriate balance between exploration and exploitation and has more superior performance when solving complex function optimization problems. Furthermore, the paper also considered the solving of three classical engineering problems using the hybrid algorithm. The comparative results show that the BOA-CE algorithm provides very competitive results when solving real problems with unknown search spaces. In future research, a discrete version of the BOA-CE algorithm will be constructed to solve combinatorial optimization problems.

**Author Contributions:** G.L. designed the algorithm, implemented all experiments, and wrote the manuscript. F.S. revised the manuscript. P.Z. conducted the literature review and analyzed the data. C.L. revised the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant No. 71761012), the Natural Science Foundation of Anhui Province (Grant No. 1808085MG224), and the Natural Science Foundation of West Anhui University (Grant No. WXZR201818).

**Acknowledgments:** The authors are grateful for the support provided by the Risk Management and Financial Engineering Lab at the University of Florida, Gainesville, FL 32611, USA.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

See Table A1.


*D*: dimensions.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
