A Hybrid Epigraph Directions Method for Nonsmooth and Nonconvex Constrained Optimization via Generalized Augmented Lagrangian Duality and a Genetic Algorithm

The Interior Epigraph Directions (IED) method for solving constrained nonsmooth and nonconvex optimization problem via Generalized Augmented Lagrangian Duality considers the dual problem induced by a Generalized Augmented Lagrangian Duality scheme and obtains the primal solution by generating a sequence of iterates in the interior of the epigraph of the dual function. In this approach, the value of the dual function at some point in the dual space is given by minimizing the Lagrangian. The first version of the IED method uses the Matlab routine fminsearch for this minimization. The second version uses NFDNA, a tailored algorithm for unconstrained, nonsmooth and nonconvex problems. However, the results obtained with fminsearch and NFDNA were not satisfactory. The current version of the IED method, presented in this work, employs a Genetic Algorithm, which is free of any strategy to handle the constraints, a difficult task when a metaheuristic, such as GA, is applied alone to solve constrained optimization problems. Two sets of constrained optimization problems frommathematics andmechanical engineering were solved and compared with literature. It is shown that the proposed hybrid algorithm is able to solve problems where fminsearch and NFDNA fail.


Introduction
We present a new version of the Interior Epigraph Directions Method published by Burachik et al. [1] for solving constrained nonsmooth and nonconvex optimization problems of the following kind: minimize  () over all  ∈  satisf ying  () = 0, () where  ⊂ R  is compact and the functions  : R  → R and  : R  → R  are continuous.Inequality constraints can be incorporated into () by using the operator  + = max{, 0},  ∈ R.
The IED method uses a Lagrangian duality technique, where each constraint function is appended to the objective function via a multiplier, or dual variable, to form the classical Lagrangian.In a duality scheme, instead of problem (), the problem solved is the one where the dual function is maximized.In this approach, the value of the dual function at a given point in the dual space is computed by minimizing the Lagrangian which is an unconstrained problem, and therefore unconstrained methods can be used.
It turns out that the optimal value of the dual problem associated with the classical Lagrangian is in general not the same as that of the primal problem (), if () is nonconvex.The difference between the two optimal values is called duality gap.
Zero duality gap and saddle-point properties without convexity assumptions are known to hold for the following family of Lagrangian function [2,Remark 2.1].
The IED method uses the family of Lagrangians (1).These Lagrangians induce dual variables (, ) ∈ R  × R + , where  is the classical multiplier associated with the linear term in the Lagrangian and  ≥ 0 is the penalty multiplier associated with the augmenting term.
A duality approach based on the Lagrangian  for solving () can be described as follows.Given a current dual iterate (  ,   ), a primal iterate is computed through the rule where  is as in (1).
The iterate   obtained in (2) can be used for computing the search direction which leads to the next dual iterate ( +1 ,  +1 ).Indeed,   is used to obtain a deflection of the subgradient direction which in turn is utilized to improve the dual values.This is the idea behind the methods studied in [2][3][4][5][6][7][8].These methods are referred to as deflected subgradient (DSG) methods.
The IED method [1] combines the Nonsmooth Feasible Directions Algorithm (NFDA) for solving unconstrained nonsmooth convex problems, firstly presented in [9] and further studied in [10], with a certain DSG method.
Roughly speaking, the IED method works as follows: starting at a point in the interior of the epigraph of the dual function, a deflected subgradient direction is used to define a linear approximation to the epigraph.An auxiliary point is obtained by solving the optimality conditions of a resulting linear problem.If the auxiliary point belongs to the interior of the epigraph, the IED step is declared serious and the auxiliary point becomes the next iterate, from which the process is repeated.If the auxiliary point does not belong to the interior of the epigraph, the IED step is declared null and a suitable version of the DSG step is applied from the original point.Thanks to the properties of the DSG method, this step stays in the interior of the epigraph.The new point takes the place of the original point, from which the process is repeated, until a serious step is performed.If an infinite number of null steps occurs, the method converges to a primal solution, as a consequence of the convergence properties of the DSG.
An important aspect of the method is the minimization of the Lagrangian function.Although the solver used for this task does not affect the convergence of the IED method, as proved in [1], solving this problem effectively is fundamental for the success of the methods based on Lagrangian duality.
The first version of the IED method [1] employs the Matlab routine fminsearch for the Lagrangian minimization.The second version of the IED method [15] uses an algorithm called Nonsmooth Feasible Directions Nonconvex Algorithm (NFDNA) [16,17] for solving nonsmooth and nonconvex unconstrained optimization problems which correspond to our problem (2).Although the previous versions of IED were able to solve many problems, they failed some test problems and did not find the desired solutions as one can observe through some numerical experiments presented in this paper and found in [1,15].We do not know why fminsearch and NFDNA fail.For fminsearch is a routine of Matlab, we think that it does not make sense trying to fix or change it.The NFDNA performance needs deeper investigation.In this work, we have used a Genetic Algorithm (GA) for minimizing the Lagrangian function.With this modification, the method gained robustness and was able to solve all the test problems analyzed here.The hybridization of the IED method and the Genetic Algorithm is the main purpose of this work and can be considered innovative.Also, it is important to notice that the GA is free of any strategy to handle the constraints, a difficult task when a metaheuristic, such as the GA, is applied alone to solve constrained optimization problems [18].
We have solved twenty problems and compared the performance of the current version of the IED method with the performance of the first version.We also solved four realworld mechanical engineering optimization problems with the current version of the IED.
This paper is organized in 8 sections.In the next section, we state the primal and the dual problems and recall how to find a subgradient of the dual function.In Section 3, we show how the search direction used by IED is obtained.In Section 4, we present NFDA, an algorithm for convex optimization problems that inspired the IED method.The IED method is then presented in Section 5.In Section 6, we present a brief description of the main characteristics of the Genetic Algorithm and its basic features adopted in this paper.Numerical results and comparisons are given in Section 7. Finally, conclusions and suggestions for future work are given in Section 8.

