1. Introduction
Structural optimization is considered as one of the most important fields in engineering because it can propose innovative designs that satisfy different multidiscipline design requirements. There are three kinds of structural optimization: geometry, size, and topology. Geometry optimization involves finding the best arrangement of the joints. In size optimization, the main goal is to find the dimensions of the elements where the locations of the joints are fixed. Structural topology optimization is used to determine the best number, location, and size of voids as well as the best connectivity of the design domain. For any structural topology optimization problem, the applied loads, the support condition, the possible design requirements or any other restrictions, and the volume of structure to be constructed are the only known quantities. Then, in an iterative manner, the material is redistributed to meet the design requirements [
1].
Structural topology optimization methods have received much attention in the relevant literature. They have been improved because they are very attractive for use in many industrial applications in civil and mechanical engineering [
2]. These methods can help designers select suitable structures and reduce the effort required by trial and error. In structural topology optimization, similar to any other type of optimization, two methods are used: the deterministic (i.e., calculus-based or gradient-based) method and stochastic (i.e., metaheuristic or nature-inspired) algorithms.
In the deterministic method, the Solid Isotopic Microstructure with Penalization (SIMP) or the power law method introduced by Bendsøe [
1] is a generally accepted and efficient tool for structural topology optimization. In the SIMP method, the design variables are the density of each element, and the material property of each element is constant [
3]. The gradient-based SIMP method is used to search for the optimum point in the neighborhood of the current estimate. It tends to converge to a local optimum with a blurry boundary. Moreover, the SIMP method cannot directly solve 0/1 problems, so it requires a filtering method to resolve the final design [
4,
5,
6]. The Canonical Duality Theory (CDT) developed by Gao [
7] is used to analytically solve structural topology optimization problems with solid or void material distribution. The numerical results of solving benchmark topology optimization problems showed that the CDT method is powerful in comparison with SIMP and Basic Evolutionary Structural Optimization (BESO) methods in terms of computational time, the final compliance, not generating grayscale areas, and fewer checkerboard patterns [
8,
9,
10]. For problems with mesh sizes of 180 × 60 and 150 × 50, the CDT method showed some checkerboard patterns. Rational Approximation of Material Properties (RAMP) developed by Stolpe and Svanberg [
11] and the generalized form of it (called GRAMP) were used to solve two benchmark topology optimization problems for minimum compliance. Using a penalization factor of value equal to 3, Dzierżanowski showed that GRAMP was able to produce an almost “0–1” layout of material [
12]. In comparison with the SIMP method, GRAMP and RAMP were able to obtain designs with fewer grayscale areas and better compliance. However, GRAMP and RAMP were not able to produce exactly void/solid material distribution. In addition to the compliance minimization topology design methods, the stress-constrained topology optimization methods have received much attention since the pioneering work of Duysinx and Bendsøe [
13]. It involves stress-related objectives or constraints to overcome three main issues in the minimizing compliance design: highly non-linear stress behavior, the local nature of stress constraints, and the stress singularity [
14,
15,
16]. Similar to the compliance topology optimization approaches, the stress-constrained topology optimization methods search locally, and final designs can have blurry boundaries or checkerboard patterns depending on an exponent factor.
The deterministic approach described in the preceding paragraph requires gradient information of functions. The other fundamental approach to topology optimization is to use metaheuristic methods [
17]. Over the last few decades, many metaheuristic optimization algorithms have been explored to solve complex optimization problems. These methods are based on simulations, and they do not require gradient information. Based on the number of solutions in each iteration, there are two types of metaheuristic: (1) single solution-based, such as the Simulating Annealing (SA) and Harmony Search (HS) algorithms, and (2) population-based, such as the Genetic Algorithm (GA). Based on the inspiration, metaheuristic algorithms can be categorized into four types: evolutionary algorithms, such as GA and the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [
18,
19], physical-based, such as Colliding Bodies (CB), human-based, such as Tabu Search (TS), and swarm intelligence, such as Particle Swarm Optimization (PSO) [
20]. All types of metaheuristic algorithms have common features: the search in the entire design space, two searching phases—exploration and exploitation, the use of random numbers to generate new designs, and no good initial designs are required. Metaheuristic algorithms are suitable for tackling many real-life optimization problems, especially those having many local minimum designs or those in which it is not possible to calculate the gradient information of the problem. Like gradient-based algorithms, metaheuristic algorithms have some drawbacks. That is, depending on the number of design variables and the range of each design, they may require many function evaluations to obtain a reasonable solution in comparison with gradient-based optimization algorithms because the number of possible combinations drastically increases. However, this computation efficiency comparison may not be meaningful because the gradient-based approach is a local search method, and the problem of global optimization requires all local optimal solutions. Thus, the deterministic approach may become as computationally expensive as metaheuristic algorithms are [
4].
The Bat Algorithm (BA) is based on observations of a natural phenomenon. It was inspired by the echo locative behavior of bats. In the BA, artificial bats (
designs of dimension
) are search agents mimicking the search nature of real bats using loudness and emission rates. Regarding some benchmark continuous design variable optimization problems, such as Rosenbrok’s (
) and Schwefel’s (
), it was proven that the BA is capable of obtaining competitive results in comparison with other well-known metaheuristic algorithms, such as GA and PSO [
21]. The BA has also been used to solve some discrete design variables optimization problems, such as the symmetric and asymmetric traveling salesman problems [
22] and steel frames subjected to American Institute of Steel Construction-Allowable Stress Design (AISC-ASD) design requirements [
23]. The binary version of the BA was developed by Mirjalili et al. [
24], who termed it the Binary Bat Algorithm (BBA). The BBA has been used to solve 22 benchmark functions. The results have shown that the BBA was able to surpass the binary PSO and GA in terms of the quality of the solutions and the total number of function evaluations required to reach the optimal design [
24].
Metaheuristic algorithms for topology optimization usually use binary code (i.e., 0/1) to represent void/solid material distribution. Thus, they can be used directly to solve structural topology optimization problems because no gray area is generated. Although a great advantage over gradient-based topology optimization algorithms, in addition to the advantages mentioned above, metaheuristic algorithms are generally less efficient than mathematical programming methods [
4,
5,
11]. GA has been used to solve topology optimization problems for more than three decades [
5,
11,
12,
13,
14]. Kaveh et al. [
15] used the Ant Colony methodology for the topology optimization of structures to find the stiffest structures based on the element’s contribution to the strain energy. Wu and Tseng [
3] modified the binary differential evolution to solve structural topology optimization problems. The finest finite element mesh size used was
(
elements). Azar et al. [
16] developed binary particle swarm optimization (BPSO) to solve topology optimization problems when finite element mesh size was
(
elements). Lee and Han [
25] proposed a topology optimization scheme based on HS to investigate static stiffness topology optimization problems. A finite element mesh size of
(
elements) was used.
This study is conducted because most of the previous work in implementing metaheuristic algorithms to solve topology optimization problems was limited to a few hundred design variables. Many real-life structural topology optimization problems require discretizing the design domain into a fine mesh to make sure the finite element method yields accurate results. The present study is an attempt to improve the quality of the solution of the topology optimization method using a large number of design variables (i.e., fine finite element mesh) by incorporating a new filtering algorithm and penalty function to suit the nature of the bit-array topology optimization. In this study, the numerical results of the proposed algorithm used to solve the benchmark structural topology optimization problems obtained are compared with those found using the SIMP method. The comparison is made based on minimum compliances and final design configurations to demonstrate the performance and viability of the proposed algorithm. The rest of this paper is organized as follows:
Section 2 describes the problem formulation for structural topology optimization. In
Section 3, the BA and the BBA are reviewed briefly. Structural topology optimization using the BBA is discussed in
Section 4. The numerical and design results are presented in
Section 5. Finally, the conclusions of the study are described in
Section 6.
2. Topology Optimization: Problem Formulation
The main aim of structural topology optimization is to maximize the stiffness of a structure subjected to a volume constraint or to minimize the material used to construct a structure subject to stress and/or deflection constraints. Specifically, the optimization problem is to minimize an objective function subjected to some constraints by gradually distrusting the material or removing the low stressed material in a predefined domain [
26]. In general, the structural topology optimization problem can be stated as follows:
where
is an objective function (e.g., strain energy),
is a volume constraint,
is a volume limit,
are other possible inequality constraints (e.g., displacement),
are other possible equality constraints,
is the density variable that describes the material distribution in the design domain (
), and
is the vector of all points in the design domain [
27]. The reader is referred to References [
28,
29,
30] for more details.
For topology optimization using the BBA (or any metaheuristic algorithm), the problem can be stated as follows:
where
is the design vector (each design variable has a value of 0 or
),
is the compliance,
is the displacement vector,
is the load vector, and
is the total volume. The formulation above (Equations (6) and (7)) is combined to define what is called the penalty (merit) function, which can be used in metaheuristic algorithms as follows:
where
is a constraint violation function,
is a penalty function exponent,
is the limit number of iterations, and
is the current iteration. Note that the penalty function exponent is a linearly increasing function (
–
) [
31]. In the first iterations of the BBA, a relatively large volume constraint violation was observed; thus,
is about the value of
. In the later iterations, the
increases to guarantee that feasible designs were obtained (not even a very small volume constraint violation). Solving the numerical design problems (
Section 5) shows that this type of merit function is very effective in comparison with penalizing the violation using a constant positive number (i.e.,
).
3. Binary Bat Algorithm
In nature, bats use sonar to navigate and hunt. They tend to increase the rate and decrease the loudness of emitted ultrasonic sounds when they chase prey. In the BA, virtual bats (possible designs) have the same two main characteristics. Namely, an artificial bat has a position vector, velocity vector, and frequency. The algorithm starts by initially generating a set of random designs from the design domain. Then, in every iteration, new designs are generated and analyzed. Mathematically, an artificial bat’s location vector, velocity vector, and frequency are updated as follows [
32]:
where
is the
-th virtual bat’s location (the
-th design vector),
is the
-th virtual bat’s velocity,
is the best design obtained so far,
is the frequency of the
-th virtual bat,
and
are the minimum and maximum frequencies respectively, and
is a random number of a uniform distribution
. Initially, each bat is randomly assigned a frequency which is drawn uniformly from [
,
].
and
are algorithmic parameters that set depending on the problem. Equations (11)–(13) guarantee the exploitability phase of the algorithm. However, the exploitation phase is performed by randomly selecting a solution among the current best solutions to perform a random wall as follows:
where
is a random number of a uniform distribution
, and
is the average loudness of emitted sound that bats use to perform exploration at the current step. The pseudocode of BA is shown in Algorithm 1.
In Algorithm 1, note that
and
are random numbers of a uniform distribution
, and
is the compliance (Equation (8))
and
(if they are not set to constant) are updated as follows:
where
and
are constants. For any
and
in the end,
and
converge to zero and
, respectively [
21].
and
are algorithmic parameters that set depending on the problem.
In the BBA, the search space is considered a hypercube. That is, artificial bats can shift to the nearer and farther corners of this hypercube by flipping bits (i.e., switching between “0” or “1” values). Therefore, the velocity and position updating expressions (Equations (11) and (12)) must be modified to be consistent with the bit representation. This is done by using a transfer function that is able to map the search process in a continuous search space to a binary search space while preserving the main characteristics of the evolutionary algorithms. Mirjalili et al. [
24] used a v-shaped transfer function and the position updating rule to perform the mapping as follows:
where
and
are the position and velocity of the
-th bat at iteration
in
-th dimension respectively,
is the complement of
, and
is a random number of a uniform distribution
. This step is a major difference between the BA and the BBA. It is used to force artificial bats to move in a binary search space. The pseudocode of the BA is shown in Algorithm 2. The difference between the two algorithms are underlined.
Algorithm 1. Pseudocode of BA [21]. |
- Inputs:
The population size, maximum number of iterations, minimum and maximum pulse frequencies, loudness, and pulse rate. Initialize the bat population positions and velocities. Define pulse frequency.
|
Outputs: The location of the best bat (design) and its fitness value. |
while |
Generate new designs by adjusting frequencies, updated velocities, and locations (Equations (11)–(13)). |
if |
Select a solution among the best solutions randomly. |
Generate a local solution around the selected best design. |
end if |
Generate a new solution by flying randomly. |
if |
Accept the new solution. |
Reduce and increase |
end if |
Rank the bats and find the current best design () |
end while |
Algorithm 2. Pseudocode of BBA [24]. |
Inputs: The population size, maximum number of iterations, minimum and maximum pulse frequencies, loudness, and pulse rate. Initialize the bat population positions and velocities . Define pulse frequency. |
Outputs: The location of the best bat (design) and its fitness value. |
while |
Generate new designs by adjusting frequency and updated velocities. |
Calculate transfer function value using Equation (17) and update positions using Equation (18). |
if |
Select a solution among the best solutions randomly. |
Change some of the dimensions of position vector with some of the dimensions of . |
end if |
Generate a new solution by flying randomly. |
if |
Accept the new solution. |
Reduce and increase |
end if |
Rank the bats and find the current best design () |
end while |
4. Topology Optimization Based on the BBA
The design domain of the two-dimensional continuous structure is discretized into small and square elements in the finite element model (FEM), where each element represents either a material or a void. The binary string or bit-array is used to store the design information. In the string, the code “0” refers to a void element, while the code “1” represents a solid element in FEM. In this study, the stiffness of the void element is considered
instead of zero stiffness to avoid the singularity of the global stiffness matrix in FEM. The decoding step of the binary string into a two-dimensional design domain (topology design) is simple, as shown in
Figure 1. In particular, for a 5 × 5 design problem, a binary string of 25 bits is generated (
Figure 1a). Then, a 5 × 5 matrix is organized as follows: the first five bits construct the first row, the bits from 6 to 10 contract the second row, and so on (
Figure 1b). Finally, the code “0” is replaced by a void element, while the code “1” is replaced by a solid element (
Figure 1c). If a matrix has
rows and
columns, the mathematical expression to map a string into a matrix is as follows:
where
represents the location of the element in the binary string,
is the integer division operation, and
is modulo division operation.
During the optimization process, the BBA generates binary strings (designs) to be evaluated. These designs are mapped to the topology designs as described in the previous paragraph. Then, FE analyses are carried out to calculate each design’s compliance to find the associated penalty function value (Equation (8)). The current iteration penalty function values are used to find the best design of the current iteration and to generate the next iteration designs (as briefly described in
Section 3).
This presentation method does not prevent the formation of the following issues (Wang and Tai [
4]):
Design connectivity must be handled properly in the structural topology design to ensure the success of the BBA. A connected structure is referred to as a group of solid elements, each of which shares at least two edges with its neighboring elements with the exception of support and load elements that can share one edge with neighboring elements. Because of the four issues described above, it was necessary to develop an algorithm to modify the BBA designs before FE analysis. The proposed algorithm is based on the algorithm developed by Wu and Tseng [
3], and new steps were added to suit the BBA and improve its convergence rate. The algorithm is stated as follows:
Elements where supports or loads should be selected as the starting elements since these elements must exist. Code “2” is assigned to the staring elements. Now elements with supports and loads have the code “2” and the rest of the elements have either the code “0” or “1”.
If an element with the code “1” shares an edge with an element with the code “2”, the code “1” is changed into the code “2”. This procedure is repeated until no element with the code “1” can be changed into the code “2”.
Elements with the code “2” represent a continuous structure, and elements with the code “1” represent a discontinuous structure(s). The code “1” is changed into the code “0” while the code “2” is changed into the code “1”. In this way, the binary string representation of the structure topology is restored.
If an element with the code “0” shares no edge with another element with the code “0” (i.e., void element surrounded by solid element), its code is changed into the code “1”.
If an element with the code “1” is adjacent to only one element with the code “1” (i.e., solid element adjacent to three void elements), its code is changed into the code “0”.
Finally, after steps 1 to 5 are applied, the modified design replaces the original design before step 1. It was observed that this step improved the convergence behavior of the algorithm because in the next iteration, the BBA tried to improve designs without isolated objects, many scattered void elements, and solid elements that shared one edge with the continuous structure. Moreover, few iterations of modification were needed (steps 1–5) in the next iteration designs because the modification of any new design occurred where there were new isolated objects, scattered void elements, and solid elements that shared one edge with the continuous structure. This step sped up the algorithm because few iterations of modification were needed.
The initial design population (a matrix of all designs) must have one vector of ones. This design is a solid structure with no voids. The rest of the designs are generated randomly.
Steps 1 to 3 solved the discontinuity problem (issue 1), step 4 solved the checkerboard and scattered void elements problems (issues 2 and 3), and step 5 solved issue 4.
Figure 2 shows that steps 4 and 5 were necessary to improve the quality of the topology designs. That is, for the first design example, at iteration 1000, the compliance was reduced by about
(there was no volume violation in both cases) by adding steps 4 and 5 to the algorithm. Moreover, the design that is shown in
Figure 2b is more practical in comparison with the design shown in
Figure 2a.
5. Numerical Examples
Four benchmark structural topology optimization problems were employed to examine the search capability of the bit-array representation with the BBA and the proposed filtering algorithm and merit function (Equation (8)).
The structures were analyzed using the finite element (direct stiffness) method, and the algorithms were coded using MATLAB. The BBA algorithmic parameters were set as follows: population size is , the maximum number of iterations () , loudness () is , and pulse rate () is . These parameters are selected depending on the problem. Thus, in this study, they are set based on the experience gained after some testing runs. The minimum and maximum pulse frequencies are 0 and 2, respectively.
Because the BBA is stochastic in nature, five independent optimization runs were performed for each example and the best topology is reported here. It is interesting to note that each independent run produced not only a different compliance but also different configuration because of the random search nature of metaheuristic algorithms.
In all problems, the objective was to find the minimum compliance subjected to a maximum volume constraint of
(i.e.,
), Young’s modulus is
unit force/unit area, and Poisson’s ratio is
. For design examples 1, 2, 4, and 5, the design domain is rectangular of
unit lengths in the
direction and
unit length in the
direction. For the finite element analysis, the domain is divided into
quadrilateral square plane stress elements. That is, the total number of design variables is
because each element is a design variable with 0 or 1 values. For design example 3, the design domain is rectangular of
unit lengths in the
direction and
unit length in the
direction. For the finite element analysis, the domain is divided into
quadrilateral square plane stress elements (7500 elements). This mesh size is chosen because the SIMP method becomes less mesh-dependent when the mesh is fine. Although a finer mesh produces a slightly better design, it is more computationally expensive. The numerical results from the present bit-array representation with the BBA are compared with those from the SIMP method using Sigmund’s 99 line topology optimization code written in MATLAB [
27]. In the power law approach, the penalization power and the filter size are set to 3 and 1.2, respectively.
5.1. The Cantilever Plate Problem (Case 1)
The problem is illustrated in
Figure 3, which shows a cantilever plate subjected to a point load. That is, the left boundary is fixed, and a unit point load is applied vertically downward at half the height of the right boundary. The problem was previously solved by Wu and Tseng [
3] and Wang and Tai [
4] using GA and a mesh size of
(
elements).
Table 1 shows the topology designs obtained by the SIMP method, the BBA without symmetry (5000 design variables), and the BBA with symmetry (2500 design variables), as well as their corresponding compliances. The symmetry is taken around an axis parallel to the
-axis and passing through
on the
-axis. Although the SIMP method with a mesh size of
produced distinct boundaries, many gray elements were found due to the smoothing filtering across the edges. The mass density is taken as a continuous design variable rather than a discrete binary variable [
3,
4]. Using the SIMP method, the minimum compliance of
is obtained at iteration 184. With a mesh size of
(9800 elements), the minimum compliance of
is obtained at iteration 557. However, the computation time increased about 33 times. This result shows that a mesh size of
is resealable for the types of problems discussed in this study. Topology optimization using the present bit-array with the BBA is performed without and with symmetry. In both cases, it can be seen that the distinct black (material) and white (void) design topology are obtained because discrete binary design variables are employed. The minimum compliances of
and
are obtained from the BAA without symmetry and the BBA with symmetry, respectively. Minimum compliances of 64.44 and 65.26 were found by Wu and Tseng [
3] and Wang and Tai [
4] using GA, respectively. Therefore, better compliances were found using the present bit-array with BBA but with different configurations (see
Table 1).
5.2. The Cantilever Plate Problem (Case 2)
The problem is illustrated in
Figure 4. This problem is similar to the problem described in
Section 5.1. However, the location of the point load is different. That is, the left boundary is fixed, and the unit point load is applied vertically downward at the bottom right corner. To the best of the authors’ knowledge, this problem has not been previously solved before using metaheuristic algorithms.
Table 2 shows the topology designs obtained by using the power law approach and the BBA without symmetry as well as their corresponding compliances. In this problem, symmetry is not applicable. Minimum compliances of
and
were obtained using the SIMP method and the bit-array with the BBA, respectively. The present bit-array with the BBA resulted in slightly less compliance than the SIMP method did but there were no blurry boundaries.
5.3. The Cantliver Plate Problem (Case 3)
The third example is the
cantilever beam.
Figure 5 demonstrates the problem definition. This problem is similar to the first design example (
Section 5.1). However, in this case, the beam is one unit-length longer. The problem was previously solved by Gao [
10] using the CDT method with a mesh size of
(
elements).
Table 3 shows the final design layout and compliance of the BBA with symmetry. The proposed bit-array with BBA converged to compliance of 177.73 where the CDT and the SIMP methods reached minimum compliances of 182.475 and 182.816, respectively [
10]. The presented result in this study shows no checkerboard patterns or grayscale areas as in those from the CDT and the SIMP methods.
5.4. The Pined-Pined Plate Problem
The problem definition is illustrated in
Figure 6. The nodes at the left and right bottom are set as fixed nodal displacements, and the unit point load is applied vertically downward at half the length of the bottom boundary. This problem was previously solved by Kaveh et al. [
33] using the Ant Colony algorithm with a mesh size of
(
elements). However, the compliance of the final design was not provided.
The results derived by the SIMP method, the BBA without symmetry, and the BBA with symmetry as well as their corresponding minimum compliances are shown in
Table 4. The symmetry is taken around the axis parallel to the
-axis passing through
on the
-axis.
Table 4 demonstrates that the proposed bit-array with the BBA converges to better solutions.
5.5. The Pined-Roller Plate Problem
The problem definition is illustrated in
Figure 7. This problem is similar to the problem in
Section 5.3. However, the pinned support on the right side is replaced by a roller support. A unit point load is applied vertically downward at the half-length of the bottom boundary. This problem was previously solved by Kaveh et al. [
33] using the Ant Colony algorithm with a mesh size of (
elements). However, the compliance of the final design was not provided.
Table 5 shows the layouts obtained after the optimization process and the minimum compliances of the SIMP method, the BBA without symmetry, and the BBA with symmetry. The symmetry is taken around the axis parallel to the
-axis passing through
on the
-axis. Again, the minimum compliances of the proposed bit-array representation are slightly better the SIMP method but with no blurry boundaries and the final design configurations are different.
According to the results discussed in
Section 5.1,
Section 5.4, and
Section 5.5, the BBA with symmetry is the best among the three methods because the number of design variables in the symmetrical problems is reduced to half. Therefore, the number of possible combinations is also reduced to half.