Implicit constraints handling for efficient search of feasible solutions

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 363 (2020) 112917 www.elsevier.com/locate/cma

Implicit constraints handling for efficient search of feasible solutions Amir H. Gandomia ,∗, Kalyanmoy Debb a

University of Technology Sydney, Ultimo, NSW 2007, Australia b Michigan State University, East Lansing, MI 48824, USA

Received 13 June 2019; received in revised form 5 February 2020; accepted 6 February 2020 Available online xxxx

Abstract Real-world optimization problems usually involve constraints and sometimes even finding a single feasible solution is a challenging task. This study introduces a new approach for implicitly handling constraints. The proposed approach reduces the consideration of infeasible solutions by directly updating variable bounds with constraints, which is called the boundary update (BU) method. Two illustrative examples are used to explain the proposed approach, followed by applying it to mathematical and engineering constrained optimization problems. Finally, a surrogate-based problem and a large-dimensional and highly constrained problem are used to evaluate the BU method on these types of problems. The BU method is coupled with seven well-known evolutionary and mathematical optimization algorithms and the results show that the proposed BU method is a practical and effective approach and leads to better solutions with fewer function evaluations in nearly all cases, particularly for population-based optimization algorithms. This study should motivate optimization researchers and practitioners to pay more attention to the direct handling of constraints, rather than constraint handling by penalty or other fix-ups. c 2020 Elsevier B.V. All rights reserved. ⃝ Keywords: Constrained optimization; Engineering optimization; Boundary update; Nonlinear constraints; Evolutionary computation; Global optimization

1. Introductions Engineering problems usually involve constraints, which may be functional, operational, geometrical, or physical. The presence of constraints makes it challenging to handle an optimization problem. In certain highly constrained problems, finding a single feasible solution is often the primary task. Thus, it is not surprising that the literature contains many studies on handling different types of constraint. Because of the complexity of constrained problems in making the search space spotted with fragmented feasible regions, a successful optimization algorithm must have a global perspective and avoid both infeasible regions and locally optimal feasible regions in order to reach the globally optimal feasible region and eventually converge with the globally optimal solution. Moreover, standard penalty function approaches for handling constraints are too rigid and parameter-dependent for them to produce an easy path to a feasible and globally optimal solution. Algorithms that use gradient information cannot differentiate between local and global optima and they are likely to trap into local optima. Also, gradient-based algorithms have issues with the discontinuity of decision ∗ Corresponding author.

E-mail address: [email protected] (A.H. Gandomi). https://doi.org/10.1016/j.cma.2020.112917 c 2020 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

variables and objective landscapes. Metaheuristic algorithms, such as evolutionary algorithms (e.g., [1]) and swarm intelligence (e.g., [2]) methods, use different types of information (usually, non-gradient types) to update solutions from one iteration to next. These algorithms are usually inspired by nature and use evolution and swarm behavior for optimization purpose. Metaheuristic optimization algorithms are global optimization algorithms, meaning that they are less likely to be trapped into the local minima. In addition, these algorithms are flexible enough to be changed with problem information; thus, they can adapt in order to negotiate different types of complexities, nonlinearities, and discontinuities. In this study, a new approach is introduced for implicitly handling constrained problems. The proposed approach narrows down variable space by directly using the constraints and is therefore called the boundary update (BU) method. This approach, by limiting the feasible search space for one or more variables, can be used to limit the search space, thus forcing the algorithm to concentrate its search in the feasible regions. The proposed approach can be used for linear or nonlinear constraints and even problems having constraints defined in a gray/black-box manner, if a defined condition is satisfied. The BU method is used with seven well-known optimization algorithms, including both mathematical and metaheuristic optimization algorithms. Several constrained and global optimization problems, including mathematical and engineering case studies, have been investigated in this study. When the BU method is not directly applicable, surrogate models can be used, which is discussed in a case study. Also, a highly constrained large-dimensional optimization problem is solved by using the proposed approach in order to show performance with a case in which even finding a feasible solution is challenging. The proposed BU approach was effective in nearly all cases. Not only does the BU approach allow the algorithms to quickly reach feasible regions, it also improves the convergence rate in reaching a better solution. In the remainder of this paper, Section 2 describes the BU method along with two illustrative examples and also explains the optimization methods used in this study. The numerical studies and their results are presented in Section 3. The detailed formulations of these case studies are presented in the supplementary material. Summary and conclusions are described in the last section of this paper. 2. Boundary update approach 2.1. Methodology A constrained optimization problem can be formulated as follows: min f (x) x∈Ω

subject to g j (x) ≤ 0 for j ∈ {1, . . . , m} where x is a vector of decision variables, x = (x1 , . . . , xn ), which is subjected to lower bound and upper bound vectors as Ω = [l b, ub], and f (x) is the objective function. Each decision variable boundary could be either finite or infinite. Generally, optimization is an iterative process, which continues until at least one of the stopping criteria is met (shown in Fig. 1(a)). As can be seen in Fig. 1, handling the boundaries and constraints usually happens after applying the search operators. Here, the proposed approach intends to handle constraints when it checks the variable boundaries, which results in more efficient searching for feasible solutions. The main philosophy of( the)proposed BU approach is that if a constraint function, i.e., g j (x) ≤ 0, can be rewritten as either xi ≥ li, j x ̸=i or xi ≤ u i, j (x ̸=i ), it can be considered as a function ( ) for the ith variable boundary in which x ̸=i represents the variable vector without the ith variable and li, j x ̸=i and u i, j (x ̸=i ) are the dynamic lower and upper boundaries for the ith variable, respectively, which are functions of variables, except x i . Therefore, the boundaries are now changing and are updated during an(optimization process. ) ( In)general, both lower and upper bounds can be updated so that xi can be bounded as li, j x ̸=i ≤ xi ≤ u i, j x ̸=i . As previously mentioned, a constraint should be solved by the ith variable in order to be considered as the lower or upper bounds of a variable. This condition can be mathematically presented for the ith variable as follows: [ ( ) ] ∃ i ∈ {1, . . . , n} : ∀ j ∈ {1, . . . , m} : xi ≥ li, j x ̸=i ∨ xi ≤ u i, j (x ̸=i ) (1) Definition 1. A variable with an updated boundary is called a repairing variable and a variable without any updated boundary is called a non-repairing variable.

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

3

Fig. 1. (a) Iterative search procedure, (b) Iterative search procedure for handling constraints and boundaries simultaneously.

Definition 2. l and u are lower and upper boundary functions, respectively. Definition 3. Considering the dynamic nature of the ith decision variable boundaries in each iteration, the original variable bounds of the ith decision variable, lbi and ubi , are now called updated lower and upper bounds (lbiu and ubiu ), respectively. If there are no original bounds for a repairing variable (in other words lbi = −∞ and ubi = +∞), the boundary functions are the only boundaries for the variable and, therefore, the updated boundaries can be defined as: ( ) lbiu = li, j x ̸=i (2a) ( ) u ubi = u i, j x ̸=i (2b) However, if the original bounds (lbi and ubi ) already exist for a repairing variable, then the updated bounds can be defined as follows to avoid violating the original bounds: ( ( ( ) ) ) lbiu = min max li, j x ̸=i , lbi , ubi (3a) ( ( ( ) ) ) u ubi = max min u i, j x ̸=i , ubi , lbi (3b) It should be noted that if both boundaries of a repairing variable are updated, the updated boundaries may conflict. In other words, ubiu could become smaller than lbiu or vice-versa. For example, if the upper bound is updated after the lower bound, the upper bound should not be less than the updated lower bound. In this case, the updated upper bound should be changed to: ( ( ( ) ) ) ubiu = max min u i, j x ̸=i , ubi , lbiu (4) The same process should be considered if the lower bound is updated after the upper bound. For the remainder of this paper, in order to simplify the formulation, it is assumed the upper bound is always updated after the lower bound. Using the updating strategy, xi is placed within the updated boundaries in order to satisfy the jth constraint. If a repairing variable, xi , is used to handle more than one constraint (e.g. the first ki constraints), the updated boundaries can be written as follows: ( [ ( ) ( ) ] ) lbiu = min max li,1 x ̸=i , . . . , li,ki x ̸=i , lbi , ubi (5a) ( [ ( ) ( ) ] ) u u ubi = max min u i,1 x ̸=i , . . . , u i,ki x ̸=i , ubi , lbi (5b) where ki ≤ m. Fig. 2 shows a schematic example of multiple functions for updating the boundaries of xi . Two lower bound functions (li,1 and li,2 ) and two upper bound functions (u i,1 and u i,2 ) were used for updating the boundaries. Based

4

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Fig. 2. A schematic of multiple functions to update decision variable boundaries.

on Eqs. (5a) and (5b), lbiu = li,1 and ubiu = u i,2 . Based on Fig. 2, the following conclusions can be made for the seven different ranges of xi : (1) If xi < li,1 , li,2 , lbi → li,1 , li,2 , lbi are violated by xi (2) If li,2 ≤ xi < li,1 , lbi → li,1 , lbi are violated by xi (3) If lbi , li,2 ≤ xi < li,1 → li,1 is violated by xi (4) If li,1 ≤ xi ≤ u i,2 → xi is a feasible solution (5) If u i,2 < xi ≤ u i,1 , ubi → u i,2 is violated by xi (6) If u i,2 , ubi < xi ≤ u i,1 → u i,2 , ubi are violated by xi (7) If u i,1 , u i,2 , ubi < xi → u i,1 , u i,2 , ubi are violated by xi It should be noted that this process still computes each constraint only once and, therefore, there is no computational cost added to the problem from this perspective. However, solving the inequality constraints with respect to one variable (a repairing variable) in the problem formulation step has an additional computational cost. More than one repairing variable may be defined for handling constraints of an optimization problem. If another repairing variable (e.g. xr ) is introduced, it cannot be involved in updating the boundaries of the first repairing variable (xi ), because xr can change after defining its updated boundaries in order to stay within the bounds. In this case, xi boundary functions should exclude xr as well, which can be formulated as follows: ( [ ( ) ( ) ] ) lbiu = min max li,1 x ̸=i,r , . . . , li,ki x ̸=i,r , lbi , ubi (6a) ( [ ( ) ( ) ] ) u u ubi = max min u i,1 x ̸=i,r , . . . , u i,ki x ̸=i,r , ubi , lbi (6b) where ki ≤ m. In this case, the first repairing variable (xi ) can be used for updating the boundaries of the second repairing variable (xr ), as its boundaries have already been checked. For example, if the next kr constraints (ki + 1 to kr constraints) are used to update the second repairing variable bound(s), the updated boundaries of xr can be written as follows: ( [ ( ) ( ) ] ) lbru = min max lr,ki +1 x ̸= r , . . . , lr,ki +kr x ̸= r , lbr , ubr (7a) ( [ ( ) ( ) ] ) u u ubr = max min u r,ki +1 x ̸= r , . . . , u r,ki +kr x ̸= r , ubr , lbr (7b) where r ∈ {1, . . . , n} and r ̸= i. Lemma 4. A repairing variable cannot be used for updating the boundaries of another repairing variable, unless it is updated beforehand. In other words, repairing variable boundaries must be checked before they can be used for updating the function of another repairing variable. Using this method, decision variables space (Ω ) can be split into repairing variable(s) space (Ωb ) and non-repairing variable(s) space (Ωn ) and, therefore, Ω is divided into Ωb and Ωn , where Ωb ∪ Ωn ≤ Ω and Ωb ∩ Ωn = ∅. Applicability condition: Based on the above discussion, in order to apply the proposed method to satisfy the jth constraint using the ith repairing variable, the following condition should be satisfied, which can be considered as a generalized form of Eq. (1): [ ( ( )) ( )] − − ∃ xi ∈ Ωb : ∀ j ∈ {1, . . . , m} : xi ≥ li, j x ∈ Ωn ∪ Ωb,i ∨ xi ≤ u i, j (x ∈ Ωn ∪ Ωb,i ) (8) − + where Ωb,i and Ωb,i denote the repairing variable (xi ) spaces before and after updating, respectively. Following this procedure, any number of variables’ boundaries may be updated to implicitly handle the constraints of a problem. By keeping a repairing variable within the updated boundaries, the involved constraints are satisfied since the new bounds intend to satisfy these constraints.

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