The Problem in Study
The problems we are interested in are of the same type of the primal problem ().
The augmented Lagrangian function  : R  × R  × R + → R associated with () is defined as in (1).The corresponding concave dual function  : The dual problem is given by max whose elements are obtained by the GA used in this work.
Figure 8 shows a flowchart that summarizes how such minimization is done.
For simplicity ( * ) will be our dual problem and we will refer to  fl − as the dual function.
We recall now an important result that shows how, after solving problem (4), we can easily obtain a subgradient  ∈ (, ).

The IED Search Direction
In this section we describe how the IED search direction is obtained.
The methodology employed here was firstly used by an algorithm called Feasible Directions Interior Point Algorithm (FDIPA) [19] for smooth nonlinear problems.
Let us consider the inequality constrained optimization problem minimize  () , where  : R  → R and  : R  → R  are continuously differentiable.
The point  such that () ≤ 0 is called a "Primal Feasible Point" and  ≥ 0 a "Dual Feasible Point".Given an initial feasible pair ( 0 ,  0 ), KKT points are found by solving iteratively the nonlinear system of equations ( 5) and ( 6) in (, ), in such a way that all the iterates are primal and dual feasible.
A Newton-like iteration to solve the nonlinear system of equations ( 5) and ( 6) can be stated as where (  ,   ) is the current point of the iteration , Λ is a diagonal matrix with Λ  ≡   , and ) is the Hessian of the Lagrangian function (, ) = () +   () or some quasi-Newton approximation.However,   must be symmetric and positive definite in order to ensure convergence.
The addition of a negative term −    produces a proportional deflexion of   into the feasible region, as one can see in Figure 1.The problem now is that   might not be a descent direction for the function .However, we can keep such property if   is properly chosen.Notice that the direction   can be obtained by solving two systems with the same matrix Thus, with   chosen as described above, we have that   =   1 +     2 is a feasible descent direction.The ideas just described are sufficient for us to understand how the search direction used by the IED method is computed.
For more on FDIPA, its features, and convergence, see [19].

