A Hybrid Global Optimization Method Based on Genetic Algorithm and Shrinking Box

This paper proposes a hybrid genetic algorithm method for optimizing constrained black box functions utilizing shrinking box and exterior penalty function methods (SBPGA). The constraints of the problem were incorporated in the fitness function of the genetic algorithm through the penalty function. The hybrid method used the proposed Variance-based crossover (VBC) and Arithmetic-based mutation (ABM) operators; moreover, immigration operator was also used. The box constraints constituted a hyperrectangle that kept shrinking adaptively in the light of the revealed information from the genetic algorithm about the optimal solution. The performance of the proposed algorithm was assessed using 11 problems which are used as benchmark problems in constrained optimization literatures. ANOVA along with a success rate performance index were used to analyze the model. Based on the results, we believe that the proposed method is fairly robust and efficient global optimization method for Constrained Optimization Problems whether they are continuous or discrete.


Introduction
Global optimization for Black-box functions, BBF, is still demanding even with the great advances in computational power of the modern computers as it is used in many applications including finite element analysis and computational fluid dynamics.Metamodels-based design optimization MBDO has been used for over two decades for optimizing such functions as they are using approximate models to optimize the BBF (Abu Bakar, Ramlan, Muhammad (2007), Bazaraa, Sherali, Shetty (2006), Belegundu, Arora (1985), Mahdavi, Fesanghary, Damangir (2007), Mallipeddi, Suganthan (2010)).Fitting response surfaces has been also used extensively in the literature to carry out MBDO (Camp (1955), Carroll & Fiacco(1961), Picheny, Wagner, Ginsbourger (2012), Pierskalla (1968)., Ramadan, S. & Ramadan, K. (2012), Regis & Shoemaker (2005)) such that the response surfaces were used to visualize the relationship between the values of the variables and the objective function values.This information was then used to estimate the location of the optimal value.Based on the estimated optimal location, some points were proposed for further expensive evaluations to improve the estimation.Mode-pursuing sampling model was also used to optimize MBDO Reference (Qiang & Changzhi (2014), Ramadan (2013), Simpson, Peplinski, Koch,Allen. (2001)).This model systematically generates more sample points around the mode of the black-box function and in the same time covers the entire sample space.Constrained Optimization Problems, COP, were also solved in literatures (Cohoon, Hedge, Martin, Richards (1987), Goldberg (1989), Hedar, & Fukushima (2002), Kazemi, Wang, Rahnamayan, & Gupta (2010), Moslem, Wang, Shahrayar, Kamal (2011), Sandgren (1990), Scholau,Welch, Jones (1998), Sharif, Wang, El Mekkawy (2008)).Jones (2001) mentioned three general categorize for constraints handling methods; penalty-based methods, biasing-based methods, and multiobjective-based methods.Of the three categorize, penalty-based methods are the most used one in the literatures because it has a simple theory and straight forward implementation.The principle behind penalty-based methods is that any violation in the constraints must finally make the objective function value infeasible.This means that the constraints must be impeded in the objective function itself using an auxiliary function known as penalty function, hence, transforming the COP into unconstrained optimization problem UCOP.These penalty functions can be further categorized into interior and exterior penalty functions (Chelouah, Siarry (2005), Chen, Qiu, Jiao (2013)).The exterior penalty function is designed such that a sequence of infeasible points is generated to lead to the optimal solution for the COP (Kramer (2010)).This method utilizes the maximum function to incorporate the constraints into the objective function and thus makes the resulting UCOP a non smooth problem.The non smoothness in the resulting UCOP calls for using derivative-free-based optimization methods (Lam (2008)).However, the advantage of this method is that with a reasonably large penalty parameter, which is not required to approach infinity, it can yield the optimum solution of the COP from the UCOP.Recently, hybrid stochastic-deterministic global optimization methods were gaining attention to solve COP (Coello, MezuraMontes (2002), Dellino, Kleijnen, Meloni (2010), Durand, Alliot (1999), Forrester,& Keane (2009), Huidae, Francisco, & Seth (2008)).These methods combine the benefits of the deterministic and stochastic methods.The deterministic methods such as Newton's method, conjugate gradient method and coordinate search method are known for their local search efficiency.The stochastic methods such as Genetic algorithm, Tabu search, and simulated annealing are known for their global search efficiency.In this study, a free sampling algorithm and easyhybrid global optimization method is proposed to solve COP using the stochastic genetic algorithm method and the deterministic shrinking box method.
Genetic algorithm GA is rooted in Darwin's theory of evolution.The theory is based on natural selection process in which the strong individuals (those who fit well in their environment) survive and their offspring prosper while the weak individuals (those who do not fit well in their environment) die out and thus their offspring eventually vanish.GA is highly immune to be trapped in local optima as it has two operators to enhance the global search, namely, crossover and mutation operators (Coello, MezuraMontes (2002), Moghaddam, Wang, Yannou, Wu (2006), Wang, Cai, Zhou, Fan (2009), Wu, Chow (1995)).The purpose of using the shrinking box method is to narrow down the search space for the GA which will enhance its ability to find the optimal solution.The proposed hybrid algorithm consists of two loops an outer loop and an inner loop.The outer shrinking box loop utilizes the results obtained from the inner genetic algorithm loop to reduce the search space and then the outer loop feeds back the new reduced search space into the inner loop.This process continues until a stopping criterion is met.The rest of the paper will be arranged as follows: section 2 introduces the proposed method.Section 3 describes the elements of the proposed method.Section 4 conducts ANOVA for the proposed method.Section 5 shows the experimentation.And finally, section 6 presents the conclusions.)

