Abstract
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 from mathematics and mechanical 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.
1. 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: where is compact and the functions and are continuous. Inequality constraints can be incorporated into by using the operator , .
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].where , , , is a real symmetric matrix, and is a continuous function such that , and .
The dual problem generated by such Lagrangians as in (1) is a (nondifferentiable) convex problem, which is usually solved by nonsmooth convex optimization techniques, such as subgradient methods and their extensions [2–8].
The IED method uses the family of Lagrangians (1). These Lagrangians induce dual variables , where is the classical multiplier associated with the linear term in the Lagrangian and 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 rulewhere is as in (1).
The iterate obtained in (2) can be used for computing the search direction which leads to the next dual iterate . 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–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 real-world 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.
2. The Problem in Study
The problems we are interested in are of the same type of the primal problem .
The augmented Lagrangian function associated with is defined as in (1). The corresponding concave dual function is defined as
The dual problem is given byand, since the dual function is concave, is equivalent to the convex problem
For convenience, we introduce the setwhose 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 as the dual function.
We recall now an important result that shows how, after solving problem (4), we can easily obtain a subgradient .
Theorem 1 (see [2], Lemma 2.1(a)). Given the dual function and a pair , take . Then
3. 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 problemwhere and are continuously differentiable.
Let be a regular point of problem . The Karush-Kuhn-Tucker (KKT) first order necessary optimality conditions are expressed as follows: If is a local minimum of then there exists such thatwhere is a diagonal matrix with .
The point such that is called a “Primal Feasible Point” and a “Dual Feasible Point”. Given an initial feasible pair , 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 aswhere 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.
Calling , the following system can be written from (9):
The solution of the system (10)-(11) provides a descent direction as proved in [19]. However, can not be employed as a search direction, since it is not necessarily a feasible direction because, in the case when , it follows from (11) that .
To avoid this drawback, a new system with unknowns and is defined by adding a negative matrix , with , to the right side of (11),
Now, (13) is equivalent to , . Consequently, if , the latest equation is reduced to , which means that is a feasible direction.
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 matrixand setting , which gives = + . In the case where we have , .
Thus, must be well chosen if . In this case, imposing , where , implies + ≤ which leads to Thus, with chosen as described above, we have that 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].
4. 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 is a convex, not necessarily differentiable function.
Problem can be replaced with the equivalent constrained problemIt is assumed that one arbitrary subgradient can be computed at any point .
NFDA starts at a point 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:where is a vector function with given by .
Problem is not solved. Indeed, it might not have a solution. NFDA uses its structure, which is similar to the structure of the problem , 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 and is a predefined parameter, since is not always finite.
A new point with is computed. If then and the step is called serious, as shown in Figure 2. Otherwise, if then and the step is called null; see Figure 3.
In any case, a new supporting hyperplane at the point is added to form a new auxiliary problem like .
These are the main ideas from NFDA; we need to understand the IED algorithm. For more on NFDA, see [9, 10].
5. 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 problemwhich, in its turn, is equivalent to
One can notice that problems 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 and a primal iterate such that .
Now, having , according to Theorem 1, we can easily find a subgradient
This subgradient defines a supporting hyperplane to the at the point . By employing this supporting hyperplane, the following auxiliary problem is definedwhere is given by .
Solving two systems, similar to (14)-(15) and (16)-(17), obtained from , the directions and 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 must be chosen so that is a descent direction for , i.e., 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., .
In other words, IED finds a direction which, at each iteration, belongs to the interior of the cone according to Figure 4.
Burachik et al. [1] proved that if satisfies then .
From the current iterate , using a suitable stepsize , an auxiliary point is computed.
If where then IED stops.
Otherwise, i.e., if and then .
If then a DSG step is performed in order to obtain the next iterate .
A new supporting hyperplane at the point is computed, a new auxiliary problem is defined, and the algorithm goes this way until a primal iterate satisfying is found.
We now describe formally the algorithm.
5.1. The IED Algorithm
Step 0. Fix a sequence and a sequence of positive definite matrices. Fix , , , where is the optimal dual value, a symmetric matrix , and an augmenting function . Choose , , , and . Take such that . Set .
Step k. Let , , and such that . Compute .
Step k.1. Set . For , find , , , and by solving systems
Find , satisfying = < < , with , and the search direction .
Compute the stepsize .
Step k.2. Compute the auxiliary point . Find . If , stop and set and .
Step k.3.a (serious step). If the auxiliary point then set , , , , and go to Step .
Step k.3.b (null step). Do . Choose where and . Compute a new point in by performing a DSG step from :
Find . If , stop. Otherwise, compute and take such that Set ; , and go to Step .
Figure 5 shows the serious step (a) and the null step (b) of the IED algorithm.
(a) Serious step:
(b) Null step:
As proved in [1], for the convergence of the method, the matrix , , must have the following structure: where is positive definite. Also, there must exist positive numbers and such that Moreover, there must exist so that for all .
To finalize this section, we present the result that justifies the stopping criterion used by the method.
Theorem 2 (see [2], Theorem 2.2). Suppose that for some we have . Then is a solution of and is a solution of if and only if .
Theorem 2 shows that a minimizer of the Lagrangian not only is the solution of but is also used to stop the algorithm.
6. The Genetic Algorithm
Genetic algorithms (GAs) [20–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 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 , that is, children tend to be close to their parents.
To maintain the average property, the offspring values are calculated as where , , and are the parents and, and are the offspring. It can be observed that . 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 , a value that works well according to the original reference [23].
6.2. 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 , 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 , a value that works well according to the original reference [24].
7. 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, , , , , , , and . 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 have been replaced with their nonsmooth equivalent so that the problems have been converted into equality constrained optimization problems.
As the stopping criterion for IED we have used . 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 - - 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.
7.1. 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. GLR-P1-1: subject to GQR-P1-1: subject to PPR-P1-2: subject to PQR-T1-7: subject to SQR-P1-1: subject to LGR-P1-1: subject to QBR-T1-1: subject to PBR-T1-1: subject to PBR-T1-2: subject to QQR-P1-3: subject to QQR-T1-6: subject to PLR-T1-1: subject to PBR-T1-3: subject to QLR-T1-1: subject to PQR-T1-1: subject to QQR-T1-3: subject to GBR-T1-1: subject to PQR-T1-4: subject to QQR-T1-2: subject to QPR-T1-1: subject to
7.1.1. Comparison of Results of the First Set of Numerical Experiments
Table 1 provides the initial guesses and 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, and 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.
7.2. 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 of active coils of the spring, the winding diameter and the wire diameter . The volume of the coil to be minimized is written as [27] and is subject to the constraints where
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 (), the module of teeth (), the number of teeth on pinion (), the length of the shaft 1 between the bearings (), the length of the shaft 2 between the bearings (), the diameter of the shaft 1 (), and, finally, the diameter of the shaft 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 bysubject to
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 , and . The objective function to be minimized is the cost of the beam given as subject to The expressions for , , , and involve the design variables and are given by
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 (in constant steps of 0.0625) and . 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.
7.2.1. Comparison of Results of the Second Set of Numerical Experiments
Table 22 provides the initial guesses and 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.
8. 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.
Data Availability
Data are available at https://www.dropbox.com/sh/5nahzdoyu37ug7z/AABYgt71fC_UB1EsP38xD08fa?dl=0 or from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors thank CNPq (Grant 306186/2017-9), FAPEMIG, and FADEPE (Grant 01/17 STRICTO SENSU).