NFDA for Convex Problems
The Nonsmooth Feasible Directions Algorithm (NFDA), firstly presented by Freire [9], has been devised for solving unconstrained nonsmooth convex problems of the following kind: where  : R  → R is a convex, not necessarily differentiable function.
It is assumed that one arbitrary subgradient  ∈ () can be computed at any point  ∈ R  .
NFDA starts at a point ( 0 ,  0 ) in the interior of the epigraph of the function .At a point (  , (  )), the method determines a supporting hyperplane to the epigraph of the function where   ∈ (  ) and defines an auxiliary constrained linear problem employing all the supporting hyperplanes computed so far, as follows: min  (, ) =  subject to   (, ) ≤ 0 where   = (  1 , . . .,    ) : R +1 → R  is a vector function with    : R +1 → R given by    (, ) = ℎ   () − .Problem ( 4.3 ) is not solved.Indeed, it might not have a solution.NFDA uses its structure, which is similar to the structure of the problem ( 3 ), to obtain a search direction by solving two systems and choosing the parameter  in the same way FDIPA does.
At the iterate (  ,   ) and having the search direction   at hand, a step   is obtained by the rule where  max = max{ |    ((  ,   ) +   ) ≤ 0} and  > 0 is a predefined parameter, since  max is not always finite.
In any case, a new supporting hyperplane at the point (  , (  )) is added to form a new auxiliary problem like ( 4.3 ).
These are the main ideas from NFDA; we need to understand the IED algorithm.For more on NFDA, see [9,10].

The IED Method
In this section we describe the Interior Epigraph Directions (IED) method for solving constrained, nonsmooth, and nonconvex optimization problems.
As we have seen in Section 2, the dual problem () is equivalent to the problem min  (, ) One can notice that problems ( 4.2 ) and ( * ) have the same structure.Therefore, we can apply to ( * ) the ideas used by NFDA.
However, NFDA can only minimize nonsmooth convex functions which have bounded level sets.The convex function , in general, does not enjoy this property.Therefore, the IED method uses directions which are suited to the function  we have.Indeed, the IED method takes advantage of both NFDA search direction and the special features of the dual function .
At the step  of the IED method, we have a point (V  ,   ) = ((  ,   ),   ) ∈ int(epi()) and a primal iterate   ∈ (V  ) such that (  ) ̸ = 0. Now, having   , according to Theorem 1, we can easily find a subgradient This subgradient defines a supporting hyperplane to the epi() at the point (V  , (V  )).By employing this supporting hyperplane, the following auxiliary problem is defined where Solving two systems, similar to ( 14)-( 15) and ( 16)-( 17), obtained from (AP  ), the directions   1 and   2 are computed.It is important to highlight that there is no duplication of the computational cost since the coefficient matrices of both systems are the same.
The search direction is thus defined by where the parameter   > 0 must be chosen so that   is a descent direction for , i.e., (  )  ∇(V  ,   ) < 0 and, at the same time, a local descent direction for the dual function .This requires   to be a descent direction for  as well, i.e., (  )  ∇(V  ,   ) < 0.
In other words, IED finds a direction  which, at each iteration, belongs to the interior of the cone according to Figure 4.
A new supporting hyperplane at the point (V +1 , (V +1 )) is computed, a new auxiliary problem is defined, and the algorithm goes this way until a primal iterate   satisfying (  ) = 0 is found.
We now describe formally the algorithm.
Figure 5 shows the serious step (a) and the null step (b) of the IED algorithm.
As proved in [1], for the convergence of the method, the matrix   ∈ R (+2)×(+2) ,  ∈ N, must have the following structure: where   ∈ R (+1)×(+1) is positive definite.Also, there must exist positive numbers  1 and  2 such that Moreover, there must exist  > 0 so that   >  for all  ∈ N.
To finalize this section, we present the result that justifies the stopping criterion used by the method.
Theorem 2 shows that a minimizer of the Lagrangian not only is the solution of () but is also used to stop the algorithm.