The Proposed Method
In this study a primary COP of the form: is considered where f: ∅ → ℝ is the objective function, ∅ is a card(θ)-dimensional sample space, θ is the set of the random variables, andg(θ) is the set of constraints.
Let g(θ) be the complementary set of functions for g(θ) such that if g(θ) is satisfied, g(θ) is not satisfied and vice versa.We can define a penalty function φ(g(θ), PenVal) as follows: φ(g(θ), PenVal) = ∑ PenVal × max 0, g̀ (θ) (1) where PenVal is the penalty coefficient corresponding to variable i and it is sufficiently large value, NC is the total number of constraints.Then This is true because g̀ (θ) and g (θ) are two mutually exclusive events due to the complementary relation between them.
Equation ( 1) can be used to penalize the primary objective function f(θ) if one or more values of the random variablesθ did not satisfy one or more constraints in g(θ).Hence an auxiliary UCOP can be deduced from the primary COP as follows: It should be noticed here that the set of primary constraints g(θ) includes the box constraints and any other constraints like integer or binary constraints on the random variables θ.Hence the constraints θ ∈ ∅ are incorporated in g(θ).Therefore these constraints will not show in the auxiliary UCOP.Lam (2008) showed that the optimal solution and the optimal value of the UCOP are the same optimal solution and optimal value for the COP if the penalty values PenVal are infinite because Unfortunately, making PenVal very large may cause serious computational difficulties (Lam (2008)).Reference (Jones, Schonlau, Welch (1998)) showed that PenVal does not need to be infinite for the UCOP to capture the optimal solution and the optimal value of the COP, it only needs to be sufficiently large.
As mentioned earlier, since the maximum function is integrated in the auxiliary UCOP, the auxiliary UCOP is not smooth.Therefore, all the derivative-based optimization methods cannot be applied directly to solve the UCOP.This matter calls for using derivative-free-based methods to solve it such as Genetic Algorithm.
In this study, a hybrid GA with Variance-based crossover VBC and Arithmetic-based mutation ABM operators along with a fitness function of the form Equation (3) will be used in the context of shrinking box method.This is a hybrid method between the stochastic nature of the GA and the deterministic nature of the shrinking box method.GA is used because it is good in global searching while shrinking box method is used to adaptively update and narrow the search space.This combination of stochastic-deterministic nature of the proposed method is aimed to increase the global and local search capabilities of the proposed method.

Elements of the Proposed Method
In this section we will present the different elements of the proposed method.