5

Lemma 5. After updating a boundary, the updated boundary is only valid until the end of that iteration, which means the updated boundary values cannot be used in other iterations. Constraint Satisfactory Condition: As previously mentioned, if the value of an updated upper bound is less than an updated lower bound (ubiu < lbiu ) for a repairing variable, this means that the updated bounds conflict. In other words, the variable does not have any feasible range. In this case, any value of xi violates at least one constraint. If the original bounds do not exist for the repairing variable, the considered constraints are satisfied, as long as the variable is placed within the non-conflicted updated bounds. If the original boundaries are already specified for the repairing variable, whether the bound functions u or l violate the original boundaries (ub and lb) should be checked first. Then, whether the updated boundaries conflict must be checked. Therefore, there is a range over which a variable repairing variable can satisfy the involved constraints if the following condition is satisfied: ) ) ( ( (9) max lbi , lbiu ≤ min ubi , ubiu It should be noted that in a population-based algorithm the same repairing variable could have different boundaries in different solutions since its boundaries are functions of the other solution variables. In other words, the updated boundaries are varied from solution to solution. In general, the variable boundaries should be checked for every new solution. Using the updated procedure described above, the feasible range for a repairing variable is likely reduced from its original range to satisfy the constraints. Updating the range happens in each iteration, which causes discontinuity among iterations. To overcome this issue and to only search within the updated boundaries, a mapping scheme is proposed in this work by replacing a repairing variable (xi ), which has updated boundaries, with a new variable ( pi ). Then, a new variable ( pi ) can be mapped to its actual variable, using the following formula: ( ) xi = lbiu + pi ubiu − lbiu (10) where 0 ≤ pi ≤ 1 and therefore pi indicates the relative distance from the updated lower bound to the updated upper bound. Definition 6. A p variable can be considered as a general form of the semi-independent variable (SIV), which means that a variable can be considered as a fraction of other variable(s) [3]. Therefore a p variable is called a generalized SIV (GSIV). Using this approach, a repairing variable (xi ) is substituted with the GSIV ( p) when search operators are applied and then the p remaps to the actual repairing variable for the rest of the search cycle, as shown in Fig. 1. Considering the first h variables as repairing variables, the new optimization framework is expressed as shown in Algorithm 1. Fig. 1(b) illustrates how the two stages of handling boundaries and constraints can merge by using the BU method.

This approach changes the formulation of the problem to narrowing down the search space to the feasible region. Moreover, based on the framework described for this approach, the formulation of the problem can be changed as

6

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

follows: min

x∈(Ωb ∪Ωn )

f (x)

subject to g j (x) ≤ 0, j ∈ {1, . . . , m} ] [ where Ωb = l bu , ubu and Ωn represents the variable space of non-repairing variables. The proposed approach still needs to be coupled to an explicit constraint handling scheme, such as penalty function, since applicability and satisfactory conditions may not be satisfied in some cases. Therefore, this approach is called an implicit constraint handling approach. The BU method narrows down the search process to the remaining feasible region with respect to the considered constraints. The proposed method can be used to handle any number of constraints. This process helps to reduce the computational burden, as it eliminates the creation of solutions in the infeasible region. The proposed approach eliminates the infeasible search space and transforms a constrained optimization problem into an unconstrained optimization problem during iterations. Unlike cutting plane methods [4], the proposed BU approach does not approximate the constraint structures, and can simply be used to handle a wide range of constraints. By using the proposed method, we can implicitly handle a constraint as long as one of the repairing variables can be either explicitly or implicitly rewritten in terms of other variables. In other words, our approach is efficient as long as a constraint, g(x) ≤ 0, can be written as li (x ̸=i ) ≤ xi and/or xi ≤ u i (x ̸=i ). Here, li and u i can have any form and can even represent a black-box function, which can occur in many real-world problems. Therefore, our proposed approach has a wider range of application potential than the existing approximated cutting plane methods. A few methods in the literature attempted to map the entire search space to the feasible search space, such as decoder-based methods. In these methods, every solution under consideration always represents a feasible solution [5]. The original variables are mapped in a space where every point represents a feasible solution in the original space, and a complicated but implicit mapping is obtained by pivoting solutions around the current best solution. The idea is theoretically promising, but the actual implementation of decoder-based mapping involves a high computational cost, which makes the overall approach challenging to use in practice [6]. We believe that our proposed approach represents a compromise between decoder-based concepts and the ease of computational efforts. 2.2. Variable selection Multiple candidates may be considered as repairing variables in a constrained optimization problem and therefore selecting the best repairing variable(s) should be studied. Here, a simple selection procedure is proposed, as shown in Algorithm 2. The purpose of this repairing variable(s) selection strategy is to handle the maximum number of constraints with the minimum number of repairing variables.

For example, if we have the following repairing variable candidates: • • • •

x1 x2 x3 x4

handles handles handles handles

g1 , g2 , g3 g2 , g4 g3 , g5 g6 , g7

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

7

Based on the proposed strategy, the first repairing variable is x4 , which handles two constraints that are not handled by other decision variables. Each of the other three variables handles only one non-conflicting variable. Based on the proposed strategy, among these variables, x1 is selected as the second repairing variable, as it handles three constraints (g1 , g2 , and g3 ) in total, which is more than the other two decision variables (x2 and x3 ). After selecting x1 as the second repairing variable, the only remaining constraints are g4 and g5 and, therefore, x2 and x3 are randomly selected as the third and fourth repairing variables. 2.3. Illustrative examples 2.3.1. Demonstration of the BU approach Rosenbrock’s Banana or valley function is one of the most well-known optimization problems that exemplifies the importance of variable-linkage in optimization (Fig. 3). The problem has a parabolically shaped valley that contains the global minimum, thereby making the problem difficult to solve. Adding constraints to this problem makes it even harder, especially when the global minimum is not the same point as the global minimum of the unconstrained problem. The 2D version of this problem with constraints can be formulated as follows: ( )2 min: f (x1 , x2 ) = 100 x2 − x12 + (1 − x1 )2 (11) s.t.:

g1 (x1 , x2 ) = 1.5 + x1 x2 − x1 − x2 ≤ 0

(12a)

g2 (x1 , x2 ) = −x1 x2 − 10 ≤ 0

(12b)

where −5 ≤ x1 ≤ 5 and −10 ≤ x2 ≤ 20, and the global optimum located at [−1.108, 1.237]. Here both decision variables can be used as repairing variables to handle both constraints; x1 is selected for this purpose. After solving the inequality constraints, g1 and g2 , with respect to x1 (repairing variable), the repairing variable can be used as the lower and upper bound functions. These bound functions for the first and second constraints can be presented as follows: ⎧ ⎨ x2 − 1.5 x2 < 1 (13a) l1,1 = x2 − 1 ⎩ −∞ other wise ⎧ ⎨ +∞ x2 ≤ 1 u 1,1 = x2 − 1.5 (13b) other wise ⎩ x2 − 1 ⎧ ⎨ −∞ x2 ≤ 0 l1,2 = −10 (13c) other wise ⎩ x ⎧ 2 ⎨ −10 x2 < 0 u 1,2 = (13d) x2 ⎩ +∞ other wise The ±∞ indicates there is no bound or limitation. Using the bound functions, the updated bounds are as follows: ( [ ] ) lb1u = min max l1,1 , l1,2 , −5 , 5 (14a) ( [ ] ) u ub1 = max min u 1,1 , u 1,2 , 5 , −5 (14b) Finally, x1 can be calculated as: ( ) x1 = lb1u + p ub1u − lb1u

(15)

where 0 ≤ p ≤ 1. The population simulation of 1000 solutions for this problem using the BU method is shown in Fig. 3. Instead of distributing the solutions throughout the entire search space, the solutions are distributed within the feasible search space. In other words, variable x1 , which was originally defined to lie within −5 to 5, is now restricted to be between two functions, lb1u and ub1u , as illustrated in Fig. 3. It should be noted that points created in the (p, x2 )

8

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Fig. 3. Contour plot of the search space with highlighted constraints and simulations of the population (p, x2 ).

space are only placed in the feasible region. This strategy can narrow down the infeasible area with respect to the updated boundaries; therefore, not only can it help an optimization algorithm avoid going into the infeasible region, but can also allow it to search the feasible area more efficiently. In this illustrative example, the BU method can eliminate the entire infeasible search space. However, for complex problems, repairing the feasible search space may occur during search iterations by involving all constraints for boundary updating. 2.3.2. Conditions assessment The BU approach works perfectly for the previous illustrative example. However, as previously mentioned, a constraint should satisfy the applicability condition in order to be handled by the proposed approach. The current example can have difficulties with satisfying this condition. The following problem was introduced by Deb [7] as a challenging constrained optimization problem: )2 ( )2 ( (16) min: f (x1 , x2 ) = x12 + x2 − 11 + x22 + x1 − 7 s.t.:

g1 (x1 , x2 ) = 4.84 − (x1 − 0.05)2 − (x2 − 2.5)2 ≥ 0

(17a)

g2 (x1 , x2 ) =

(17b)

x12

+ (x2 − 2.5) − 4.84 ≥ 0 2

where 0 ≤ x1 ≤ 6 and 0 ≤ x2 ≤ 6, and the global optimum located at [2.2468, 2.3819]. The contour plot of the objective function along with the constraint lines is shown in Fig. 4. As shown in Fig. 4, the feasible area has a narrow crescent-shape. The total feasible search region of this problem is approximately 0.7% [7]. If x2 is considered as the repairing variable, it will have two separate ranges for x1 < 2.5. For example, for x1 = 1.5 (dashed line), there are two separate feasible areas for x2 (around x2 = 1 and 4). Therefore, the applicability condition cannot be satisfied unless a more sophisticated strategy is used with a book-keeping or a linking of these feasible ranges. Here, x1 can be used as a repairing variable since a function of x2 can limit x1 to

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

9

Fig. 4. Upper and lower functions of x1 based on x2 .

lie within the crescent feasible region, as follows: √ √ 2 0.05 − 4.84 − (x2 − 2.5) ≤ x1 ≤ 0.05 + 4.84 − (x2 − 2.5)2 √ √ 4.84 − (x2 − 2.5)2 ≤ x1 and x1 ≤ − 4.84 − (x2 − 2.5)2

(18a) (18b)

where Eqs. (18a) and (18b) represent the solutions of g1 (x1 , x2 ) and g2 (x1 , x2 ), respectively, which are illustrated √ in Fig. 4, which shows that both constraints have two solutions. For Eqs. (18b), it is trivial that x1 ≤ − 4.84 − (x2 − 2.5)2 is not within the search space, since 0 ≤ x1 and, therefore, can be eliminated. This shows the importance of considering the variable bounds when an equation is solved, particularly the constraints with more than one answer. Therefore, the new bounds can be defined as follows: ( ) ( ) √ √ √ max 0.05 − 4.84 − (x2 − 2.5)2 , 4.84 − (x2 − 2.5)2 , 0 ≤ x1 ≤ min 0.05 + 4.84 − (x2 − 2.5)2 , 6 (19) A population simulation of 1000 solutions for this optimization problem using x1 as the repairing variable is illustrated in Fig. 5. Initially, the simulated points are mostly located within the feasible search space, which is an interesting result for a problem with such a small feasible area. As can be seen in Fig. 5, there is no feasible area for the low and high x2 and therefore all solutions are placed in the infeasible region. During boundary updating of these values, the updating function is less than the original lower bound and therefore all x1 values in these ranges reflect the original lower bound. In other words, the updated function violates the lower bound of x1 and therefore all the x1 values in these two ranges stick to the lower bound. The performance of the proposed BU approach was investigated in Section 3 by using several complex constrained optimization problems. 2.4. Optimization algorithms Optimization algorithms can be classified based on different features. For example, they can be classified as trajectory-based (single solution-based) or population-based optimization algorithms, or they can be classified as mathematical and nature-inspired (or metaheuristics) algorithms. In general, mathematical optimization algorithms

10

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Fig. 5. (a) Global optimum of the constrained problem, (b) Contour plot of the search space with highlighted constraints and simulation of population.

are trajectory-based and use derivative information in order to improve the solution. Therefore, these gradient-based algorithms have a greater potential of being trapped in a local minimum as compared to metaheuristic algorithms, which use other sources of information such as evolution, swarm behavior, and annealing. Furthermore, these gradient-based optimization algorithms are more sensitive to the nonlinearity and discontinuity of the variables and landscape. Metaheuristics have been widely used during the last two decades and have remained a popular research topic. These techniques have been applied to many complex real-world problems [8], particularly to constrained problems (e.g., [9]). Metaheuristics are a subset of artificial intelligence; however, they are slightly different from the classical methods in that their intelligence comes from nature. The efficiency of metaheuristics is due to their significant ability to imitate the best features of nature, which have evolved by natural selection over millions of years. Most metaheuristics are population-based algorithms and gradient-free algorithms. In the BU method, the search operator(s) do not directly apply to the repairing variable(s). Instead, they apply to the GSIV(s), which subsequently are mapped to repairing variable(s) to calculate the objective and constraint values. This mapping makes the BU method challenging to implement for the optimization algorithms that use gradient information and therefore gradient-based optimization algorithms may not be suitable for the proposed framework. However, metaheuristics can be well-adopted tools for such a framework, as they can adapt with the mapped variables and discontinuities in the search space. To investigate the performance of different optimization algorithms, seven different optimization algorithms, including both mathematical and metaheuristic algorithms, were used, including: • • • • • • •

Interior-Point (IP) Sequential Quadratic Programming (SQP) Direct Search (DS) Simulated Annealing (SA) Particle Swarm Optimization (PSO) Genetic Algorithm (GA) Differential Evolution (DE)

The first four algorithms (IP, SQP, DS, and SA) are trajectory-based optimization algorithms and the PSO [10], GA [11], and DE [12] algorithms are population-based optimization algorithms.

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

11

2.4.1. Implementation of optimization algorithms To evaluate the proposed BU method for the first six optimization algorithms (IP, SQP, DS, SA, PSO, and GA), MATLAB® optimization and global optimization toolboxes were used, which include well-tuned and wellstabilized algorithms. The default and suggested values of MATLAB® were used as parameters for the algorithms. A developed standard code was used for the DE algorithm, which is based on the basic algorithm. The DE parameters differential weight (F) and crossover probability (Cr ) were r assigned as 0.7 and 0.9, respectively, based on suggestion from the literature [13]. For all seven algorithms, the same maximum number of function evaluations were set; for the population-based algorithms, the same population sizes were used. It should be noted that all the algorithms were stopped if the relative change in the best objective function value over the last 20 iterations was less than 1 × 10−6 . Each of these seven algorithms runs 51 times with random initial point(s), because of the stochastic nature of the algorithms. As previously mentioned, our approach represents an implicit constraint handling approach and therefore a penalty method is coupled with the BU method as an explicit constraint handling scheme to penalize a solution that still violates constraint(s). Here, a static penalty method was implemented so that objective values of infeasible solutions were penalized based on their violations, defined as follows: f itness (x) = f (x) + µ ×

m ∑

violation i

(20)

i=1

where violation i is the absolute value of the jth constraints violation and 1 ≤ µ is the penalty coefficient that should be sufficiently large, depending on the quality of the solution needed. It should be noted that the static penalty will only be applied to constraints not satisfied by the BU method. 3. Numerical study The defined BU method was applied to several constrained problems that satisfy a defined condition. First, two constrained mathematical problems were solved. The first mathematical problem has linear constraints and the performance of BU was evaluated for three different strategies. Unlike the first problem, the second, called Himmelblau’s problem, has non-linear constraints. Next, eight well-known constrained engineering optimization problems, borrowed from the literature, were considered, including pressure vessel design, welded beam design, a tension–compression spring design, a speed reducer design, a corrugated bulkheads design, a rolling element bearing design, and a multiple disk clutch brake design. Following these case studies, the proposed approach was used to optimize a surrogate model for a car side impact design problem. This case study can be considered as an example of using the proposed approach for a black-box problem. Finally, a high-dimensional and highly constrained problem (a 100D problem with 100 constraints) was solved to further assess our proposed method. All problem formulations are presented in the supplementary materials. The case studies along with their sections and references are listed below: • Section 3.1: Constrained Mathematical Problems: ◦ G1 Problem [14] ◦ Himmelblau’s problem [15] • Section 3.2: Constrained Engineering Problems: ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦

Welded Beam [16] Pressure Vessel [17] Tension–Compression Spring [18] Speed Reducer [19] Corrugated bulkheads design [20] Heater exchanger [21] Multiple disk clutch brake [22] Rolling element bearing [23]

12

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

• Section 3.3: Surrogate model optimization: ◦ Car side design [24] • Section 3.4: Highly Constraint and Large-dimensional Problem: ◦ Stepped beam design problem [25] 3.1. Constrained mathematical problems Two mathematical problems were considered. The constraints of the first problem are linear and the second problem has nonlinear constraints. These two problems were used to investigate the effect of linear and nonlinear constraints. The maximum number of function evaluations (FEs) for these two problems were selected to be 5000 and 20,000, respectively, which enabled the population-based algorithms to perform 100 and 250 maximum iterations, respectively. 3.1.1. Linear constraints: G1 problem This problem is a well-known mathematical constrained problem with 13 variables. The objective function is quadratic, but all nine constraints of this problem are linear. The ratio of the feasible search space over the entire search space (hit ratio) is approximately 0.0111% for this problem. Therefore, this is a good case for studying a highly constrained problem with linear constraints. In this problem, six of nine constraints are active in the global optimum. The BU method is applied in three different steps in which it intends to handle all nine constraints in three steps, with three constraints in each step. One, two, and three repairing variables are used in steps 1 to 3, respectively. These steps are defined as follows: Step 1 (1p): x10 as a repairing variable for handling three constraints (g1 , g4, and g7 ) Step 2 (2p): x10 and x11 as repairing variables for handling six constraints (g1 , g2 , g4 , g5 , g7, and g8 ) and Step 3 (3p): x10 , x11 and x12 as repairing variables for handling all nine constraints; The problem is optimized by using IP, SQP, DS, SA, PSO, GA, and DE algorithms without BU and the three defined steps with 1p, 2p, and 3p. The statistical results of the algorithms without the proposed approach and with the three steps are presented in Table 1, which shows that the IP, SQP, and DS methods found feasible solutions only if all constraints were handled by the proposed approach (3p). The Wilcoxon rank-sum test at a significance level of 5% was performed to compare the results obtained step-by-step. For all metaheuristics except the SA algorithm, it can be seen that with handling more variables, the statistical results are significantly improved, which is consistent for 1p to 3p. The SA algorithm is not suitable for this case study, as it does not stop until the maximum number of FE is reached. The population-based algorithms were successful in finding a feasible area in nearly all cases. The mean convergence histories of these algorithms are presented in Fig. 6(a) to (c). These convergence histories were improved at almost any point during the search by handling more constraints. Based on Fig. 6(a), not only was the PSO convergence rate significantly improved after handling more constraints, but the number of function evaluations (FE) was reduced from p1 to p3. For the GA, the results were improved after coupling the BU method with the constraint handling scheme. After handling more constraints using the BU method, the results were slightly improved. For the DE algorithm, the results showed similar trends, when more constraints are handled, a better convergence rate is achieved. To better assess BU performance, the percentage of constraint violations during iterations were tracked, and the results are shown in Fig. 6(d) to (f) for the PSO, GA, and DE algorithms, respectively. Based on the results, it can be seen that with increasing the number of handled constraints using the BU method, a smaller percentage of constraint violations is obtained. This shows that the BU method eliminates more infeasible search space by handling more constraints. An interesting point is that after handling all the constraints by the BU method (3p), no violation occurred for this problem, which has nine linear constraints. Therefore, the BU method changed a problem with a hit ratio of 0.0111% into an unconstrained problem. In other words, when all the constraints are considered in this process, the BU method eliminated 99.9889% of the search space, which are infeasible regions from the first iteration.

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

13

Table 1 Statistical results of the different methods on the G1 problem (Columns 4-7 indicate statistics of objective function values and the best performing method is marked in bold) Method

BU Step

Ave No. FE

Best

Mean

Median*

Worst

St. Dev.

IP

Nonea Step 1 (1p) Step 2 (2p) Step 3 (3p)

527 557 666 429

−14.7098 −14.6176 −14.8938 −15

91.57562 N/A 97.82933 −13.6044

−4.4103343 N/A −6.9422484 −13.589442

N/Ab N/A N/A −10.4378

N/A N/A N/A 1.432381

SQP

None Step 1 (1p) Step 2 (2p) Step 3 (3p)

40 41 62 131

−1.80E−5 −2.00002 −15 −15

−2.40E−6 −0.35294 26.13126 −13.8995

0 −1.78E−15 0 −13.828125

0 0 N/A −11.8281

5.65E−6 0.770028 N/A 1.171343

DS

None Step 1 (1p) Step 2 (2p) Step 3 (3p)

1670 1714 1924 1904

−11.7281 −15 −15 −15

N/A N/A N/A −14.7539

N/A N/A 13.347739 −14.999956

N/A N/A N/A −12.4531

N/A N/A N/A 0.684101

SA

None Step 1 (1p) Step 2 (2p) Step 3 (3p)

20,000 20,000 20,000 20,000

−15 −15 −15 −14.2721

−14.1886 −13.2305 −13.3812 −11.4669

−15 −15 ↓ −15 ≈ −11.472473 ↓

−8.98457 −8.48028 −9.45394 −9.77542

1.802651 2.576221 2.183331 0.823801

PSO

None Step 1 (1p) Step 2 (2p) Step 3 (3p)

5300 4150 2450 1750

−15 −15 −15 −15

−9.62745 −11.5882 −12.8235 −14.9412

−9.0000001 −12 ↑ −12 ↑ −15 ↑

−6 −9 −12 −12

2.591992 1.991746 1.352122 0.420084

GA

None Step 1 (1p) Step 2 (2p) Step 3 (3p)

16,750 16,900 16,850 17,900

−15 −15 −15 −15

−5.3217 −14.2595 −14.7541 −14.9976

−14.999351 −14.99952 ↑ −14.999642 ↑ −14.999963 ↑

N/A −10.3225 −12.053 −14.9762

N/A 1.248775 0.760859 0.004485

DE

None Step 1 (1p) Step 2 (2p) Step 3 (3p)

5150 6250 5900 5450

−13.8496 −14.8215 −15 −15

−10.2621 −12.8151 −14.1068 −14.5677

−10.529502 −13.083049 ↑ −14.785518 ↑ −14.942472 ↑

−4.59663 −9 −11.8068 −11.1851

1.851014 1.581668 1.153452 0.799997

a None

means BU is not used and the explicit constrained handling scheme (here static penalty method) works alone. means not available, as the algorithm is not converged *The Wilcoxon rank-sum test at a significance level of 5% was performed to compare the results obtained by different approaches. ↑, ≈ and ↓: The approach performs statistically significantly better, equal, and worse than the former approach. b N/A

3.1.2. Non-linear constraints: Himmelblau’s problem After solving the G1 problem with linear constraints, the BU method was benchmarked using a well-known problem with six nonlinear inequality constraints. This problem was initially proposed by Himmelblau [11], and it has been widely used as a benchmark nonlinear constrained optimization problem. In this problem, there are five design variables with ten defined boundary conditions. In this problem, two constraints are active at the global optimum, and this optimization problem has a hit ratio equal to 51.1230%, which means approximately half of the search space is in the infeasible region. The fifth variable (x5 ) is the repairing variable and its bounds were updated using the six nonlinear constraints. The Himmelblau problem was optimized by using the seven optimization algorithms and the results are presented in Table 2, both with and without the BU method. For the IP method, the results were improved after using the proposed approach. Results of the SQP method show that by using only the proposed approach, all the runs can find feasible solutions. Results of all other methods (DS, SA, PSO, GA, and DE algorithms) were also improved after using the proposed approach with the nonlinear constraints for this case study. Based on the Wilcoxon rank-sum test, the results of population-based algorithms (PSO, GA, and DE algorithms) were significantly improved after using BU approach. The convergence histories of the best solutions for the population-based algorithms (PSO, GA, and DE) are presented in Fig. 7(a) to (c). It can be seen that the convergence histories were improved in all cases after using the BU approach. It can be seen that, although the PSO algorithms finally converged to the same value, the convergence

14

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Fig. 6. Convergence histories of (a) PSO, (b) GA, and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA, and (f) DE algorithms with and without the BU approach for G1 problem.

Fig. 7. Convergence histories of (a) PSO, (b) GA, and (c) DE algorithms and percentage of constraint violation histories of (d) PSO, (e) GA, and (f) DE algorithms with and without the BU approach for Himmelblau’s problem.

15

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917 Table 2 Statistical results of the different methods used on Himmelblau’s problem (Columns 4-7 indicate statistics on objective function values) Method

BUa

Ave No. FE

Best

Mean

Median*

Worst

St. Dev.

IP

No Yes

151 175

−30,372.344 −30,665.276

N/A N/A

−28,999.384 −29,671.289

N/A N/A

N/A N/A

SQP

No Yes

93 135

−30,648.680 −30,514.974

N/A −27,779.518

−26,010.718 −29,160.920

N/A −23,068.751

N/A 2598.321

DS

No Yes

327 467

−29,990.349 −30,213.624

−29,129.981 −29,281.557

−29,135.778 −29,548.490

−28,242.639 −27,461.713

424.846 651.177

SA

No Yes

10,000 10,000

−30,586.169 −30,659.787

−30,207.756 −30,225.788

−30,302.689 −30,232.027 ≈

−29,441.442 −29,672.082

313.040 301.932

PSO

No Yes

5440 2440

−30,665.538 −30,665.539

−30,665.348 −30,665.538

−30,665.529 −30,665.539 ↑

−30,661.131 −30,665.534

0.804 0.001

GA

No Yes

4880 4880

−30,638.021 −30,660.049

−30,033.872 −30,411.344

−30,085.975 −30,481.412 ↑

−29,371.713 −29,813.671

305.275 231.796

DE

No Yes

2800 3040

−30,660.41 −30,665.54

−30,540.543 −30,660.433

−30,564.525 −30,664.733 ↑

−30,221.391 −30,624.297

100.929 9.180

a In

this column “No” means the BU method was not used and the explicit constrained handling scheme (here static penalty method) worked alone, and “Yes” means the BU method is coupled with the explicit constrained handling scheme. *The Wilcoxon rank-sum test at a significance level of 5% was performed to compare the results obtained by the approach with and without the BU approach.↑, ≈ and ↓: Method with the BU approach performed statistically significantly better, equal, and worse.

rate of the PSO algorithm with BU was improved and PSO-BU needed a lower number of FE to reach the final solution. The improvements of the GA and DE algorithms were more significant in this optimization problem, as they converged to a better value after using the BU in the optimization process. The percentages of the constraint violations during iterations are presented in Fig. 7(d) to (f) for the PSO, GA, and DE algorithms, respectively. From the results, it can be seen that the BU approach significantly decreased the constraint violations in both the PSO and DE algorithms. From Fig. 7(b), it can be seen that the percentage of violations is almost the same before and after using the BU in GA. However, as shown in Fig. 7(b), the convergence history of GA improved significantly, even with no particular changes in the percentage of constraint violations. Comparing this problem with the former problem, it can be seen that after implicitly handling all the constraints in the Himmelblau’s problem, there are still some constraint violations. This is because cutting the entire infeasible region cannot be done in one step; it is an iterative process for such problems with nonlinear constraints. 3.2. Constrained engineering problems Eight constrained engineering problems were optimized in this section. The first three problems involved pressure vessel design, welded beam design, and tension–compression spring design, which are well-known benchmark engineering constrained problems that have been used by many researchers to show their algorithms’ capability of handling the constraints. The other problems, including speed reducer design, corrugated bulkheads design, rolling element bearing design, and multiple disk clutch brake design, are constrained engineering problems that have been previously solved in the literature. General information about variables and constraints for these problems is presented in Table 3. All of these problems involve minimization except for the last problem, which is a maximization problem. The detailed formulation of these problems can be found in the Appendix. In addition to the general parameter settings mentioned for each optimization algorithm, population size and the maximum number of iterations should be set for these algorithms. The last two columns of Table 3 show these numbers for the selected problems. 3.2.1. Welded beam design problem Minimizing the overall cost of welded beam fabrication is the objective of this well-known benchmark. This problem has four variables: the thickness of the weld, the length of the welded joint, the width of the beam, and

16

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Table 3 Main features of the constrained problems and basic setting of the optimization algorithms for the problems. Problem name

Problem characterization

Welded beam Pressure vessel Tension–compression spring Speed reducer Corrugated bulkheads design Heater exchanger Multiple disk clutch brake Rolling element bearing a LC

Optimization setting

No. Var.

Var. Type

No.

4 4 3 7 4 8 5 10

Continuous Mixed Continuous Continuous Continuous Continuous Discrete Mixed

4 3 1 2 4 3 7 7

LCa

No.

NLCb

1 1 3 7 2 3 1 4

Pop. Size

Max No. Iter.

40 50 50 40 40 100 40 25

250 200 300 200 250 1000 25 50

is Linear Constraint. is Non-linear Constraint.

b NLC

Table 4 Statistical results of the different methods used on the welded beam design problem (Columns 4-7 indicate statistics on objective function values) Method

BU

Ave No. FE

Best

Mean

Mediana

Worst

St. Dev.

IP

No Yes

208 275

2.867425 2.512943

N/A N/A

N/A 5.92662

N/A N/A

N/A N/A

SQP

No Yes

235 196

4.96 2.485397

N/A N/A

N/A 5.86

N/A N/A

N/A N/A

DS

No Yes

479 459

3.160056 2.720765

N/A N/A

7.2276645 6.5092949

N/A N/A

N/A N/A

SA

No Yes

10,000 10,000

2.412406 2.41625

2.602433 2.540633

2.5631256 2.5342211 ≈

3.295827 2.865976

0.161404 0.090123

PSO

No Yes

9000 6360

2.381136 2.381134

2.63765 2.381139

2.4132258 2.3811353 ↑

5.989648 2.381245

0.571996 1.56E−5

GA

No Yes

5640 5480

2.59891 2.388155

N/A 2.985064

4.767209 2.8597056 ↑

N/A 4.659896

N/A 0.511872

DE

No Yes

2560 3480

2.55119 2.381612

3.680576 2.498382

3.3646944 2.4456959 ↑

7.915949 3.128179

1.056895 0.152148

a ↑,

≈ and ↓: Method with BU approach performs statistically significantly better, equal, and worse.

the thickness of the beam. All six constraints of this problem were handled by updating the boundaries of beam thickness. The hit ratio of this problem is 2.686%, and the BU approach represents an effort to cut the remaining 97.314% (the infeasible search space) for this problem. The seven optimization algorithms were applied to this problem and the results are presented in Table 4. From this table, it can be seen that most of the IP and SQP runs did not find even a single feasible solution if they did not use the proposed approach. After using BU in IP and SQP algorithms, not only did most of them reach the feasible area but, also, they found more feasible solutions. Results of the DS algorithm show that the BU approach can improve the solutions. For the SA algorithm, all statistics were improved, except for the best solution of all runs. Based on the Wilcoxon rank-sum test, the results of population-based algorithms (PSO, GA, and DE algorithms) were significantly improved after using BU approach. Statistical results of the population-based algorithms PSO, GA, and DE were significantly improved. Based on the GA results, some runs did not find even a single feasible solution if the BU was not integrated into this algorithm. The convergence histories of the population-based algorithms are shown in Fig. 8. Fig. 8(a) to (c) show the convergence histories of the PSO, GA, and DE algorithms, respectively. For the PSO and DE algorithms, using the BU approach significantly improved the algorithms’ convergence rate and the final results. Although the number of FE in the DE is less than the DE with the BU approach, it converged to a larger number for this minimization

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

17

Fig. 8. Convergence histories of (a) PSO, (b) GA, and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA, and (f) DE algorithms with and without the BU approach for the welded beam design problem.

problem. For the GA results, it should be mentioned that one out of the 51 runs was not able to find the feasible region. Therefore, this run’s convergence history is not considered here. Therefore, the GA may not be able to find a feasible area without the BU approach. Based on the GA convergence results (Fig. 8(b)), even after removing the infeasible run history, GA-BU performed much better than the GA. Fig. 8(d) to (f) show the percentage of constraint violation histories for the PSO, GA, and DE algorithms, respectively. Based on these three figures, the proposed BU approach effectively decreased the constraint violations in this problem. 3.2.2. Pressure vessel design problem In this design problem, the cylindrical pressure vessel is capped at both ends by hemispherical heads. This problem is a mixed variable problem and the total cost should be minimized by satisfying four constraints: three linear and one nonlinear. The problem variables are thickness, thickness of the head, the inner radius, and the length of the cylindrical section of the vessel. The hit ratio is equal to 39.676% for this problem, and the boundaries of the inner radius and the length of the cylindrical section of the vessel were updated to handle the four constraints. The problems solved by the IP, SQP, DS, SA, PSO, GA, and DE algorithms with and without the BU approach and the statistical results are presented in Table 5. The IP results were significantly improved; however, the IP methods cannot converge to good solutions and they are both trapped in local optima. The results of the SQP algorithm show that, although the median of the final results was improved, the best solution got worse. It should be noted that the required number of FEs in SQP decreased after using the proposed approach. For this case study, the DS results became worse after implementing the BU approach. Implementing the BU approach in the SA algorithm was very effective, and all runs found feasible and better solutions. Again, it can be seen that all the statistical results of these three population-based algorithms were notably improved. Among these algorithms, the average number of FEs in the PSO algorithm is decreased by approximately 70% after using the proposed approach. Based on the Wilcoxon rank-sum test, the results of metaheuristics were significantly improved after using the BU approach.

18

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Fig. 9. Convergence histories of (a) PSO, (b) GA, and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA, and (f) DE algorithms with and without the BU approach for pressure vessel design problem.

Table 5 Statistical results of the different methods on the pressure vessel design problem. Method

BU

Ave No. FE

Best

Mean

Mediana

Worst

St. Dev.

IP

No Yes

222 110

11,149.38 8265.363

N/A N/A

54,226.344 41,382.257

N/A N/A

N/A N/A

SQP

No Yes

153 57

9112.831 7026.748

N/A N/A

48,953.879 69,917.424

N/A N/A

N/A N/A

DS

No Yes

440 395

6434.457 6059.715

7397.758 N/A

7421.4887 6059.7195

8025.892 N/A

384.527 N/A

SA

No Yes

10,000 10,000

6079.8 6059.714

N/A 8886.518

6856.1786 6059.7143 ↑

N/A 21,929.2

N/A 4629.473

PSO

No Yes

7100 2200

6059.714 6059.714

6479.912 6086.195

6371.1388 6059.7143 ↑

7544.493 7410.217

445.0283 189.1082

GA

No Yes

5450 5450

6103.559 6059.714

6868.061 6088.245

6873.7118 6059.7222 ↑

7573.466 6633.037

434.9137 104.8163

DE

No Yes

3200 3350

6182.118 6059.714

8450.549 6311.195

7696.4462 6059.7143 ↑

14,890.34 8530.48

2019.808 474.182

a ↑,

≈ and ↓: Method with the BU approach performed statistically significantly better, equal, and worse.

The convergence histories of the population-based algorithms are presented in Fig. 9. The convergence results of the PSO, GA, and DE algorithms for this problem are presented in Fig. 9(a) to (c), respectively. From the convergence plots, it is clear that not only does the BU approach improve the convergence rate, but the BU approach also guides the population-based optimization algorithms to better solutions that may not be reached without the BU approach.

19

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Table 6 Statistical results of the different methods used on the tension–compression spring design problem (Columns 4-7 indicate statistics on objective function values) Method

BU

Ave No. FE

Best

Mean

Mediana

Worst

St. Dev.

IP

No Yes

311 222

0.012949 0.012754

N/A N/A

0.030951 0.022764

N/A N/A

N/A N/A

SQP

No Yes

15 127

0.175088 0.012666

N/A N/A

N/A 0.030084

N/A N/A

N/A N/A

DS

No Yes

510 200

0.013161 0.012759

N/A 0.030254

N/A 0.022194

N/A 0.061299

N/A 0.017098

SA

No Yes

15,000 15,000

0.012685 0.012671

N/A 0.012817

0.012845 0.012812 ↑

N/A 0.013218

N/A 0.000114

PSO

No Yes

3950 3050

0.012666 0.012665

N/A 0.012684

0.012719 0.012677 ↑

N/A 0.012734

N/A 0.000019

GA

No Yes

10,950 3600

0.012674 0.012672

N/A 0.013301

0.015008 0.013199 ↑

N/A 0.015827

N/A 0.000642

DE

No Yes

2500 2400

0.012801 0.012675

0.013755395 0.012877987

0.01334 0.012831 ↑

0.017706 0.013246

0.000951 0.000168

a ↑,

≈ and ↓: Method with the BU approach performed statistically significantly better, equal, and worse.

Fig. 9(d) to (f) show the percentage of constraint violation histories for the population-based optimization algorithms. For the PSO and GA algorithms, the advantages and reduced constraint violations are clear. For the DE algorithm, although there is a significant improvement in the convergence rate and final solution results, the percentage of constraint violations is almost the same during iterations for the algorithms with and without the BU approach. 3.2.3. Tension–compression spring design problem The tension–compression spring design is another well-known benchmark constrained problem. This engineering problem has three design variables – wire diameter, the mean coil diameter, and the number of active coils – that minimize the weight of the spring. This problem has three nonlinear and one linear constraint, which make the problem highly constrained, with a hit ratio of 0.754%. Therefore, in this problem, even finding a feasible region is a major achievement. All of the constraints are handled by updating the bounds of the first variable, wire diameter. The problems solved by the PSO, GA, and DE algorithms, with and without the BU approach, and statistical results are presented in Table 6, which shows that the final statistical results of all algorithms were improved. Most SQP runs and all DS, SA, PSO, and GA runs were able to find the feasible area only by using the proposed BU approach. Although the DE algorithm could find at least a single feasible solution in each run, the ultimate solutions were improved. For GA, the average number of FEs was decreased by approximately 30% after using the proposed approach. Based on the Wilcoxon rank-sum test, the results of metaheuristics were significantly improved after using the BU approach. The convergence histories of the population-based algorithms are presented in Fig. 10. The convergence of the best solutions of the PSO, GA, and DE algorithms for this problem are presented in Fig. 10(a) to (c), respectively. First, it should be noted that PSO and GA had 4 and 6 infeasible runs, respectively, which were removed from the convergence plots. These infeasible runs show that GA and PSO may not be able to find a feasible area without the BU approach. The convergence history results of the PSO and GA algorithms are interesting, in that their results are still far worse than the PSO-BU and GA-BU results, even without considering their infeasible runs. This shows that the BU not only found a feasible area in all cases but it also helps the optimization process to find better solutions. For the DE algorithms, the proposed BU approach is still a major improvement in terms of the convergence rate and the ability to reach a better solution. Based on the convergence plots, not only did the BU approach improve the convergence rate but also guided the population-based optimization algorithms to better solutions that may not have been reached without the BU approach.

20

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Fig. 10. Convergence histories of (a) PSO, (b) GA and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA and (f) DE algorithms with and without BU approach for tension–compression spring design problem.

Fig. 10(d) to (f) show the percentage of constraint violation histories for the population-based optimization algorithms. From these plots, it is clear that after using the BU method, the percentage of constraint violations decreased significantly for all PSO, GA, and DE algorithms. 3.2.4. Corrugated bulkheads design problem The objective of this problem is the minimization of weight design of the corrugated bulkheads for a tanker with varying width, depth, length, and plate thickness. All six defined constraints for this problem were handled using the proposed approach. For this purpose, the lower bounds of the length and thickness were updated using the six constraints. The statistical results of the seven optimization algorithms for this case study are presented in Table 7. From the table, it is clear that IP, SQP, and SA found a feasible solution in all runs only if they use the proposed approach. For this case study, statistical results were significantly improved for all optimization algorithms after using the proposed approach for implicitly handling the constraints. The convergence and constraint violation histories of the three population-based optimization algorithms are presented in Fig. 11. Fig. 11(a) shows that the convergence rate of the PSO algorithm was improved along with a slight improvement in the final results. However, for the GA and DE algorithms (Fig. 11(b) and (c), the BU approach significantly improved both convergence rates and final results. The percentage of constraint violations histories for the PSO, GA, and DE algorithms are shown in Fig. 11(d) to (f), respectively. Based on these figures, the proposed approach can handle all the constraints at the early iterations and can keep all solutions in the feasible region afterward. 3.2.5. Speed reducer design problem The speed reducer design problem has seven variables: the face width, the module of teeth, number of teeth on the pinion, lengths of the first and second shafts between bearings, and diameters of the first and second shafts. Minimizing the total weight of the speed reducer is the objective in this design problem that is subject to eight constraints on the stresses and deflections. To handle the constraints using the proposed approach, three repairing variables were used as the face width and the two shaft diameters. Both lower and upper bounds of these variables were updated in this problem to reduce the entire infeasible search space.

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

21

Table 7 Statistical results of the different methods on the corrugated bulkheads design problem (Columns 4-7 indicate statistics on objective function values) Method

BU

Ave No. FE

Best

Mean

Mediana

Worst

St. Dev.

IP

No Yes

197 505

7.995253 5.886933

N/A 8.501321

15.643122 5.8907472

N/A 45.4883

N/A 8.071991

SQP

No Yes

199 298

6.81 5.886732

N/A 8.079758

12.136814 5.89

N/A 20.1267

N/A 4.605882

DS

No Yes

343 804

8.527781 5.887325

11.19954 5.931027

10.469909 5.9183374

18.64175 6.037309

2.01554 0.039332

SA

No Yes

10,000 10,000

6.277661 5.965082

N/A 6.505472

7.8439923 6.3442512 ↑

N/A 8.839613

N/A 0.553489

PSO

No Yes

8120 4440

5.886888 5.886733

5.937081 5.887928

5.9326612 5.8873569 ↑

6.022314 5.891192

0.031512 0.001211

GA

No Yes

2840 3680

5.944572 5.886907

6.827663 5.961474

6.7240928 5.9172349 ↑

8.597889 6.296148

0.653139 0.103952

DE

No Yes

3320 4000

5.916377 5.88676

6.710859 5.937565

6.4460257 5.9110398 ↑

9.145381 6.330802

0.821435 0.083645

a ↑,

≈ and ↓: Method with BU approach performs statistically significantly better, equal, and worse.

Fig. 11. Convergence histories of (a) PSO, (b) GA, and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA, and (f) DE algorithms with and without the BU approach for corrugated bulkheads design problem.

IP, SQP, DS, SA, PSO, GA, and DE algorithms were applied to the optimum design of the speed reducer with and without using the proposed approach, and the statistical results are presented in Table 8, which shows that all the methods can find the feasible region. It can also be seen that all statistical results of the final results were improved, except the best and worst solutions of the SQP method. Based on the Wilcoxon rank-sum test, the results of metaheuristics are significantly improved after using BU approach.

22

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Table 8 Statistical results of the different methods on the speed reducer design problem (Columns 4-7 indicate statistics on objective function values) Method

BU

Ave No. FE

Best

Mean

Mediana

Worst

St. Dev.

IP

No Yes

333 244

3038.491 2994.472

3554.457 3256.598

3298.1091 3066.3287

5636.63 5277.263

584.2696 462.1228

SQP

No Yes

230 118

3000.759 3003.014

3179.043 3120.643

3227.9257 3007.442

3561.528 4589.43

112.864 324.7745

DS

No Yes

1531 1159

2994.472 2994.472

2994.482 2994.48

2994.4787 2994.4779

2994.495 2994.555

0.007145 0.011272

SA

No Yes

8000 8000

3004.615 2994.491

3268.775 2997.866

3129.1346 2996.7567 ↑

4427.868 3013.968

363.4283 3.76164

PSO

No Yes

5760 2040

2994.471 2994.471

2995.496 2994.727

2994.4717 2994.4712 ↑

3033.749 3007.45

5.757177 1.82

GA

No Yes

7960 7200

2994.72 2994.479

2996.395 2995.642

2996.2411 2994.9535 ↑

3000.309 3009.629

1.20314 2.485562

DE

No Yes

3760 4160

2996.762 2994.471

3022.575 2994.654

3012.8139 2994.4943 ↑

3091.73 2997.492

23.63895 0.527972

a ↑,

≈ and ↓: Method with BU approach performs statistically significantly better, equal and worse.

Fig. 12. Convergence histories of (a) PSO, (b) GA, and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA, and (f) DE algorithms with and without the BU approach for speed reducer design problem.

Among the population-based algorithms, the average number of FEs in the PSO algorithm decreased by approximately 65% after using the proposed approach. The convergence and constraint violation histories of the three population-based optimization algorithms are illustrated in Fig. 12. The convergence results show that coupling the BU with the population-based optimization algorithm notably increased the convergence rate of the algorithms for this problem and all the algorithms quickly converged to the global maximum after using the proposed BU approach.

23

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Table 9 Statistical results of the different methods used on the heat exchanger design problem (Columns 4-7 indicate statistics on objective function values) Method

BU

Ave No. FE

Best

Mean

Mediana

Worst

St. Dev.

IP

No Yes

566 728

18 913.31 8954.532

N/A N/A

N/A N/A

N/A N/A

N/A N/A

SQP

No Yes

336 290

N/A 12 404.07

N/A N/A

N/A N/A

N/A N/A

N/A N/A

DS

No Yes

1687 1092

N/A 7278.164

N/A N/A

N/A 10 110.92

N/A N/A

N/A N/A

SA

No Yes

100 000 100 000

7104.19 7076.451

N/A N/A

10 351.832 8130.68 ≈

N/A N/A

N/A N/A

PSO

No Yes

42 300 22 600

7123.604 7056.142

N/A 7240.345

N/A 7250.602 ↑

N/A 7626.203

N/A 86.434

GA

No Yes

97 900 100 000

7715.965 7374.931

N/A 9693.205

N/A 9587.343 ↑

N/A 1.26E+4

N/A 1209.594

DE

No Yes

5400 6200

12 592.53 8036.613

N/A 10 134.94

N/A 9792.719 ↑

N/A 15 197.8

N/A 1613.002

a ↑,

≈ and ↓: Method with BU approach performs statistically significantly better, equal and worse.

Percentage of constraint violation histories for the population-based optimization algorithms are shown in Fig. 12(d) to (f). From these figures, it can be seen that the percentage of constraint violations decreased significantly in all PSO, GA, and DE algorithms after they were coupled with the proposed approach. 3.2.6. Heat exchanger design problem In this problem, all six constraints are binding, and the hit ratio is equal to 0.001%. Therefore, the heat exchanger design problem is one of the most challenging constrained problems. Of the eight variables in this problem, only three contribute to the objective function, but all of them contribute to the constraints. The fourth and eighth variables were selected as the repairing variables in the BU approach, and their lower and upper bounds were updated using the constraints to reduce the infeasible region, which was approximately 99.999% of the entire search space. The seven algorithms were applied to the optimal design of the speed reducer, with and without using the proposed approach, and the statistical results are presented in Table 9, which shows that all results were improved after using the BU approach in these algorithms. It can also be seen that finding a single feasible solution is difficult for most methods if they do not use the proposed approach, which is expected, as the problem has a very small feasible region. Based on Table 9, SQP and DS methods were stuck in the infeasible region in all their simulations if they failed to use the proposed constraint handling scheme. As shown in Table 9, the population-based algorithm without the proposed approach could not find a single feasible solution in most runs. However, using the proposed approach, all their runs could find feasible solutions. Table 9 shows that this approach is effective for all the seven optimization algorithms and particularly in population-based algorithms. Based on the Wilcoxon rank-sum test, the results of population-based metaheuristics were significantly improved by the BU approach. The best solutions for population-based algorithms are presented in Fig. 13(a) to (c). It should be mentioned that all the algorithms without BU have several infeasible runs. In 52 runs, the number of infeasible runs was 18, 14, and 15 for the PSO, GA, and DE algorithms, respectively. Therefore, it is clear that all the algorithms are associated with a small probability of finding a feasible region when they do not use the BU approach. To visualize the convergence results, all the infeasible runs were excluded from the convergence plots. In all three convergence histories, the population-based algorithms with the BU approach can always find the feasible region very quickly. Furthermore, the algorithms using the BU always converge faster and to a better solution as compared to feasible runs of the algorithms without BU. Fig. 13(d) to (f) respectively show the percentage of constraint violations histories for the PSO, GA, and DE algorithms. For the PSO and GA algorithms, the advantages and reduced constraint violations are clear. However, for

24

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Fig. 13. Convergence histories of (a) PSO, (b) GA, and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA, and (f) DE algorithms with and without the BU approach for the heat exchanger design problem. Table 10 Statistical results of the different methods on the multiple disk clutch brake design problem (Columns 4-7 indicate statistics on objective function values) Method

BU

Ave No. FE

Best

Mean

Mediana

Worst

St. Dev.

DS

No Yes

207 232

0.317577 0.313657

0.367132 0.320483

0.352864 0.313657

0.441080 0.406773

0.037786 0.020226

SA

No Yes

1000 1000

0.313657 0.313657

0.388835 0.331623

0.380995 0.325419

0.537235 0.409101 ↑

0.062294 0.021566

PSO

No Yes

800 575

0.313657 0.313657

0.318195 0.313657

0.313657 0.313657

0.339999 0.313657 ↑

0.009371 0.000000

GA

No Yes

1000 1000

0.313657 0.313657

0.364087 0.335573

0.345022 0.317577

0.536034 0.476366 ↑

0.051066 0.035331

DE

No Yes

950 875

0.313657 0.313657

0.322554 0.313810

0.313657 0.313657

0.345022 0.317577 ↑

0.010696 0.000769

a ↑,

≈ and ↓: Method with BU approach performs statistically significantly better, equal, and worse.

the PSO algorithm, the percentage of constraint violations drops to zero, but the percentage of constraint violations drops to about half of the population, and this percentage remains nearly constant until the end of runs. For the DE algorithm, although there is a significant improvement in the convergence rate and final solution results (Fig. 13(c)), the percentage of constraint violations is not improved after using the BU approach. 3.2.7. Multiple disk clutch brake design problem This problem is a minimization optimization problem with five discrete design variables: inner and outer radii, the thickness of the disc, actuating force, and the number of friction surfaces. The BU approach is applied here to cut the infeasible solutions by updating the upper bound of the inner radius and updating the upper and lower bounds of the actuating force and number of friction surfaces. The IP and SQP methods performed poorly and are quickly trapped in a local minimum (average number of FE is equal to 6) as they are discrete optimization algorithms, and, therefore, their results are not presented. This discrete

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

25

Fig. 14. Convergence histories of (a) PSO, (b) GA, and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA, and (f) DE algorithms with and without the BU approach for multiple disk clutch brake design problem.

optimization problem is solved by using DS, SA, PSO, GA, and DE optimization algorithms, with and without using the proposed approach, and the statistical results are presented in Table 10. From this table, all statistical results of the final results for the seven optimization methods were improved after implementing the BU method. Based on the Wilcoxon rank-sum test, the results of metaheuristics are significantly improved after using the BU approach. The convergence results of the PSO, GA, and DE algorithms are presented in Fig. 14(a) to (c), respectively. As shown in the figures, the proposed approach effectively cut the infeasible search space in the population-based optimization algorithms. Fig. 14(d) to (f) show the percentage of constraint violations histories for the populationbased optimization algorithms. From these plots, it is clear that after using the BU method, the percentage of constraint violations decreases significantly in all PSO, GA, and DE algorithms. Based on Fig. 14(f), the DE performs the whole searches in a feasible area after it couples with the proposed BU approach. 3.2.8. Rolling element bearing design problem This problem involves the maximization of the dynamic load-carrying capacity of a rolling element bearing. This problem has 10 variables: pitch and ball diameters, a number of balls, and inner and outer raceway curvature coefficients. The problem has 11 constraints and the BU is applied to handle them. For this purpose, four repairing variables were used as pitch and ball diameters, a number of balls and a coefficient for outer ring strength consideration. The upper bounds of the pitch diameter and the coefficient and both bounds of the ball diameters and number of balls are updated. The seven algorithms are applied to the optimum design of the speed reducer with and without using the proposed approach, and the statistical results are presented in Table 11. First, it can be seen that all the statistical results of the seven optimization methods were improved after implementing the BU method for this maximization problem. From the IP and SQP method results, it can be seen that most of these runs cannot find the feasible region if the BU approach is not used along with these algorithms. However, after using the proposed approach in both IP and SQP methods, not only can all the runs reach the feasible region, but they also find better fitness values for this maximization problem. From Table 11, it can also be seen that SA and GA can find a feasible area in all runs only

26

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Table 11 Statistical results of the different methods on rolling element bearing design problem (Columns 4-7 indicate statistics on objective function values) Method

BU

Ave No. FE

Best

Mean

Mediana

Worst

St. Dev.

IP

No Yes

412 700

18 713.95 82 339.84

N/A 52 975.23

N/A 52 124.48

N/A 19 737.79

N/A 18 753.11

SQP

No Yes

338 88

26 271.23 82 113.3

N/A 56 152

N/A 58 085.23

N/A 20 416.98

N/A 14 762.44

DS

No Yes

1119 866

85 491.56 81 236.88

78 469.92 69 415.03

77 231.89 76 477.45

70 159.9 41 567.4

4074.62 11 923.85

SA

No Yes

1250 1250

73 186.62 85 533.15

N/A 83 835.11

35 239.82 85 241.86

N/A 75 061.51 ↑

N/A 2.84E+3

PSO

No Yes

1250 1200

85 526.12 85 533.18

79 821.99 85 464.94

84 459.58 85 533.18

43 553.74 84 459.13 ↑

11 248.22 255

GA

No Yes

1250 1250

63 448.96 84 669.21

N/A 74 622.69

32 449.49 76 185.59

N/A 58 852.9 ↑

N/A 7531.45

DE

No Yes

1175 1175

84 171.18 85 463.66

76 092.67 84 527.23

78 213.85 84 454.4

43 239.9 82 352.89 ↑

7367.43 577.38

a ↑,

≈ and ↓: Method with BU approach performs statistically significantly better, equal, and worse.

Fig. 15. Convergence histories of (a) PSO, (b) GA, and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA, and (f) DE algorithms with and without BU approach for rolling element bearing design problem.

if they use the proposed approach. Based on the Wilcoxon rank-sum test, the statistical results of metaheuristics are significantly improved after using BU approach. The convergence results of the population-based methods are shown in Fig. 15. The convergence histories of this maximization problem for the PSO, GA, and DE algorithms are shown in Fig. 15(a) to (c), respectively. From the results, it can be seen that by using the BU approach, both convergence rates and the final results of all populationbased algorithms were significantly improved. It should be noted that five of 51 runs of GA without BU did not

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

27

Fig. 16. Schematic of car side design problem [24].

even converge, which means they are excluded from Fig. 15(b). All these results confirm that the BU approach notably enhances the results of the population-based algorithm for this problem. Fig. 15(d) to (f) show the percentage of constraint violations histories for the PSO, GA, and DE algorithms, respectively. For GA and DE algorithms, the advantages and reduced constraint violations are clear. However, for the PSO algorithm with the BU approach, the percentage of constraint violation dropped quickly in early iterations; however, it remains constant afterward and will converge to the same constraint violations of the PSO without the BU approach. 3.3. Surrogate model optimization Sometimes the constraint(s), such as black-box constraints, are so complex that they cannot be solved by one variable. Additionally, evaluating the actual constraints can be very time-consuming. One effective and popular method to reduce the complexity and simulation time with such problems, especially to reduce the simulation time, is using surrogate models instead of the actual constraint. If a considered surrogate model can be solved using the one variable, then the BU approach could be applied to handle that surrogate model of the constraint. Therefore, by choosing a solvable surrogate model, in addition to replacing a complex and tedious constraint with a simpler model, the BU approach can be applied to the problem. The car side impact problem is one of the most complex and time-consuming problems (shown in Fig. 16). The objective function is a minimization of the car weight, which is a summation function. However, computing constraints need a full-vehicle finite element analysis (FEA) model, an FEA side impact dummy model, and an FEA deformable side impact barrier model [24]. This problem is a mixed variable problem with nine variables: thicknesses of B-Pillar inner, B-Pillar reinforcement, floor side inner, cross members, door beam, door beltline reinforcement and roof rail, materials of B-Pillar inner and floor side inner and barrier height and hitting position. The problem formulation is presented as follows: Minimize f(x) = Weight; Subject to g1 (x) = Fa (load in abdomen) ≤ 1 kN g2 (x) = V × Cu (dummy upper chest) ≤ 0.32 m/s g3 (x) = V × Cm (dummy middle chest) ≤ 0.32 m/s g4 (x) = V × Cl (dummy lower chest) ≤ 0.32 m/s g5 (x) = ∆ur (upper rib deflection) ≤ 32 mm g6 (x) = ∆mr (middle rib deflection) ≤ 32 mm g7 (x) = ∆lr (lower rib deflection) ≤ 32 mm g8 (x) = Fp (Pubic force) ≤ 4 kN g9 (x) = VMBP (Velocity of V-Pillar at middle point) ≤ 9.9 mm/ms

28

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Table 12 Statistical results of the different methods used on the car side design problem (Columns 4-7 indicate statistics on objective function values) Method

BU

Ave No. FE

Best

Mean

Mediana

Worst

St. Dev.

IP

No Yes

406 669

24.17905 24.26806

29.76170 N/A

29.38639 26.89433

39.32257 N/A

3.77191 N/A

SQP

No Yes

205 297

26.74241 24.18723

32.58760 28.61597

31.51611 28.29241

37.46901 36.66473

3.88078 3.58368

DS

No Yes

1420 1918

23.56836 24.20571

27.55274 25.76193

27.79711 25.30845

31.20660 29.48342

1.69953 1.28207

SA

No Yes

8000 8000

23.89347 23.83496

24.78709 24.51510

24.67856 24.42886 ↑

25.91154 25.83053

0.43392 0.41689

PSO

No Yes

7240 6880

23.46559 23.46651

23.56380 23.54601

23.54421 23.54705 ≈

24.13825 23.59425

0.10109 0.03675

GA

No Yes

7960 7920

23.56262 23.60405

24.66953 24.07140

24.36477 23.90356 ↑

27.29189 25.28348

0.83735 0.45051

DE

No Yes

3720 3280

23.78042 23.67979

24.39080 24.36401

24.25911 24.23688 ≈

26.29569 25.92218

0.54243 0.53612

a ↑,

≈ and ↓: Method with BU approach performs statistically significantly better, equal and worse.

g10 (x) = VFD (Velocity of front door at V-Pillar) ≤ 15.7 mm/ms All the constraints are modeled using response surface methodology as the surrogate models [26]. All the explicit models of constraints are solved based on the thicknesses of floor side inner, and its bounds are updated in order to handle all the variables. The maximum number of FE for this problem is set as 8000 and population-based algorithms were applied to this problem, with an initial population equal to 40. IP, SQP, DS, SA, PSO, GA, and DE algorithms were applied to this problem with and without using the proposed approach, and the statistical results are shown in Table 12. The statistical results of IP become worse after using the proposed algorithm, and only its median is improved. In Table 12, it can be seen that all of the SQP and SA results are improved. For the DS algorithm, the statistical results are enhanced, except for the best solution. Although the mean results of the PSO and GA are improved, their best solutions are slightly worse after using the BU approach. For the DE algorithm, all statistical results are boosted after using the proposed approach. Based on the Wilcoxon rank-sum test, the statistical results of SA and DE are significantly improved after using the BU approach. The simulation results for the population-based optimization algorithms are shown in Fig. 16. Fig. 17(a) to (c) show the convergence histories of the PSO, GA, and DE algorithms, respectively. From the results, it can be seen that although the BU has a slight improvement over the PSO and DE convergence histories, it improves the convergence histories of the GA significantly, which leads to better solutions. The percentages of constraint violations histories of this problem for the PSO, GA, and DE algorithms are shown in Fig. 17(d) to (f), respectively. From the GA and DE results, it can be seen that using the BU approach slightly decreases the constraint violations. However, the BU has a more significant effect on the constraint violations of the PSO algorithm and significantly decreases it during iterations. It should be noted that building the most effective surrogate models for the constraints is one of the biggest challenges in surrogate modeling [27]. By choosing a surrogate model that can satisfy the applicability condition, we can apply the proposed BU method to satisfy any constraints. This could be very useful for current real-world problems since not only are they time-consuming but also usually involve many constraints, such as mechanical and geometrical constraints. 3.4. Highly constraint and large-dimensional problem This problem is a benchmark problem for large-scale structural optimization [25]. The problem is also highly constrained, and its hit ratio is almost zero for a large number of segments. The schematic of this problem is presented in Fig. 18; the variables of this problem are the height (h) and the width of the cross-sectional area (b) of each segment. Fifty segments (N = 50) are considered here and, therefore, this problem has 100 dimensions.

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

29

Fig. 17. Convergence histories of (a) PSO, (b) GA and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA and (f) DE algorithms with and without BU approach for car side design problem.

Fig. 18. Schematic of the stepped beam design problem.

The objective is the minimization of the beam volume, which can be formulated as follows:

V =

N ∑

lbi h i

(21)

i=1

where l is the length of each segment. Each segment has one geometrical and mechanical constraint and, therefore, this case is subject to 100 constraints. Geometrical constraints: keep the aspect ratio between height and width of the cross-section of each segment less than 20, which can be simply formulated for the ith section as: hi −1≤0 20bi

(22)

30

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

Table 13 Statistical results of the different methods on the stepped beam design problem (Columns 4-7 indicate statistics on objective function values) Method

BU strategy

Ave No. FE

Best

Mean

Mediana

Worst

St. Dev.

IP

None BU1 BU2

2924 2536 2248

87 208.71 85 135.94 56 763.36

N/A N/A 79 855.44

106 861.17 93 077.31 80 114.62

N/A N/A 96 114.29

N/A N/A 8349.65

SQP

None BU1 BU2

3619 2842 8801

59 906.50 65 040.00 55 571.69

68 654.26 72 196.36 57 943.64

67 154.87 70 316.67 56 735.79

93 944.82 88 748.73 72 455.14

7472.95 5878.03 3399.89

DS

None BU1 BU2

100 000 100 000 100 000

58 296.72 57 444.91 55 569.66

59 876.03 58 684.22 55 604.62

59 858.76 58 643.48 55 603.64

61 293.94 60 100.69 55 648.97

688.28 695.88 16.80

SA

None BU1 BU2

100 000 100 000 100 000

81 013.99 82 634.74 74 009.79

N/A N/A 78 137.12

88 436.40 N/A 78 092.08 ↑

N/A N/A 81 870.91

N/A N/A 1830.10

PSO

None BU1 BU2

99 200 98 400 89 000

61 499.81 58 862.61 55 521.30

66 214.17 61 838.27 55 867.71

65 830.80 61 428.08 ↑ 55 965.11 ↑

72 106.51 74 806.33 57 501.82

2051.64 2303.30 403.97

GA

None BU1 BU2

99 000 100 000 100 000

69 908.97 62 276.63 60 142.12

74 955.93 64 226.44 62 421.96

74 815.98 64 295.30 ↑ 62 321.42 ↑

81 671.74 66 067.65 65 785.15

2224.53 832.13 1092.87

DE

None BU1 BU2

8900 8700 6700

N/A 95 941.05 83 219.58

N/A N/A 94 006.06

N/A 103 559.02 ↑ 93 726.05 ↑

N/A N/A 105 273.19

N/A N/A 3640.24

a ↑,

≈ and ↓: The approach performs statistically significantly better, equal and worse than the former approach.

Mechanical constraints: the bending stresses in each beam segment (σi ) should be less than the allowable stress (σa ), which is formulated as: σi −1≤0 (23) σa The BU approach for this problem is applied in two different strategies as: BU1: where only 50 geometrical constraints are handled; and BU2: All 100 geometrical and mechanical constraints are handled. In both strategies, only the lower bound of the width of each cross-sectional area is updated using the section constraint(s). The population is considered to be 100 and the maximum number of iterations is set to 1000 for all simulations. The optimization algorithms, with and without the BU strategies, are applied for this problem and the statistical results are presented in Table 13. First, it can be seen that using the proposed approach is effective for this highly constrained and large-dimensional problem. Handling the 50 constraints using BU1 significantly improves the statistical results of all optimization methods, except SA. After using the BU2, all the algorithm results are even more enhanced. For the SA algorithm, not only are the statistical results notably improved, but all SA simulations can find a feasible solution for this highly-constrained case study. From DE results, it can be seen that this algorithm cannot find the feasible area except when the proposed BU strategies are used for handling the constraints. Based on the Wilcoxon rank-sum test, the statistical results of metaheuristics are significantly improved after using BU approaches for implicitly handling more constraints. The convergence histories of the PSO, GA, and DE algorithms are presented in Fig. 19(a) to (e), respectively. In each figure, three different strategies compare algorithms without BU, with BU1, and with BU2. For the PSO convergence results, at first, it can be seen that both BU strategies find feasible solutions after a few FEs. However, PSO without BU approach may take many FE to even find the feasible search space. From Fig. 19(a), it is clear that after handling more constraints using the proposed approach, both convergence rate and final results are notably improved. From the GA convergence results (Fig. 19(b)), there are remarkable improvements after using BU strategies, in both convergence rates and final results. None of the DE runs without the BU approach can

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

31

Fig. 19. Convergence histories of (a) PSO, (b) GA and (c) DE algorithms and percentage of constraint violations histories of (d) PSO, (e) GA and (f) DE algorithms with and without BU1 and BU2 approaches for stepped cantilever beam design problem.

find the feasible region. When only geometrical constraints are considered in the BU, 16 runs still cannot find the feasible region. The DE algorithm seems not to work very well for this problem and may not be even able to find the feasible area after using DE with BU1 and without the BU approach. However, the only strategy that can find the feasible region in all runs is when DE is coupled with the BU and used to handle all the constraints. Fig. 19(d) to (f) show the percentages of constraint violations histories for the PSO, GA, and DE algorithms, respectively. Based on Fig. 19(d), the more constraints that are handled with the BU method, the fewer constraint violations. However, it is clear that without using the proposed approach, nearly all PSO solutions are stuck in the infeasible area. Using the BU1 for PSO, the infeasible solutions are more than 20% and, using the BU2, the violation percentage drops quickly. The GA can find a feasible area at its early generations and decreases the violation effectively. However, after using the CR approach, the number of violations significantly dropped and the violations are less than the GA without the CR. For the GA with the BU2 approach, the violation is almost zero after the initial stages, which shows that the CR approach can effectively cut the infeasible search space. From the violation histories in DE, it can be seen that the mean of violation histories is almost identical for the DE with BU1 and without the BU approach. This shows that only handling the geometrical constraints using the proposed approach are not very effective. However, it seems considering all constraints in the BU method (BU2), can significantly reduce the violations and effectively cut the infeasible region. This is because all the DE-BU2 could find a feasible area in the early stages. 4. Summary and conclusions In the current study, a new approach has been proposed for implicitly handling constraints in an optimization problem. The method is named boundary update (BU), as it has been designed to repair bounds of one or more selected variables in order to satisfy the constraints. In the BU method, boundaries of certain variables (termed repairing variables) have been restricted by directly using constraints. New solutions have been created within the updated boundaries in the search space. Therefore, this method cuts (or eliminates) the infeasible region in order to

32

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

speed up the search effort. Properly choosing repairing variables is an important step for the proposed BU method. However, we have found a way to determine repairing variable sets, which intend to select the repairing variable to handle the maximum number of constraints with the minimum number of repairing variables. Although the process involves pre-processing constraints to classify all variables into repairing and non-repairing variable sets, the ease of operation makes the approach a good compromise between suggested decoder-based approaches and generic constraint handling methods. Thus, the proposed approach is capable of handling a wide range of constrained problems with linear and nonlinear constraints, even grey or black-box type constraints. The approach has been first tested on two benchmark constrained numerical problems having linear and nonlinear constraints. Then, it has been tested on eight popular constrained engineering optimization problems. Next, a case study of surrogate models was presented as a remedy to when the BU method cannot be applied directly to solve constraints. Finally, the performance of the BU method was evaluated on a highly constrained, large-dimensional problem in which even finding a feasible solution is challenging. The results have clearly showed the advantages of the proposed BU method on all the above-mentioned problems. This approach is an implicit constraint-handling approach, which should be coupled with an explicit constraint-handling approach. Therefore, a study of a suitable explicit approach to couple with the BU method should be a future direction of this research. Automating and optimizing the proposed approach are other future directions of this research. In some cases, only using the BU approach has led the algorithm to feasible solutions or to the global optimum when the constraint handling method alone was unable to find a satisfactory solution. Thus, despite the challenges, the results have amply shown that the proposed BU method can be applied as a first approach to quickly obtain a feasible solution. Acknowledgments This material is based in part upon work supported by the National Science Foundation (NSF) under Cooperative Agreement No. DBI-0939454. Any options, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of NSF. Funding body There is no Funding support for this research. Ethical statement The paper has been submitted with full responsibility, following due ethical procedure, and there is no duplicate publication, fraud, plagiarism. None of the authors of this paper has a financial or personal relationship with other people or organizations that could inappropriately influence or bias the content of the paper. Appendix A. Supplementary data Supplementary material related to this article can be found online at https://doi.org/10.1016/j.cma.2020.112917. References [1] D.C. Charmpis, N.D. Lagaros, M. Papadrakakis, Multi-database exploration of large design spaces in the framework of cascade evolutionary structural sizing optimization, Comput. Methods Appl. Mech. Engrg. 194 (30–33) (2005) 3315–3330. [2] S. Gholizadeh, E. Salajegheh, Optimal design of structures subjected to time history loading by swarm intelligence and an advanced metamodel, Comput. Methods Appl. Mech. Engrg. 198 (37–40) (2009) 2936–2949. [3] A.H. Gandomi, K. Deb, R.C. Averill, S. Rahnamayan, M.N. Omidvar, Using semi-independent variables to enhance optimization search, Expert Syst. Appl. 120 (2019) 279–297. [4] J. Elzinga, T.G. Moore, A central cutting plane algorithm for the convex programming problem, Math. Program. 8 (1) (1975) 134–145. [5] S. Koziel, Z. Michalewicz, Evolutionary algorithms, homomorphous mappings, and constrained parameter optimization, Evol. Comput. 7 (1) (1999) 19–44. [6] E. Mezura-Montes, C.A.C. Coello, Constraint-handling in nature-inspired numerical optimization: past, present and future, Swarm Evol. Comput. 1 (4) (2011) 173–194. [7] K. Deb, An efficient constraint handling method for genetic algorithms, Comput. Methods Appl. Mech. Engrg. 186 (2) (2000) 311–338. [8] G.G. Dimopoulos, Mixed-variable engineering optimization based on evolutionary and social metaphors, Comput. Methods Appl. Mech. Engrg. 196 (4–6) (2007) 803–817.

A.H. Gandomi and K. Deb / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112917

33

[9] J. Rubio-Loyola, C. Aguilar-Fuster, G. Toscano-Pulido, R. Mijumbi, J. Serrat-Fernández, Enhancing metaheuristic-based online embedding in network virtualization environments, IEEE Trans. Netw. Serv. Manag. 15 (1) (2017) 200–216. [10] R. Eberhart, J. Kennedy, A new optimizer using particle swarm theory, in: Proceedings of the Sixth International Symposium on Micro Machine and Human Science, 1995. MHS’95, IEEE, 1995, pp. 39–43. [11] J.H. Holland, Adaptation in Artificial and Natural Systems, The University of Michigan Press, Ann Arbor, 1975. [12] R. Storn, K. Price, Differential evolution – A simple and efficient heuristic for global optimization over continuous spaces, J. Global Optim. 11 (4) (1997) 341–359. [13] A.H. Gandomi, X.S. Yang, S. Talatahari, S. Deb, Coupled eagle strategy and differential evolution for unconstrained and constrained global optimization, Comput. Math. Appl. 63 (1) (2012) 191–200. [14] Z. Michalewicz, Genetic algorithms, numerical optimization, and constraints, in: Proceedings of the Sixth International Conference on Genetic Algorithms, Vol. 195, Morgan Kaufmann, San Mateo, CA, 1995, pp. 151–158. [15] D.M. Himmelblau, Applied Nonlinear Programming, McGraw-Hill Companies, 1972. [16] S.S. Rao, Engineering Optimization: Theory and Practice, John Wiley & Sons, 2019. [17] E. Sandgren, Nonlinear integer and discrete programming in mechanical design optimization, J. Mech. Des. 112 (2) (1990) 223–229. [18] J.S. Arora, Introduction to Optimum Design, McGraw-Hill, New York, 1989. [19] J. Golinski, An adaptive optimization system applied to machine synthesis, Mech. Mach. Theory 8 (4) (1973) 419–436. [20] D. Kvalie, Optimization of Plane Elastic Grillages (Ph.D. thesis), Norges Teknisk Naturvitenskapelige Universitet, Norway, 1967. [21] X.S. Yang, A.H. Gandomi, Bat algorithm: a novel approach for global engineering optimization, Eng. Comput. 29 (5) (2012) 464–483. [22] A. Osyczka, Evolutionary algorithms for single and multicriteria design optimization, 2002. [23] W. Changsen, Analysis of Rolling Element Bearings, Mechanical Engineering Publications Ltd., London, 1991. [24] B.D. Youn, K.K. Choi, A new response surface methodology for reliability-based design optimization, Comput. Struct. 82 (2–3) (2004) 241–256. [25] G.N. Vanderplaats, Very Large Scale Optimization, National Aeronautics and Space Administration (NASA), Langley Research Center, 2002, NASA/CR-2002-211768. [26] L. Gu, R.J. Yang, C.H. Tho, M. Makowskit, O. Faruquet, Y.Y.L. Li, Optimisation and robustness for crashworthiness of side impact, Int. J. Veh. Des. 26 (4) (2001) 348–360. [27] K. Deb, R. Hussein, P.C. Roy, G. Toscano-Pulido, A taxonomy for metamodeling frameworks for evolutionary multiobjective optimization, IEEE Trans. Evol. Comput. 23 (1) (2019) 104–116.