The Genetic Algorithm
Genetic algorithms (GAs) [20][21][22] were inspired by the Darwinian principles of natural selection and evolution.They are considered the most famous bioinspired metaheuristic and the basis of innumerous algorithms found in the related literature.GAs encode all the variables, using an adopted alphabet, as binary code, for example, corresponding to a candidate solution (defined as the phenotype) in a data structure (defined as the genotype).A population of genotypes evolves mimicking the evolutionary process found in the Nature.
The candidate solutions are submitted to the evaluation of their quality and they are selected by a stochastic process favoring better solutions.The genetic materials of these are recombined and mutated by means of genetic operators generating a new population.Basically, a GA presents the following five steps [21]: (i) a genetic coding of solutions to the problem; (ii) a procedure for creating an initial population of solutions; (iii) an evaluation function that returns the quality of each individual (fitness function); (iv) genetic operators manipulating genetic material of parents during the reproduction process giving rise to new individuals; and (v) parameters to be used in the algorithm during the recombination and mutation operations.
The GAs are free-derivative algorithms and massively used to solve unconstrained and constrained optimization problems.
The link between IED and Genetic Algorithm recognized through the notation   ∈ (V  ) can be seen in Figures 6 and  7.
The pseudocode depicted in Algorithm 1 describes the main steps of a generational Genetic Algorithm for minimizing the Lagrangian.
In this paper, a Genetic Algorithm with a real code is adopted and the Simulated Binary Crossover (SBX) [23] and  the polynomial mutation [24] have been set as the genetic operators.The next subsections describe these operators.The baseline real-coded generational GA uses a rank-based selection and elitism (the best element is always copied into the next generation avoiding the loss of good solutions).
6.1.SBX Crossover.SBX is a crossover based on the properties presented by the one-point crossover in binary-coded.These properties are as follows: (i) Average: the average of the parents is equal to the average of the children after the crossover operation; (ii) Spread factor : most of the crossover events result in a spread factor of  ≈ 1, that is, children tend to be close to their parents.
To maintain the average property, the offspring values are calculated as where  = (1/2)( 1 +  2 ),  2 >  1 ,  1 and  2 are the parents and, ch 1 and ch 2 are the offspring.It can be observed that ch = .The parameter  is a random variable obtained from the following probabilistic distribution function: where small values for  generate offspring which are more similar to their parents and, on the other hand, great values generate offspring which are quite different from their parents.In this work, we used  = 2, a value that works well according to the original reference [23].

Polynomial Mutation.
Polynomial mutation is an operator that uses polynomial probability distribution to perturb a solution in a parent's proximity.Given a parent solution  ∈ [, ], the polynomial mutation solution   is created by using a random number  ∈ [0, 1], for a given variable, as follows: where   and   can be calculated as follows: It is important to note that no value outside the range [, ] is created by the polynomial mutation.In this work, we used   = 100, a value that works well according to the original reference [24].

Numerical Experiments
In this section, we present the results obtained by the current version of the IED method which uses genetic algorithms for minimizing the Lagrangian.
Regarding the function  (see ( 1)), we have used seven norms in our first set of numerical experiments, namely, 2 ).Usually, those norms are used in order to check the efficiency of algorithm [1,2,25].Therefore, we have kept the norms to compare the results obtained by our algorithm.
With respect to the matrix , we have considered  =  which corresponds to the augmented Lagrangian case.
Inequality constraints () ≤ 0 have been replaced with their nonsmooth equivalent max{(), 0} = 0 so that the problems have been converted into equality constrained optimization problems.
As the stopping criterion for IED we have used ‖()‖ ≤ 10 −6 .A total of 10 independent runs have been executed by the IED + GA, extracting the best solution, the worst solution, the mean of the solutions, and the standard deviation (sd), for the first set of experiments.For the second set, we have executed 30 independent runs.
The algorithms have been implemented in Matlab -(2012) -in a microcomputer core (TM) i7 of 2.60 GHz with 8.00 GB of RAM.
This section is divided into two sets of test problems.

