Memetic-Based Biogeography Optimization Model for the Optimal Design of Mechanical Systems

1. Introduction

Engineering problems have become increasingly complex, making their analytical resolution practically impossible to complete in many situations. Numerical solution methods have therefore been given greater focus. Optimization, fitting into this category, has become one of the most pivotal subjects in the day-to-day of a modern engineer. Specifically, because there is not an area in engineering where optimization isn't involved in some way. The optimal design of engineering systems has therefore undergone many important developments in recent years. One of the biggest developments was the introduction of the genetic algorithm (GA) by John Holand [1], the first evolutionary algorithm with structured operators, an event that revolutionized the entire optimization environment. Since then, the scientific community has been perfecting this method, creating a large number of proposals for new operators, based on all kinds of schemes and heuristics, all to solve optimization problems. The genetic algorithm uses some ideas from other scientific fields, such as genetics and biology, and applies them mathematically, in the form of an algorithm, to different areas of engineering. There are also other methods inspired by other areas, such as memetic algorithms (MA), which are based on biology (genetic evolution) [2] and learning processes (cultural evolution) [3] - [4], and the algorithm that will be presented in this article, Biogeography-Based Optimization (BBO).

Biogeography is the study of the distribution of species and ecosystems in geographical space over geological time. Biogeographical research combines information and ideas from different fields, from the physiological and ecological constraints on the spreading of organisms to the geological and climatological phenomena that occur on a global scale in evolutionary periods.

Biogeography-based optimization is an evolutionary algorithm (EA), introduced in 2008 by Dan Simon [5], which optimizes a function by iteratively and stochastically improving its solutions according to a quality measure, or fitness function. Typically, this method is used to study multivariate problems and, since it does not require the calculation of a gradient function, it can also be used to optimize discontinuous functions. BBO optimizes a problem by maintaining a population of candidate solutions and creating new solutions through the migration process.

The main objective of this work is to develop and implement the BBO algorithm, following the readily available literature, and, later on, modify and hybridize this method with other concepts of optimization, such as crossover and local search operators or other propositions used in other evolutionary algorithms. The BBO method will be applied to several problems and examples from the area of engineering and optimal design so that it is possible to study its behavior and the effects of the new modifications, introduced in the process of hybridization.

2. BBO Approach Model

The science of biogeography, created in the 19th century, was considered to be descriptive and historical, but in 1967 Robert MacArthur and Edward Wilson published "The Theory of Island Biogeography" [6], where they tried to describe this phenomenon in mathematical equations. These models have tried to describe how species move from one island to another, how new species appear, and how they become extinct. In this case, the term island should not be taken at face value, but as a place that is geographically isolated from other habitats. Therefore, the term habitat should be used to better describe these locations.

Sites with good life retention characteristics are associated with habitats with a high HSI (habitat suitability index). Factors such as the presence of drinking water, which are related to the HSI, are referred to as suitability index variables (SIVs). In mathematical terms, the SIV can be called the independent variables and the HSI the dependent variables. Then, the habitats with high HSI will have a lot of species, while the habitats with low HSI will have less species. In real life, the emigration of a species may not lead to its extinction, but, for this case, the emigration of a species will lead to the extinction of the species in that habitat. Low HSI habitats will have a high immigration rate, due to their lesser population numbers and a non-competitive environment. On the contrary, high HSI habitats will have a high emigration and low immigration rate.

Good solutions are associated with habitats with high HSI, while habitats with low HSI are associated with bad solutions. Solutions with high HSI are more resistant to change, so solutions with high HSI tend to share characteristics with low HSI solutions. The diversity of individuals is related to HSI, so when an individual arrives in a habitat with low HSI, the habitat's HSI tends to increase. This approach has been given the name Biogeography-Based Optimization (BBO), which was created and presented by Dan Simon [5]. Figure 1 presents the flow diagram of the method.

2.1. Mathematical Model

Figure 2 shows a model of the abundance of species in a single habitat. The emigration rate $\mu$ and the immigration rate $\lambda$ are functions of the number of species in the habitat. Looking at the immigration curve, it is possible to note that its maximum will occur when the migration rate of the habitat is $I$, that happens when the number of species of the habitat reaches zero. As the number of species increases, the immigration process will become more competitive, which will lead to a decrease in the number of species able to complete the process, resulting in a reduction in the immigration rate. The maximum number of species present in a habitat is defined as $S_{\max }$, at which point the immigration rate will be zero. The emigration curve, on the other hand, will have its maximum $E$ when the population of the habitat reaches its maximum and will be zero when there are no more species in the habitat.

The equilibrium point, $k_0$, is obtained when the rates of emigration and immigration are equal. This point can occasionally differ due to temporary effects such as diseases. Figure 2 represents both rates as straight lines, but in reality, they can be complex curves. However, this representation is a good approximation of the migration process.

Now let's consider a probability $P_s$ that the habitat contains exactly $S$ species. The probability changes with time, from $t$ to $t+\mathsf{\Delta }t$, following the equation below,

$$P_s\left(t+\mathsf{\Delta }t\right)=P_s\left(t\right)\left(1-\lambda \mathsf{\Delta }t-{\mu }_s\mathsf{\Delta }t\right)+P_{s-1}{\lambda }_{s-1}\mathsf{\Delta }t+P_{s+1}{\mu }_{t+1}\mathsf{\Delta }t$$

(2.1)

where ${\lambda }_s$ and ${\mu }_s$ are the immigration and emigration rates when there are $S$ species in the habitat.

Assuming that the time interval $\mathsf{\Delta }t$ is small enough, the probability of more than one emigration or immigration can be ignored. Using a limit in the above equation, with the time interval tending to zero, the following equation emerges,