Shrinking Box
Let θ ∈ ∅ be the box constraints of the COP where∅ is the set of upper and lower limits of the random variablesθ such that where the vectors l and u ∈ ℝ represent the upper and lower bounds for θ respectively and n = card(θ).
The Shrinking box operator will adaptively change l and u according to the best solution found (best chromosome) as follows: where l , l are the new lower limit and the old lower limit of the random variable i, u , u are the new upper limit and the old upper limit of the random variable i, andge is the gene value in the best chromosome corresponds to the random variable i.

Initial Population
The genes that constitute the initial population chromosomes will be generated randomly from∅.The number of chromosomes in the initial population is equal to the population size PopSize which is related to the total generation number TGN as follows: where γ is a parameter between 0 and 1 corresponding to population size proportion and it must be tuned.The initial population will be generated in two steps as follows: While j < PopSize +1 For i=1: n Step 1 Discretize the values between l and u into stp values, where stp is a discretization parameter that must be set through tuning.
Give each discretized value a sequential number starting from 1 and ending at stp.
Step 2 Randomly select a value between 1 and stpand trace it back to its original value between l and u .
Assign this value to ge .

End
Assign the vector ge to Chromosome j, i.e.Chr j.

Variance-Based Crossover
The set of chromosomes (population) at any generation can be represented as a matrix M in which its row numbers represent the chromosomes numbers and the column numbers represent the random variables (genes) numbers.This means that a general gene can be represented as ge where j stands for chromosome j and i stands for gene i, i.e. variable i. Hence Chrij stands for gene i in chromosome j.The genotype representation of the population is as follows: Figure 1.Genotype representation for the population The VBC depends on the variance for the variable values in the different chromosomes.The steps involved in VBC calculations are as follows: Step 1 Assign a zero-one score to the values of the variables as follows: Step 2 Calculate the dispersion in the Sc values for the each variable (gene) over the population using the variance formula as follows: where Sc the average of variable i values over the populationand is calculated as: Step 3 Find the indices for the (n/2) genes corresponding to the smallest (n/2) variances (LowVar) and assign the remaining indices to HighVar Step 4 Generate NOffSpr offspring from the best chromosome found in the generation (BGChr) where NOffSpr is the number of offspring generated in the generation.The generation of the offspring OffSpr is as follows: The number of offspring generated in any generation is changing with the generation number as follows: NOffSpr= IniOffSpr+GN (12) where GN is the generation number, IniOffSpr is initial number of offspring and it is a parameter that must be tuned.This dynamic change in NOffSpr is aimed to reduce the variability toward the end of the generations to allow for a delicate search while keeping the variability high at the beginning to explore the search space adequately.

Arithmetic-Based Mutation
The ABM consists of the following 3 steps: Step 1 Choose MuGe = floor(r × n) random numbers between 1 and n, where r is the mutation rate and it is taken to be 0.1.
Step 2 Choose κ = floor(r × PopSize) random numbers between 0 and 1 Step 3 Calculate the mutated offspring as follows:

Immigration
The aim of the immigration operator is to keep the population diverse (Ramadan (2013)) and thus mitigate the effect of premature convergence.The immigration operator can be performed as follows: generate (PopSize-NOffSpr-1) chromosomes as explained in section 2.2.2.The number of immigrants reduces as the generation number increases because the offspring number increases according to Equation (4).This dynamic change in immigrant numbers aims to increase the variability in the population at the beginning to increase the exploration power of the method and mitigate the premature convergence problem.In addition, it aims to decrease the variability toward the end of the generations.
At this point, the number of chromosomes constituting the existing population is PopSize-1.The final step in the population preparation is the step of including the best chromosome found in the previous generation BGChr in the existing generation.This can be done by demanding that the chromosome number PopSize be the same as BGChr.This will make the population size equals to PopSize.

Stopping Criteria
The stopping criterion for this proposed method is to stop the algorithm when the objective function enhancement in the outer loop is zero while the stopping criterion for the GA i.e., the inner loop is when the generation number reaches TGN generations.
The flow diagram representing the proposed method is given in Figure 2 as follows:

Analysis of Variance Study for the Proposed Method
In the proposed SBPGA, the parameters that need tuning is 3, namely: stp, γ, and IniOffSpr.ANOVA study was conducted on the proposed method to assess the importance of the 3 parameters on its performance.Three levels were taken for each parameter and the average of 100 replications for the fitness function value of the six-hump camel-back SC benchmark problem was recorded as the response variable.The expression of the SC is giving below: The ANOVA results are shown in Table 1.shows the source of variability and column 6 shows the p-value for the F statistic.The test is conducted at 5% significance level.The results table shows that the discretization step size stp and the population size proportion γ significantly affect the best objective function value found.These results are interesting as it says that the discretization step size and the population size proportion are the only two parameters among the three parameters that significantly affect the performance of the model.The table also shows that the interaction between the stp and γis significant.This result is important as it suggests that discretization step size and population size proportion must be tuned simultaneously as they both affect each other and affect the best value found.On the other hand, the results show that the initial size of offspring IniOffSpr and both of its interactions are not significant."

Experimentation
The proposed method has been tested using 11 benchmark problems.The formulations of the problems are described in this paper for convenience.Each problem was solved 100 times and the average of these replications was used to construct a 95% confidence interval on the mean solution along with other central tendency and dispersion measures.A comparative analysis is also shown that compares the results of the proposed method with the results obtained from literatures.Moreover, the success rate for the different problems was calculated according to [44] and used as a performance measure for the proposed method.The success rate has the form:

Success rate = ,
where andf * andf * * is the obtained best solution and the best known solution, respectively.The machine used to solve the problems has the following specifications: Manufacturer HP, Model HPE-500f, Processor AMD phenon (tn) IIX6 1045T processor 2.70GHz, RAM 8.0 GB, system 64-bit operating system.
The descriptions and the results of the problems are as follows: A simple QF function with two variables.Reference [28] reported the mathematical expression for this function is as follows: The rate of success for this problem was 100%.The 95% confidence interval on the mean of the best value was [0.981, 4.834]E-5while the best value found was 3.8964e-12 corresponding to the best solution of [-1.000, 1000].The variance of 9.66E-9 indicates that the dispersion in the values of the best values is very low and thus the method was robust in this problem.The average time needed for each replication was 3.47s which is reasonable.Hence the method was efficient in terms of computational time in this problem.Table 2 shows a comparison between the results obtained for this problem using SBPGA and the results obtained using Reference [28].[28,29,36] among other references used this problem, the mathematical expression for this function is as follows: The rate of success for this problem was 100%.The 95% confidence interval on the mean of the best value was [-1.0315 -1.0315] while the best value found was -1.0315corresponding to the best solution of [0.0899, -0.7127].
The variance of 2.3295E-08 indicates that the dispersion in the values of the best values is very low and thus the method was robust in this problem.The average time needed for each replication was 3.51s which is reasonable.
Hence the method was efficient in terms of computational time in this problem.Table 3 shows a comparison between the results obtained for this problem using SBPGA and the results obtained using Reference [28].
Table 3.Comparison between SBPGA and other literatures for SC problem Goldstein-Price function (GP) with two variables.Refs.[28,37] among other references used this problem, the mathematical expression for this problem as follows: The rate of success for this problem was 81%.The 95% confidence interval on the mean of the best value was [3.0144, 3.0383] while the best value found was 3.0264 corresponding to the best solution of [0, -1].The variance of 0.3541indicates that the dispersion in the values of the best values is somewhat high and thus the method was not robust in this problem.The average time needed for each replication was 3.65s which is reasonable.Hence the method was efficient in terms of computational time in this problem.The rate of success for this problem was 100%.The 95% confidence interval on the mean of the best value was [-3.2778, -3.2538] while the best value found was -3.3224 corresponding to the best solution of [0.20170.15000.47690.27530.31170.6573].The variance of 0.0037 indicates that the dispersion in the values of the best values is somewhat low and thus the method was reasonably robust in this problem.The average time needed for each replication was 47.13s which less than a minute and thus it is considered reasonable.Hence the method was fairly efficient in terms of computational time in this problem.Table 5 shows a comparison between the results obtained for this problem using SBPGA and the results obtained using Reference [28].
Table 5.Comparison between SBPGA and other literatures for HF problem Griewank function (GN) with two variables.Refs.[28,39] among other references used this problem,the mathematical expression for this function as follows: Reference [28] -3.322 x ∈ −100,100 The rate of success for this problem was 98%.The 95% confidence interval on the mean of the best value was [-1.531, 9.604]E-4 while the best value found was 0 corresponding to the best solution of [0, 0].The variance of 8.0687E-06 indicates that the dispersion in the values of the best values is low and thus the method was robust in this problem.The average time needed for each replication was 3.40s which is low.Hence the method was efficient in terms of computational time in this problem.Table 6 shows a comparison between the results obtained for this problem using SBPGA and the results obtained using Reference [28].
Table6.Comparison between SBPGA and other literatures for GN problem Design of a Pressure Vessel (PV) with 4 variables and 3 constraints.This problem was solved in many references including [28,29,35,40].The pressure vessel is shown in Figure 3   The objective is to minimize the total cost given by: min Z = 0.6224x x x + 1.7781x x + x (3.1661x + 19.84x ) Subject to ASME boiler and pressure code for wall thickness x and x constraints given by: and tank volume The radius and the length are positive integers, while the shell thickness and the spherical head thickness are positive real numbers.The ranges for the design variables are as follows: x ∈ 25, 150 ,x ∈ 0.625, 1.0 ,x ∈ 1.0, 1.375 , andx ∈ 25, 240 The rate of success for this problem was 100%.The 95% confidence interval on the mean of the best value was [ 7012.77014.4]while the best value found was 7007.3 corresponding to the best solution of [51.82290.62501.000284.5169].The variance of 1.9182e-05 indicates that the dispersion in the values of the best values is low and thus the method was robust in this problem.The average time needed for each replication was 62.9s which is a little bit higher than a minute.Hence the method was fairly efficient in terms of computational time in this problem.Table 7 shows a comparison between the results obtained for this problem using SBPGA  The objective is to minimize difference between the teeth ratios and the standard parameter as follows: min Z = 0.14428 − x x x x s.t x ∈ 12, 60 ∀i ∈ 1, 2, 3, 4 The rate of success for this problem was 100%.The 95% confidence interval on the mean of the best value was [5.5988.708]E-10while the best value found was 2.7009E-12 corresponding to the best solution of [19164943].
The variance of 6.2916e-19 indicates that the dispersion in the values of the best values is low and thus the method was robust in this problem.The average time needed for each replication was 6.3s which is reasonable low.Hence the method was reasonably efficient in terms of computational time in this problem.Table 8 shows a comparison between the results obtained for this problem using SBPGA and the results obtained using Reference [47].The design variables are; wire diameter x1, the mean coil diameter x2, and the number of active coils x3. Figure 5 shows a schematic diagram for the coil spring.The model is given as: Subject to: The rate of success for this problem was 100%.The 95% confidence interval on the mean of the best value was [0.01268 , 0.012716] while the best value found was 0.012626 corresponding to the best solution of [0.05170.357311.2284].The variance of 6.4127E-09 indicates that the dispersion in the values of the best values is low and thus the method was robust in this problem.The average time needed for each replication was 28.4s which is reasonable.Hence the method was fairly efficient in terms of computational time in this problem.Table 9 shows a comparison between the results obtained for this problem using SBPGA and the results obtained using Refs.[48,49].
The rate of success for this problem was 100%.The 95% confidence interval on the mean of the best value was [0.06230.4474]E-3while the best value found was 3.6732E-24 corresponding to the best solution of [3.00000.5000].The variance of 9.6504E-07 indicates that the dispersion in the values of the best values is low and thus the method was robust in this problem.The average time needed for each replication was 30.4s which is reasonable.Hence the method was fairly efficient in terms of computational time in this problem.Table 11 shows a comparison between the results obtained for this problem using SBPGA and the results obtained using Reference [44].
The rate of success for this problem was 79%.The 95% confidence interval on the mean of the best value was [0.00350.0972]while the best value found was 2.4068E-07 corresponding to the best solution of [-0.20310.24420.15150.11640.07780.0407-0.2179-0.1487-0.06320.1475]E-3.The variance of 0.0570 indicates that the dispersion in the values of the best values is not low and thus the method was not robust in this problem.The average time needed for each replication was 166.6s which is high.Hence the method was not efficient in terms of computational time in this problem.Table 12 shows a comparison between the results obtained for this problem using SBPGA and the results obtained using Reference [44].