The First Set of Numerical
Experiments.This set of numerical experiments corresponds to the analysis of 20 constrained optimization problems extracted from the Hock and Schittkowski collection [26] and all of them refer to minimization problems.They are described below were the objective () and the constraints are presented.

Comparison of Results of the First Set of Numerical
Experiments.Table 1 provides the initial guesses  0 and  0 used by IED + fmin and IED + GA, the size of the population, the number of generations, and bound constraints used by the Genetic Algorithm.The table also shows the optimal values.For solving the first six problems, we have kept the initial guesses used by the IED in its previous version as described in [1].For the rest of the problems,  0 and  0 were randomly chosen.
One can notice from Tables 2-21 that IED + GA has found the solution for all the problems regardless of the norms used whereas IED + fmin failed some of them.
It is interesting to highlight that the standard deviation of all the problems is small.This shows robustness of the version of the IED method proposed in this work.
Tables 2-21 show the results of the problems obtained by both versions of the algorithm.For the IED + GA algorithm, we present the best and the worst results, the mean and the standard deviation.Cases where a solution has been not found are indicated by a line struck through the entry of the table.

The Second Set of Numerical Experiments.
In this section the algorithm proposed is applied to solve 4 real-world mechanical engineering optimization problems, largely used as a test-bed in the literature.
The Tension/Compression String Design.In this problem, the objective is to minimize the volume  of a coil spring, depicted in the Figure 8, under a constant tension/compression load.There are three design variables to be considered in this problem: The number  1 =  of active coils of the spring, the winding diameter  2 =  and the wire diameter  3 = .The volume of the coil to be minimized is written as [27] and is subject to the constraints Table 2: Comparison of the results obtained by both versions of the IED for problem GLR-P1-1.

𝜎(⋅)
LGR-P1-     The Speed Reducer Design.The objective of this problem is to minimize the weight  of the speed reducer [27] shown in the Figure 9.The design variables are the face width ( 1 = ), the module of teeth ( 2 = ), the number of teeth on pinion ( 3 = ), the length of the shaft 1 between the bearings ( 4 =  1 ), the length of the shaft 2 between the bearings ( 5 =  2 ), the diameter of the shaft 1 ( 6 =  1 ), and, finally, the diameter of the shaft 2 ( 7 =  2 ).The third variable is integer and all the others are continuous.The constraints include limitations on the bending and surface stress of the gear teeth, transverse deflections of the shafts 1 and 2 generated by the transmitted force, and, finally, the stress in the shafts 1 and 2. The weight of the speed reducer, to be minimized, is given by subject to (78) The Welded Beam.This problem corresponds to the design of the welded beam [27] depicted in the Figure 10.The design variables are {ℎ, , , }, with bounds 0.125 ≤ ℎ ≤ 10, and 0.1 ≤ , ,  ≤ 10.The objective function to be minimized is the cost of the beam given as  (ℎ, , , ) = 1.10471ℎ 2  + 0.04811 (14.0 + ) subject to The expressions for , ,   , and  involve the design variables and are given by  = √  ) . (81) The Pressure Vessel.This problem corresponds to the weight minimization of a cylindrical pressure vessel with two spherical heads [27] as shown in Figure 11.The objective function involves four variables: the thickness of the pressure vessel (  ), the thickness of the head ( ℎ ), the inner radius of the vessel (), and the length of the cylindrical component ().Since there are two discrete variables (  and  ℎ ) and two continuous variables ( and ), one has a nonlinearly constrained mixed discrete-continuous optimization problem.The bounds of the design variables are 0,0625 ≤   ,  ℎ ≤ 5 (in constant steps of 0.0625) and 10 ≤ ,  ≤ 200.The design variables are given in inches and the weight is written as to be minimized subject to the constraints The first two constraints establish a lower bound to the ratios   / and  ℎ /, respectively.The third constraint corresponds to a lower bound for the volume of the vessel and the last one to an upper bound to length of the cylindrical component.