$${\dot{P}}_s=\left\{\begin{array}{c}-\left({\lambda }_s+{\mu }_s\right)P_s+{\mu }_{s+1}{P}_{s+1},S=0\\ -\left({\lambda }_s+{\mu }_s\right)P_s+{\lambda }_{s-1}{P}_{s-1}+{\mu }_{s+1}{P}_{s+1},1<S\le S_{\max }-1\\ -\left({\lambda }_s+{\mu }_s\right)P_s+{\lambda }_{s-1}{P}_{s-1},S=S_{\max }\end{array}\right.$$

(2.2)

Aiming to simplify the notation Simon [5], decided to perform the following changes, $n=S_{\max }$ and $P={\left[P_0\dots P_n\right]}^T$. Then it is possible to write the previous equation as,

$$\dot{P}=AP$$

(2.4)

where the $A$ matrix is given by

$$A=\left[\begin{array}{ccccc}-\left({\lambda }_0+{\mu }_0\right)& {\mu }_1& 0& \dots & 0\\ {\lambda }_0& -\left({\lambda }_1+{\mu }_1\right)& {\mu }_2& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ ⋮& \ddots & {\lambda }_{n-2}& -\left({\lambda }_{n-1}+{\mu }_{n-1}\right)& {\mu }_n\\ 0& \dots & 0& {\lambda }_{n-1}& -\left({\lambda }_n+{\mu }_n\right)\end{array}\right]$$

(2.3)

Considering the lines from Figure 2 the migration rates can be described as,

$${\mu }_s={\displaystyle \frac{Es}{n}}$$

(2.5)

$${\lambda }_s=I\left(1-{\displaystyle \frac{s}{n}}\right)$$

(2.6)

The equations above can be presented in a different form if the number of species is replaced by the value calculated in the fitness function [8],

$${\lambda }_i=I{\displaystyle \frac{f_{\max }-f\left(H_i\right)+\epsilon }{f_{\max }-f_{\min }+\epsilon }}$$

(2.7)

$${\mu }_i=E{\displaystyle \frac{f\left(H_i\right)-f_{\min }+\epsilon }{f_{\max }-f_{\min }+\epsilon }}$$

(2.8)

where $H_i$ represents the solution in study.

Following the same ideas, the calculation of the migration rates can also be given by the ranking of the solutions, after they have been placed in descending order as [8].

$${\lambda }_i=I{\displaystyle \frac{i}{n}}$$

(2.9)

$${\mu }_i=E\left(1-{\displaystyle \frac{i}{n}}\right)$$

(2.10)

Now considering the case where $E=I$, we have,

$${\lambda }_s+{\mu }_s=E=I$$

(2.11)

and the matrix A transforms into,

$$A=E\left[\begin{array}{ccccc}-1& {\displaystyle \frac{1}{n}}& 0& \dots & 0\\ {\displaystyle \frac{n}{n}}& -1& {\displaystyle \frac{2}{n}}& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ ⋮& \ddots & {\displaystyle \frac{2}{n}}& -1& {\displaystyle \frac{n}{n}}\\ 0& \dots & 0& {\displaystyle \frac{1}{n}}& -1\end{array}\right]=EA\prime$$

(2.12)

It is possible to note that zero is an eigenvalue of matrix $A\prime$ with the corresponding eigenvectors,

$$v=\left[v_1\dots v_{n+1}\right]$$

(2.13)

$$v_i=\left\{\begin{array}{c}{\displaystyle \frac{n!}{\left(n-1-i\right)!\left(i-1\right)!}},(i=1,\dots ,i^{\prime })\\ v_{n+2-i},\left(i=i^{\prime }+1,\dots ,n+1\right)\end{array}\right.$$

(2.14)

The value, in steady state, of the probability of the number of each species is given by,

$$P\left(\infty \right)={\displaystyle \frac{v}{{\sum }_{i=1}^{n+1}{v}_i}}$$

(2.15)

computation of the probability, which depended on the migration rates, previously calculated, as shown in equation 2.16. This formulation solves a problem that restricted the number of possible initial solutions.

$$P_i=\left\{\begin{array}{c}{P}_0={\displaystyle \frac{1}{1+{\sum }_{k=1}^n\frac{{\lambda }_0{\lambda }_1\dots {\lambda }_{k-1}}{{\mu }_1{\mu }_2\dots {\mu }_k}}},i=0\\ P_i={\displaystyle \frac{{\lambda }_0{\lambda }_1\dots {\lambda }_{i-1}}{{\mu }_1{\mu }_2\dots {\mu }_i\left(1+{\sum }_{k=1}^n\frac{{\lambda }_0{\lambda }_1\dots {\lambda }_{k-1}}{{\mu }_1{\mu }_2\dots {\mu }_k}\right)}},1\le i\le n\end{array}\right.$$

(2.16)

Simon [5] recommends the use of $E=I=1$, $N=50$ and $m_{\max }=0.01$ as the initial parameters for the method.

2.2. Migration Operator

Supposing that the population of solutions can be represented as a vector of integers. Each of the numbers is considered a SIV. The HSI, which is analogous to GAs fitness, can be used in other processes of optimization based on populations. High HSI solutions represent the habitats with several species and low HSI solutions represent the habitats with lesser amounts of solutions. It is possible to assume that each solution has an identical migration curve (with $E=I$), but the value of $S$ represented by the solution remains dependent only on its HSI.

The migration rates of each solution are used to, probabilistically, share information between habitats. The immigration rate is therefore used to determine whether the SIV of each solution will change. If, in a given solution $S_i$, an SIV is randomly selected to be modified, then the emigration rate of the other solutions will decide from which of the solutions an SIV, located in the same position, will migrate to solution $S_i$. This process is quite similar to other reproduction processes from GAs, but this method differs from them in a very important aspect. In the GA, reproduction by recombination is used to create new solutions, while BBO migration is used to modify existing solutions. It is therefore possible to say that, unlike the reproduction process of the GA, the migration of the BBO is an adaptive process. In algorithmic terms, the migration operator can be described in the following way:

Algorithm 1: Migration’s Pseudocode
1	For $d=1$ to $D$ do
2		If $\mathrm{rand}<{\lambda }_i$ then
3			Select another habitat ($H_j$) with a probability $\propto {\mu }_j$
4			Define $H_i\left(d\right)=H_j\left(d\right)$
5		End If
6	End For

where $d$ represents the variable under study and $D$ the size of the problem.

Normally, the use of linear models (Figure 2 and equations 2.5-2.10) to describe the migration curves is sufficient, but there are also other different options to describe this operator. Ma [7] presented several other options for the migration operator, such as, the trapezoidal model. As described in the following equation and represented graphically in Figure 3,

$${\lambda }_s=\left\{\begin{array}{c}I,s\le i\prime \\ 2I\left(1-{\displaystyle \frac{s}{n}}\right),i^{\prime }<s\le n\end{array}\right.$$

(2.17)

$${\mu }_s=\left\{\begin{array}{c}2·{\displaystyle \frac{E}{n}},s\le i\prime \\ E,i^{\prime }<s\le n\end{array}\right.$$

(2.18)

where $i\prime$ is the smallest integer close to $\left(n+1\right)/2$.

Ma [7] also presented some nonlinear models, such as the quadratic (equations 2.19-2.20 and Figure 4) and sinusoidal (equations 2.21-2.22 and Figure 5) models. These models aim to provide a more accurate simulation of real-life migration patterns, by incorporating nonlinearities.

$${\lambda }_s=I{\left(1-{\displaystyle \frac{s}{n}}\right)}^2$$

(2.19)

$${\mu }_s=E{\left({\displaystyle \frac{s}{n}}\right)}^2$$

(2.20)

$${\lambda }_s={\displaystyle \frac{I}{2}}\left(\cos \left({\displaystyle \frac{s\pi }{n}}\right)+1\right)$$

(2.21)

$${\mu }_s={\displaystyle \frac{E}{2}}\left(-\cos \left({\displaystyle \frac{s\pi }{n}}\right)+1\right)$$

(2.22)

Migration, being one of the main points of interest in the method, has been severely modified and improved. Mainly, Ma and Simon [9] presented a new hybrid algorithm, called Blended BBO (B-BBO), that uses some ideas from GAs’ crossover operators, creating a migration operator, that mixes each component of an immigrating solution ($H_i$) with the corresponding component of the emigrating solution ($H_j$), as follows,

$$H_i^{\prime }\left(d\right)=\alpha H_i\left(d\right)+\left(1-\alpha \right)H_j\left(d\right)$$

(2.23)

where $\alpha$ is a number between $\left[\mathrm{0,1}\right]$, which can be a predetermined or random value.

With the standard migration operator, it is possible to improve poor solutions by using the features of the best solutions. But the blending operator makes it possible not only to do this, but also to avoid the degradation of good solutions, since they retain some of their original properties. This operator also allows the creation of new components, promoting population diversity.

2.3. Mutation Operator

Cataclysmic events can drastically change the HSI of a natural habitat, which can make it so that the number of species varies from the equilibrium point $k_0$, with a mutation of SIV. Each element of the population has an associated probability, which indicates the probability of expecting a priori that there is a solution to the problem in question. Solutions with a very low or very high HSI are equally unlikely to mutate and solutions with a medium HSI are relatively likely. If a solution $S$ has a low probability $P_i$, then the probability of it being mutated is high. On the other hand, a solution with a high probability will have fewer opportunities to mutate. This can be implemented, using the following formula,

$$m_i=m_{\max }\left({\displaystyle \frac{1-P_i}{P_{\max }}}\right)$$

(2.24)

where $m_{\max }$ is a parameter defined by the user, that represents the maximum mutation rate and $P_{\max }$ is the maximum of $P_i$ in the population

This method makes it possible to improve solutions with low HSI. Furthermore, with the use of elitism, which guarantees that the best solutions are kept in the study space, mutation adds some guarantee that the algorithm won't get stuck in a local optimum. Algorithmically, this process can be described as shown in the pseudocode below (Algorithm 2). With $u_d$ and $l_d$ representing the upper and lower bounds of the variable, respectively.

Algorithm 2: Pseudocode for BBO’s standard mutation
1	For $d=1$ to $D$ do
2		If $\mathrm{rand}<m_i$ then
3			Define $H_i\left(d\right)=l_d+\mathrm{rand}·\left(u_d-l_d\right)$
4		End If
5	End For

Recently, the mutation process has been simplified so that only habitats with low HSI can be mutated, eliminating the need to calculate the species count probability and mutation rate for each species. Like migration, the mutation operator has undergone several modifications. Hao's taxonomic mutation [10] is one of these modifications. With this operator, solutions with a high HSI are excluded from the process and the type of mutation differs if the HSI of the solutions is medium or low. For medium HSI solutions, since the solution is somewhat valuable, a differential mutation process will be used, using a $P_C$ crossover rate to judge how the solution will be mutated. For the solutions with low HSI a random mutation, similar to BBO’s Standard mutation, but instead of using the problem bounds ($u_d$ and $l_d$) it will use the population’s best and worst solutions, will be used. For solutions with low HSI, a random mutation will be used, similar to the standard BBO mutation, but instead of using the problem limits ($u_d$ and $l_d$), it will use the best and worst solutions in the population. The following algorithm (Algorithm 3) represents the pseudocode for the differential mutation.

Algorithm 3: Pseudocode for differential mutation
1	For $i=1$ to $N_{\mathrm{pop}}$ do
2		Select three solutions ($H_1,H_2,H_3$) and a dimension ($d_r$) randomly
3		For $d=1$ to $D$ do
4			If $\mathrm{rand}<m_i$ then
5				If $\mathrm{rand}\le P_C$ or $d==d_r$ then
6					$H_i\left(d\right)=H_1\left(d\right)+F·\left(H_2\left(d\right)-H_3\left(d\right)\right)$
7				Else
8					$H_i\left(d\right)=H_{\min }\left(d\right)+\mathrm{rand}·(H_{\max }\left(d\right)-H_{\min }\left(d\right)$
9				End If
10			Else
11				If $\mathrm{rand}\le P_C$ or $d==d_r$ then
12					$H_i\left(d\right)=H_1\left(d\right)$
13				Else
14					$H_i\left(d\right)=H_{\mathrm{best}}\left(d\right)+F·\left(H_2\left(d\right)-H_3\left(d\right)\right)$
15				End If
16			End If
17		End For
18	End For

3. Proposed Approach for BBO

3.1. Operators, Modifications and Hybridization

Initially, only the concepts presented by Dan Simon [5], such as the linear migration model and basic mutation, were used in some very basic example functions, to carry out an initial study of the method. Later, the migration models presented here by Ma [7], the taxonomic mutation operator by Hao [10], and the concept of blending in migration by Ma and Simon [9] were added. This completed the initial standard BBO algorithm, which was then applied to the first three cases presented below, allowing the algorithm's behavior to be studied further.

Once this phase had been completed, the algorithm was modified to be able to remedy some of its shortcomings. To do this, a mutation cycle, in which the mutation operators are used several times in the same generation, only allows for the addition of improved solutions. A solution injection process was also added, where the last 10 % of the population was replaced by new solutions at the end of each generation. All these processes were added in an initial iteration of the studies to improve the diversity of the algorithm's population.

Finally, in the hybridization section, two concepts from the genetic and memetic algorithms were implemented to improve the method's performance. Initially, binary crossover operators were added to improve the exchange of genetic information, something that had already been explored in the BBO world with the concept of blending [9]. For the implementation of these operators, it was also necessary to implement a binary decoder and encoder since most of the operators and heuristics of the BBO method are in the real domain. Recombination processes are based on the selection of a couple for each child that will be generated. Three operators, the single-point crossover (SP) [11], multi-point crossover (MP) [11], and the uniform crossover (U) [12], which were subjected to two different mating selection mechanisms (MSM).

The first strategy for MSM denoted by Previous Elites-based Selection (PEbS), is proceeded as follows:

The first parent is selected from the current population with a probability equal to the migration rate (proportional to rank position). The second parent is selected randomly from the best solutions stored in previous generations. Then the offspring created by crossover will be evaluated. The new population will be composed using the best solutions from the current population and the offspring. The enlarged population is sorted by fitness and the worst individuals are eliminated to restore the original size of the population.

The second strategy for MSM denoted by Independent Selection-based Fitness (ISbF), is proceeded as follows:

The current population is divided into two groups according to the fitness of each solution, the best (elite) group and the worst group of solutions [13]. Each group provides a parent, which is selected through a roulette method (equation 3.1). The two stochastic procedures are independent. Then the offspring created will all be placed in the population, replacing every solution not belonging to the elite.

$$P\left({\mathit{x}}_j^k\right)={\displaystyle \frac{f\left({\mathit{x}}_{\mathit{j}}^{\mathit{k}}\right)}{{\sum }_{j=1}^{N_{\mathrm{pop}}}f\left({\mathit{x}}_j^k\right)}}$$

(3.1)

where $f\left({\mathit{x}}_{\mathit{j}}^{\mathit{k}}\right)$ represents the fitness value of the solution and $N_{\mathrm{pop}}$ the size of the population.

To conclude the hybridization, it was also added a Hooke-Jeeves algorithm as a search optimizer, to improve the final generations’ local search, something that has been well documented as a flaw of the BBO method. This operator will allow for an improvement of the search velocity, allowing for the algorithm to reach the end solution with greater ease. The Hooke-Jeeves method is common in the world of optimization, since it allows for an efficient search for spaces of great dimensions, even if the objective function is very complex. This method and its parameters were altered, so that its performance and synergy with the other operators was greatly improved. The following algorithm (Algorithm 4) shows the pseudocode for the local search method.

Algorithm 4: Hook-Jeeves local search method pseudocode
1	Compute the step size ($h^1$)
2	For $i=1$ to $N_{\mathrm{pop}}$ do
3		Select point (solution) to be explored, $H_i$
4		Perform the first step iteration (equation (3.2))
5		Evaluate the new solution, $H_{\mathrm{new}}$ in BBO algorithm
6		If $H_{\mathrm{new}}$ is better than $H_i$ then
7			$H_i$ $\leftarrow H_{\mathrm{new}}$
8		End If
9		Save the solution information, $H_{\mathrm{old}}$
10		/* Local search loop */
11		While stopping criteria are not verified do
12			Compute the new step value (equation (3.3))
13			Compute the new solution $H_{\mathrm{new}}$
14			If $H_{\mathrm{new}}$ is better than $H_{\mathrm{old}}$ then
15				$H_{\mathrm{old}}\leftarrow H_{\mathrm{new}}$
16				If $H_{\mathrm{new}}$ better than $H_i$ then
17					$H_i\leftarrow H_{\mathrm{new}}$
18				End If
19			End If
20		End While
21	End For

The first step, $h^1$, is calculated as follows,

$$h^1\left(j\right)=H_1\left(j\right)-H_{N_{\mathrm{pop}}}\left(j\right)$$

(3.2)

where $j$ represents the dimension, $H_1$ is the best solution from current population and $H_{N_{\mathrm{pop}}}$ the worst solution. The step reduction is given by,

$$h^i\left(j\right)=\left({\displaystyle \frac{1}{i}}\right)·h^1\left(j\right)$$

(3.3)

where $i$ represents the number of cycles of the local search loop. It is important to mention that the steps will be subtracted and added so that a complete search can be carried out.

For the practical cases presented here, the algorithm was tested using the migration operator in the standard and blending versions and using the mutation operator in the standard and taxonomic versions. These operators were then tested with a local optimizer and a crossover operator. First, each approach separately and after using a hybrid procedure. The studies focused on the best solutions obtained and the evolution of the elite average over the generations.

3.2. Optimization Algorithm

The algorithm was implemented in FORTRAN, with the intel compiler for visual studio. Initially, the algorithm randomly creates the initial population and calculates the mutation and migration rates. After that, the algorithm will evaluate the solutions and rank them according to their fitness. Then the solutions will go through the migration process, where it was defined that the top two solutions will not receive information from other solutions, only share. Next, the algorithm moves on to the mutation operators, where all solutions can be mutated, but due to the mutation cycle, the newly generated solutions will only enter the population if they represent an improvement. At each one hundred generations, the solutions will be subjected to the crossover operators. After this procedure, the algorithm enters the local search phase. Here, the new population is classified, so that at the end of this phase the solution injection process can replace the worst solutions. Figure 4 shows the flow diagram of the hybridized algorithm.

Figure 6. Flow diagram of the hybridized BBO [14].

Figure 6. Flow diagram of the hybridized BBO [14].

In terms of simulation parameters, it was used 50 initial solutions, 10% of which will represent the evaluated elite, an emigration ($E$) and immigration ($I$) values equal to one, a maximum mutation ($m_{\max }$) of 0.01 and an initial number of mutation cycles of three, increased by one after every 100 generations, until to maximum value of 20.

For the Hook-Jeeves local search, it was used a constant selection probability equal to 0.1 to select solutions from BBO population, except for the two best solutions, which were always selected. The stopping criteria of Hook-Jeeves local search is based on the number of cycles, set at one hundred.

4. Case Studies

Case Study 1: Welded Beam Design

The welded beam, shown in Figure 7, is a classical optimization problem collected from [15], where its objective is to minimize the cost of the welded joint. This is a non-linear problem, not only in its objective function but also in its constraints.

From a geometrical standpoint, the objective function can be given as,

$$f\left(x_1,x_2,x_3,x_4\right)=\left(C_1+C_3\right)·V_{welding}+C_2·V_{\mathrm{beam}}$$

(4.1)

where $\left(x_1,x_2,x_3,x_4\right)=\left(h,l,t,b\right)$ and where $C_1=0.10471,C_2=0.04811,C_3=1$ and,

$$V_{\mathrm{welding}}=x_1^2·x_2$$

(4.2)

$$V_{beam}=x_3·x_4$$

(4.3)

The constraints imposed in the problem are,

$$g_1\left(\mathit{x}\right)=\tau \left(\mathit{x}\right)-{\tau }_{\max }\le 0$$

(4.4)

$$g_2\left(\mathit{x}\right)=\sigma \left(\mathit{x}\right)-{\sigma }_{\max }\le 0$$

（4.5）

$$g_3\left(\mathit{x}\right)=x_1-x_4$$

(4.6)

$$g_4\left(\mathit{x}\right)=C_1·x_1^2+C_2·x_3·x_4·\left(L+x_2\right)-5\le 0$$

(4.7)

$$g_5\left(\mathit{x}\right)=0.125-x_1\le 0$$

(4.8)

$$g_6\left(\mathit{x}\right)=\delta \left(\mathit{x}\right)-{\delta }_{\max }\le 0$$

(4.9)

$$g_7\left(\mathit{x}\right)=P-P_c\left(\mathit{x}\right)\le 0$$

(4.10)

$$g_8\mathrm{to}{g}_{11}:0.1\le x_i\le 2,i=1,4$$

(4.11)

$$g_{12}\mathrm{to}{g}_{15}:0.1\le x_i\le 10,i=\mathrm{2,3}$$

(4.12)

where,

$$\tau \left(\mathit{x}\right)=\sqrt{{\left({\tau }^{\prime }\right)}^2+2·{\tau }^{\prime }·{\tau }^{″}{\displaystyle \frac{x_2}{2R}}+{\left(t^{″}\right)}^2}$$

(4.13)

$${\tau }^{\prime }={\displaystyle \frac{P}{\sqrt{2}·x_1{x}_2}}$$

(4.14)

$${\tau }^{″}={\displaystyle \frac{M·R}{J}}$$

(4.15)

$$M=P\left(L+{\displaystyle \frac{x_2}{2}}\right)$$

(4.16)

$$R=\sqrt{{\displaystyle \frac{x_2^2}{4}}+{\left({\displaystyle \frac{x_1+x_3}{2}}\right)}^2}$$

(4.17)

$$J=2\left[x_1{x}_2\sqrt{2}\left({\displaystyle \frac{x_2^2}{12}}+{\left({\displaystyle \frac{x_1+x_3}{2}}\right)}^2\right)\right]$$

(4.18)

$$\sigma \left(\mathit{x}\right)={\displaystyle \frac{6PL}{x_4{x}_3^2}}$$

(4.19)

$$\delta \left(\mathit{x}\right)={\displaystyle \frac{4P{L}^3}{E{x}_4{x}_3^2}}$$

(4.20)

$$P_c\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{4,013}E·\sqrt{\frac{x_3^2{x}_4^6}{36}}}{L^2}}\left(1-{\displaystyle \frac{x_3}{2L}}·\sqrt{{\displaystyle \frac{E}{4G}}}\right)$$

(4.21)

where $P=6000\mathrm{lb},E=30·{10}^6\mathrm{psi},L=14\mathrm{in},G=12·{10}^6\mathrm{psi},$ $\sigma =30000\mathrm{psi},{\tau }_{\max }=13600\mathrm{psi},{\delta }_{\max }=0.25.$

For this problem, the fitness function is defined as,

$$\mathrm{HSI}\left(\mathit{x}\right)=C-\left(f\left(\mathit{x}\right)+\lambda ·\mathsf{\Gamma }\left(\mathit{x}\right)\right)$$

(4.22)

where $f\left(\mathit{x}\right)$ is the objective function $\mathsf{\Gamma }\left(\mathit{x}\right)$ the penalty function and $C$ and $\lambda$ are two constants, respectively equal to 10 and 100. The penalty function is then given as,

$$\mathsf{\Gamma }\left(\mathit{x}\right)=\sum _{i=1}^7{u}_i{\left(\mathit{x}\right)}^n$$

(4.23)

where $n=2$ and with $u_i\left(\mathit{x}\right)$ defined as,

$$u_i\left(\mathit{x}\right)=\left\{\begin{array}{c}{g}_i\left(\mathit{x}\right)\mathrm{i}\mathrm{f}{g}_i\left(\mathit{x}\right)>0\\ 0\mathrm{i}\mathrm{f}{g}_i\left(\mathit{x}\right)\le 0\end{array}\right.$$

(4.24)

The optimal solution of [16] will be used as a reference optimal solution, which will be compared with the results obtained using the hybrid BBO method.

$${\mathit{x}}^{\mathit{*}}=\left(0.20572963,3.4704886,9.0366239,0.20572963\right)$$

(4.25)

$$f\left({\mathit{x}}^{*}\right)=1.7248523085$$

(4.26)

Results for the Case Study 1

This case was used to study and compare the performance of the algorithm and its different operators. The results were obtained after the program reached its stopping criteria, the number of generations, which was set at two thousand generations.

Starting with the analysis of the best solutions obtained, Table 1 shows some of the best objective results for the different tests. It is possible to note that the fully hybridized algorithm, with both mating selection mechanisms of recombination, PEbS and ISbF, and the algorithm with local search can improve the optimal solution, reaching a value lower than the one used as a reference [17]. The worst results were obtained using taxonomic mutation without blending, which was expected due to the lack of diversity in the population. Indeed, since both taxonomic mutation and unblended migration use information from the rest of the population to operate, and without a guaranteed source of greater diversity, solutions tend to have a hard time leaving local optima. Furthermore, it is important to note that the impact of the crossover operator in this case was not very important, since both solutions presented here have a value similar to that obtained in the tests without crossover. In terms of migration models, the linear model was the best, having consistently reached good solutions. The other models also performed well, especially when using other operators from the hybridization.

The evolution of the elite’s average fitness is represented, graphically, in Figure 8 and Figure 9. Highlighting Figure 8, it is possible to note one of the best benefits of the method, its fast global search, with the steep drop of the elite’s average in the first stages of the program. But that graph also highlights the biggest weakness of the standard BBO algorithm, its behavior at final generations. As it is possible to see, the algorithm spent a lot of generations stuck in a single, not optimal, solution. A problem that was somewhat solved by the use of the local search method, as it is possible to note in Figure 9, the curve here presents a more common curve, that progressively evolves into the best solution, while still preserving some of the qualities from the original method, such as the steep drop in the initial generations.

In summary, the following conclusions can be drawn:

-

Fully hybridized model shows the best results. It is even possible to obtain better results than the reference value [17];

-

The Hooke-Jeeves local optimizer was the inclusion with the most impact. In this example, the influence of the crossover was not very significant;

-

Standard mutation had the best performance, with and without the use of the blending strategy in migration;

-

Blending had an impact, also showing good synergy with the local optimizer and crossover operators;

-

The linear migration model obtained the best performance, proving to be more consistent.

-

It should be noted that the sinusoidal migration model also obtained good results.

Case Study 2: Welded Gusset Optimal Design

A welded Gusset with the geometry as shown in Figure 10 is considered in this problem [15]. A profile $L$ $100\times 65\times 11\left[\mathrm{mm}\right]$ with a section of $A=1.7\times {10}^{-3}{\mathrm{m}}^2$, must be welded to a $1.4\times {10}^{-3}\mathrm{m}$ thick plate so that the section of the profile can be subjected to tension in its totality. Considering that the base material is St37 steel (DIN 17100) with ${\sigma }_{\mathrm{adm}}\le 137.34\mathrm{MPa}$ and the metal deposited in the fillet welds has ${\sigma }_{\mathrm{adm},\mathrm{weld}}\le 107.91\mathrm{MPa}$.

The welding cost includes the cost of preparing the joint to be welded, the cost of the material deposited, and the cost associated with the welding technology. Supposing that the cost can be expressed as a function of the welded joint dimensions, it is obtained the following objective function.

$$\mathrm{Cos}\mathrm{t}=K\left(l_1·a_1+l_2·a_2\right)$$

(4.27)

where $K$ is associated with the price of the manpower and materials cost, $l_1\mathrm{and}{l}_2$ are the lengths of the fillet welds, $a_1\mathrm{and}{a}_2$ are the respective widths of the projected weld throats according to the rules of the welding code design. In Figure 10 let’s consider the following data: $b=0.1\mathrm{m}$; $e_1=0.034\mathrm{m}$ and $e_x=0.0237\mathrm{m}$ .

To evaluate the normal stress caused by the bending moment it is necessary to calculate the flexural strength module, $Z_f$,

$$Z_f=Z_{f1}+Z_{f2}$$

(4.28)

$$Z_{f1}={\displaystyle \frac{1}{6}}{l}_1^2{a}_1$$

(4.29)

$$Z_{f2}={\displaystyle \frac{1}{6}}{l}_2^2{a}_2$$

(4.30)

The normal stresses caused by the bending moment is defined as follows,

$$n={\displaystyle \frac{M_{\max }}{Z_f}}$$

(4.31)

where $M_{\max }=P\times e_x=5533.428\mathrm{Nm}$.

The stress in the fillet welds (parallel to its longitudinal axis) is given as follows,

$$t_{//}={\displaystyle \frac{P}{l_1{a}_1+l_2{a}_2}}$$

(4.32)

where $P=A\times {\sigma }_{\mathrm{adm}}=233.47\mathrm{kN}$.

Finally, according to the Portuguese regulation code on steel structures for buildings (REAE) [18], the equivalent stress ${\sigma }_{\mathrm{eq}}$, is given by:

$${\sigma }_{\mathrm{eq}}=\sqrt{1.4n^2+1.8t_{//}^2}$$

(4.33)

The stress constraint of the welded Gusset design problem is given by,

$${\sigma }_{eq}\le \overline{\alpha }{\sigma }_{adm,weld}$$

(4.34)

where $\overline{\alpha }$ is a safety factor defined according to REAE (artº77) [18] successively defined as follows,

$${\alpha }_i=\mathrm{0,8}\left(1+{\displaystyle \frac{1}{a_i}}\right),i=\mathrm{1,2}$$

(4.35)

with

$$a_i\ge 3\left[\mathrm{mm}\right],i=\mathrm{1,2}$$

(4.36)

$$\overline{\alpha }=\min \left({\alpha }_i,i=\mathrm{1,2}\right)$$

(4.37)

The fitness function for this example is going to be calculated using the formulation from equation (4.22) but using $C={10}^4,\lambda ={10}^{16}\mathrm{and}n=2$, so that no obtained solutions violate the constraints. The optimal reference point, obtained from [17], is:

$${\mathit{x}}^{*}=\left(5.13333,5.13333,400.0,248.196\right)$$

(4.38)

$$f\left({\mathit{x}}^{*}\right)=3228.13$$

(4.39)

Results for the Case Study 2

The second example used the same stopping criteria as the first example. In terms of the best solutions obtained, presented in Table 2, it is possible to note, once again, that the best solutions are improved when compared to the reference value. In this case, the use of the fully hybridized algorithm proved to be superior to the other approaches, as it repeatedly achieved the best solutions, even when using different crossover operators. The study is also a particularly good example of the benefits of adding crossover operators, since in both cases, with and without local search, their implementation helped to improve the final solution. The local search, once again, also proved its value by significantly improving the solution. On the downside, when the blending in migration was not used, the solution fell short of all the other cases, even with the help of the local optimizer. This is more than likely due to the non-blended migration since it does not introduce any diversity to the population, by only replacing the bad solutions with the better ones. This feature, when combined with a complex function, can lead to a case, like this one, where the algorithm was led to a local optimum and was never able to create a new solution outside of that point, which led to it being stuck at that point for several generations. This aspect is one of the biggest weaknesses of standard BBO algorithms, which had already been documented by other authors, referred to as a lack of speed. Dan Simon [5], himself mentioned this lack of speed in his original paper. Hopefully, the concept of blending solves this problem, even if partially.

Using the evolution of the elite's average fitness, Figure 11 once again shows the rapid decline in the first few generations of the algorithm. This figure, when compared with the previous figures, is also a good example of some of the problems that have occurred with the trapezoidal model, namely its inconsistency. This migration model was either one of the fastest or the slowest, the latter being more frequent. It also doesn't help that this model also shows poor synergy with the other modifications, specifically with taxonomic mutation. Figure 12, demonstrates a little bit of the influence of the crossover operator, with the curve of the graph being completely different from what is expected from this method. The conclusions are:

-

The fully hybridized model with the different options originally not included in the BBO continued to achieve the best results;

-

The inclusion of the crossover operator proved to be a good option, showing good results both with and without the use of the Hooke-Jeeves local optimizer. Both strategies achieved good results;

-

The Hooke-Jeeves local optimizer continued to be an advantage, showing several robust results;

-

The migration model with blending was an asset. In this example, we can see its impact on synergy with the other hybridization options, especially the Hooke-Jeeves local optimizer;

-

The linear migration model was the most prevalent in this case;

-

Both mutations achieved good results, but the taxonomic mutation showed worse synergy when integrated with other options, slowing down the algorithm in some cases.

Case Study 3: Design of Springs and Weights System

Figure 13 shows a spring system supporting weights at the connections between the springs [19]. First, the system analysis is performed to determine the equilibrium position by minimizing the potential energy. For this, five weights and six springs were considered.

The $i$-th spring deformation is given by the following equation,

$$\mathsf{\Delta }{L}_i={\left[{\left(x_{i+1}-x_i\right)}^2+{\left(y_{i+1}-y_i\right)}^2\right]}^{{\displaystyle \frac{1}{2}}}-L_i^0$$

(4.40)

with $L_i^0=10\mathrm{m}$.

There are $N+1$ springs, with $N$ being the number of weights ($N=5$ for this case). The stiffness of the $i$-th spring is,

$$K_i=500+200{\left({\displaystyle \frac{N}{3}}-i\right)}^2\left[\mathrm{N}/\mathrm{m}\right]$$

(4.41)

and the weight applied at the $j$-th node is given as follows,

$$W_j=50j\left[\mathrm{kg}\right]$$

(4.42)

The potential energy of the system is then given by,

$$PE={\sum }_{i=1}^{N+1}{\displaystyle \frac{1}{2}}{K}_i{\left(\Delta L_i\right)}^2+{\sum }_{j=1}^K{W}_j·y_{j+1}$$

(4.43)

The design variables are the nodes position in the spring-deformed state. Therefore, the design variables are $x_i,i=2,\dots ,6$ and $y_i,i=2,\dots ,6$, for a total of ten design variables. The fitness function used is going to be equation (4.43) and the only constraints applied to the problem will be the size constraints, so that the solution obtained ensures that no nodes interfere with each other in a physically impossible manner. The optimal solution described in the book [19] is given as,

$$PE\left({\mathit{x}}^{\mathit{*}}\right)=-4416\mathrm{J}$$

(4.44)

Results for the Case Study 3

Analyzing the best solutions obtained, shown in Table 3, it is important to note that only the fully hybridized method was able to reach the reference solution. The worst solution came, once again, from the tests without blending and with taxonomic mutation, assuming that the same reasons that plagued the two previous problems also hold true here. The addition of the crossover operators was, once again, very beneficial, their addition allowed for an improvement of solutions. Furthermore, both strategies performed well and were able to provide valuable solutions equally. The linear model of migration also seems to continue its dominance having a lot of placements, but this time the best solutions came from different models.

In terms of the evolution of the elite’s average fitness with the generations, Figure 14 shows more of the same behavior. The steep initial drop and the intermediate and final slump, which remain in a non-optimal solution. Figure 15 on the other hand showed a behavior similar to other GAs, this mainly shows the influence that the crossover operators had in this example, demonstrating, once again, its benefits.

However, this example highlighted a problem that affected the algorithm's performance, slowing it down. This error arose with the use of taxonomic mutation together crossover operators. In fact, in the initial generations, due to the generation of bad solutions in the recombination process, the performance of the differential taxonomic mutation worsened significantly. Using these bad solutions from the crossover worsens the quality of the population, thus reducing the convergence rate. However, this problem can be easily solved by using a solution selection process after the generation of offspring during the crossover, preventing very bad solutions from entering the new population. Or adopt a hybrid mutation strategy based on the BBO method's standard mutation in the initial generations and taxonomic mutation in the intermediate and final generations. The conclusions are:

-

The algorithm presented more difficulties in this example, but managed to obtain solutions close to the reference value when fully hybridized;

-

Both options included (crossover and local optimizer) had an impact, with the local optimizer being the most important;

-

The linear model continued to have the best performance, but the other non-linear models obtained the best results in the end;

-

Standard mutation obtained the best results, due to a problem that differential mutation had with the bad solutions coming from the crossover operators;

-

Migration with blending has continued to be an impactful modification, and this example demonstrates it well.

Case Study 4: The 10-bar Truss

The objective of this problem is to minimize the mass of the structure by varying the cross-sectional area of the elements while respecting the limitations and restrictions based on stresses and displacements. Figure 16 shows the geometry that is going to be studied in this example from [20].

The structure is loaded on nodes 2 and 4, by two vertical loads (one per node) of $P_y=-4.45\times {10}^5\mathrm{N}$. The Young’s modulus of the material is $E=6.9\times {10}^8\mathrm{N}/{\mathrm{m}}^2$.

The objective function of the problem is,

$$f\left(x\right)=\rho \sum _{i=1}^{10}{A}_i{l}_i$$

(4.45)

where $\rho =2.77\times {10}^3\mathrm{kg}/{\mathrm{m}}^3$ represents the density, $l_i$ the length of the element and $A_i$ the cross-section area of the element.

The constraint values for the stresses and displacements are $\left|\sigma \right|\le 1.72\times {10}^8\mathrm{N}/{\mathrm{m}}^2$ (traction /compression) and $\left|\delta \right|\le 5.08\times {10}^{-2}\mathrm{m}$, respectively. The cross-sectional areas of members must belong to the interval from $6.45\times {10}^{-5}{\mathrm{m}}^2$ to $2.26\times {10}^{-2}{\mathrm{m}}^2$.

The reference value from [20], which was obtained using an AGA (Augmented Genetic Algorithm), is presented below,

$$f\left({\mathit{x}}^{*}\right)=2261.42\left[\mathrm{kg}\right]$$

(4.46)

It is also important to note that the linear equation system was solved using gaussian elimination with partial pivoting.

Results for the Case Study 4

The best solutions presented in Table 4, show that in the majority of the tests, the algorithm did not reach a value close to the reference value, this may be due to the combination of the complexity of the problem and the limitation created by the stopping criteria, that, once again, was placed as 5000 generations. But on a positive note, there were also many results better than the reference value. In terms of migration models, it is important to note the good synergy that the linear model continues to show, having great results with all of the modifications. The sinusoidal model, which obtained the best value, and the trapezoidal model, which consistently obtained the best values for the respective test, should also be highlighted. This example also proved that with a more complex migration model, it is possible to obtain better results. On the side of the mutation operator, both operators (standard and taxonomic) were able to reach good values, but it is important to highlight that when the local search was used, the standard mutation is better. This factor, combined with the fact that linear migration is the model that behaves best with taxonomic mutation, shows some of the synergy problems that have already been exposed here. Finally, the crossover operator also proved to be a positive influence in this case, allowing for the improvement of the solutions.

Analyzing the evolution of the elite’s average fitness, presented in Figure 17 and Figure 18, it is possible to witness the same phenomenon mentioned before. Highlighting the first graph, it is possible to see the influence of the local search, with the steady increase in value on the trapezoidal model curve. In addition, this graph shows the dominance of the sinusoidal model, which, without that huge spike in the value of the trapezoidal model, would be the best in these tests. The second graph also shows some interesting behaviors, especially the uniform crossover (U) that shows a very steady increase in value elite’s average fitness through the generations. The curves relating to the use of the crossover operator show better algorithmic performance.

-

The inclusion options (crossover and local optimizer) had quite an impact. The fully hybridized algorithm obtained the best results;

-

This example has shown that for more complex problems, more complex migration models obtain better results. However, the linear model manages to overcome this advantage in some cases with its good synergy with other techniques;

-

The model, without hybridization, struggled to achieve reasonable results;

-

The migration model with blending maintained its good performance, especially during tests including other hybridization options;

-

The model with standard mutation obtained the best results, but taxonomic mutation also performed well.

5. Conclusions

An evolutionary algorithm based on biogeography (BBO) was implemented to be applied to various problems in the topic of optimization of mechanical systems and to study their behavior. This work adapted not only the linear migration model from Dan Simon's original work [5] but also the trapezoidal, quadratic, and sinusoidal models presented in Ma's work [7]. In terms of mutation operators, the standard BBO mutation and Hao's taxonomic mutation were implemented [10]. The initial algorithm was later hybridized with other concepts and heuristics from other methods, such as the crossover operators Single-Point (SP), Multiple-Point (MP), and Uniform (U) used in genetic algorithms, and local search-based learning from memetic algorithm’s idea.

The standard non-hybridized BBO showed its best performance in the initial generations, due to its good global search. However, the results obtained fell short of the benchmarks in the literature using other optimization methods. This poor performance was due to the lack of diversity of solutions in the populations at the final generations of the optimization process. This feature is associated with the replacement of solutions of the migration operator an intrinsic property of BBO method. The addition of the blending strategy, applied during the migration process, introduced by Dan Simon the original author of BBO, proved to be very beneficial, improving the algorithm's performance in several respects. Even so, some of BBO's initial problems persisted. Thus, other processes, such as the injection of new solutions [21], [13] before BBO hybridization, have proved effective in improving population diversity. The results obtained after these inclusions proved to be positive.

In terms of migration, the models that stood out the most during the tests were the sinusoidal and linear models, which managed to produce several optimal results. The quadratic model, although less consistent, also managed to produce satisfactory results, outperforming the other models in some cases. The worst model, the trapezoidal model, showed reliable results but was the slowest and most inconsistent when compared to the other models. In addition, this model also showed poor synergy with the other operators implemented, such as taxonomic mutation and the local optimizer. The standard mutation and taxonomic mutation operators were able to produce good results when migration with blending was used. Standard mutation prioritized the speed of the algorithm, generating good solutions quickly, while taxonomic mutation showed greater precision, at the cost of some speed.

The hybridization process began with the implementation of the crossover operators, which required the implementation of binary coding of solutions. The three implemented crossover operators were subject to two different mating selection strategies. The inclusion of these operators proved to be positive, with an improvement in results, in general. The second mating selection strategy (ISbF) appeared to have better behavior, especially when used in conjunction with the MP (Multiple-Point) crossover operator, but the first mating selection strategy (PEbS) was also quite positive. The three crossover operators obtained satisfactory, fairly similar results and showed good synergy with the other operators and concepts, especially with blending migration.

Then, the Hooke-Jeeves local optimizer was implemented, which presented excellent results, even in the tests without the inclusion of the crossover. However, the best results were obtained when the fully hybridized algorithm was used. Thus, the algorithm was able not only to quickly reach the optimum points but also, in some cases, to surpass them.

When the algorithm was exposed to a more complex problem, it presented some difficulties, initially without the use of hybridization. However, with the inclusion of crossover operators and local search, it managed to achieve good values, albeit with some struggle. This example also showed that complex migration models have an easier time dealing with complex problems. But it also showed that the greatest advantage of the linear migration model, its great synergy with the other modifications, allowed it to compete with the more complex models.

In conclusion, the initial behavior of the BBO original method suggested that it needed some improvements to reach its true potential. With the proposal innovations, the method showed excellent results, achieving its objective with relative ease. It can be said that these innovations complement the method positively and that their inclusion should be further strengthened.