Conclusions
This paper, presents a simple, fast and forward global optimization method for COP based on shrinking box and GA methods.The proposed method showed a good degree of robustness in 8 out of the 11 problems, hence the method is considered fairly robust.In addition, the proposed method showed a good degree of computational time efficiency as it was efficient in 9 out of the 11 problems.Furthermore, the proposed method was able to have 100% success rate in 7 out of the 11 problems.The success rates for the remaining 4 problems are as follows: the GN problem has 98% success rate which is very high.Problems GP and ZV have 81% and 79% success rates respectively which are reasonable.Unfortunately, Ay problem has a low success rate of 44%.This low success rate is a result of the dimensionality problem.
Even though there is no guarantee that the best values found by this method are feasible due to using the auxiliary UCOP instead of the primary COP, all the best solutions found in the 11 problems are feasible.This is due to using high penalty value in the auxiliary function the matter that helped in eliminating the infeasible solutions.Unfortunately the high penalty value reduced the computational efficiency of the proposed method.
Based on these results, we believe that the proposed method is fairly robust and efficient global optimization method for COP whether they are continuous or discrete even though it has a stochastic element in it.Moreover, the proposed method is a standalone method that does not need any other algorithms such as sampling algorithms.Like most of the evolutionary algorithms, this method suffers from the curse of dimensionality and premature convergence problems.The shrinking box strategy was used in the proposed method to mitigate the effect of these problems.

Figure 2 .
Figure 2. Flow diagram for SBPGA with the corresponding dimensions.The shell consists of two welded halves of rolled carbon steel ASME SA 203 grade B steel plates.There are 4 design variables: radius (x1), spherical head thickness (x2), Shell thickness (x3), and length (x4).

Figure 3 .
Figure 3.A pressure Vessel [28] obtained using Refs.[46,48].Table 7.Comparison between SBPGA and other literatures for PV problem f * * = 7006.8GearTrain Problem (GT) with 4 variables and integer variables constraints.This problem was solved in Reference[29].The objective of this design problem is to find the closest gear ratio to the standard parameter given by the designer in [41] which is 0.14428.There are 4 design variables corresponding to the numbers of teeth of the four gears in the setup as shown in Figure4.

Table 1 .
ANOVA results for the parameters of the proposed method

Table 2 .
Comparison between SBPGA and other literatures for QF problem

Table 8 .
Comparison between SBPGA and other literatures for GT problemWeight of a spring problem (WS) with 3 variables and 3 constraints.This problem was solved in many references including Refs.[35, 42].The objective of this design problem is to minimize the weight of the spring.

Table 9 .
Comparison between SBPGA and other literatures for WS problemThe rate of success for this problem was 44%.The 95% confidence interval on the mean of the best value was [0.00560.3646]while the best value found was 2.1620E-11corresponding to the best solution of[0.2306-0.02550.20740.11610.34160.31360.30940.01590.0161-0.3845]E-5.The variance of 0.8387 indicates that the dispersion in the values of the best values is extremely high and thus the method failed in this problem.The average time needed for each replication was 167.5s which is high.Hence the method failed also in terms of computational time in this problem.Table10shows a comparison between the results obtained for this problem

Table 10 .
Comparison between SBPGA and other literatures for Ay problem

Table 11 .
Comparison between SBPGA and other literatures for Be problem

Table 12 .
Comparison between SBPGA and other literatures for ZV problem f * * = 0