Comparison of Results of the Second Set of Numerical
Experiments.Table 22 provides the initial guesses  0 and  0 used by IED + GA, the size of the population, and the number of generations used by the Genetic Algorithm.Table 23 presents the best and the worst results, the mean, and the standard deviation.Tables 24-27 show the results found by the algorithm.

Conclusions
We have proposed a version of the IED method for nonsmooth and nonconvex optimization problems that employs a Genetic Algorithm for minimization the Lagrangian function.
We have solved twenty test problems for comparison between the first version and the current version to verify the effectiveness of the new approach proposed.
We also solved four classical engineering problems to show that the method can be applied to more complex problems.
We have evidence that shows that the more effectively we minimize the Lagrangian the better the IED method will be.Therefore, more powerful methods must be sought for this task.We also believe that the IED method can be bettered by establishing rules to update the matrices   and employing different types of augmented Lagrangians with different matrices  and different functions .
Finally, the hybridization presented in this work combines two algorithms, namely, the IED and a Genetic Algorithm.This approach frees the Genetic Algorithm of handling with constraints which is troublesome for these type of methods and, at the same time, provides IED with a good technique for minimizing the Lagrangian.

Figure 6 :
Figure 6: A flow chart for a step of the IED method with a Genetic Algorithms.

Figure 7 :Figure 8 :
Figure 7: A flow chart of the Genetic Algorithm used inside the IED.

Table 1 :
Parameters of the IED method.

Table 3 :
Comparison of the results obtained by both versions of the IED for problem GQR-P1-1.

Table 4 :
Comparison of the results obtained by both versions of the IED for problem PPR-P1-2.

Table 5 :
Comparison of the results obtained by both versions of the IED for problem PQR-T1-7.

Table 6 :
Comparison of the results obtained by both versions of the IED for problem SQR-P1-1.

Table 7 :
Comparison of the results obtained by both versions of the IED for problem LGR-P1-1.

Table 8 :
Comparison of the results obtained by both versions of the IED for problem QBR-T1-1.

Table 9 :
Comparison of the results obtained by both versions of the IED for problem PBR-T1-1.

Table 10 :
Comparison of the results obtained by both versions of the IED for problem PBR-T1-2.

Table 11 :
Comparison of the results obtained by both versions of the IED for problem QQR-P1-3.

Table 12 :
Comparison of the results obtained by both versions of the IED for problem QQR-T1-6.

Table 13 :
Comparison of the results obtained by both versions of the IED for problem PLR-T1-1.

Table 14 :
Comparison of the results obtained by both versions of the IED for problem PBR-T1-3.

Table 15 :
Comparison of the results obtained by both versions of the IED for problem QLR-T1-1.

Table 16 :
Comparison of the results obtained by both versions of the IED for problem PQR-T1-1.

Table 17 :
Comparison of the results obtained by both versions of the IED for problem QQR-T1-3.

Table 18 :
Comparison of the results obtained by both versions of the IED for problem GBR-T1-1.

Table 19 :
Comparison of the results obtained by both versions of the IED for problem PQR-T1-4.

Table 20 :
Comparison of the results obtained by both versions of the IED for problem QQR-T1-2.

Table 21 :
Comparison of the results obtained by both versions of the IED for problem QPR-T1-1.

Table 22 :
Parameters of the IED method.

Table 23 :
Statistical Analysis obtained by IED + GA in the engineering problems.

Table 24 :
Comparison for The Tension/Compression String problem.

Table 25 :
Comparison for the speed reducer problem.

Table 26 :
Comparison for the welded beam problem.

Table 27 :
Comparison for the pressure vessel problem.