Abstract

Bilevel optimization, as the name reflects, deals with optimization at two interconnected hierarchical levels. The aim is to identify the optimum of an upper-level leader problem, subject to the optimality of a lower-level follower problem. Several problems from the domain of engineering, logistics, economics, and transportation have an inherent nested structure which requires them to be modeled as bilevel optimization problems. Increasing size and complexity of such problems has prompted active theoretical and practical interest in the design of efficient algorithms for bilevel optimization. Given the nested nature of bilevel problems, the computational effort (number of function evaluations) required to solve them is often quite high. In this article, we explore the use of a Memetic Algorithm (MA) to solve bilevel optimization problems. While MAs have been quite successful in solving single-level optimization problems, there have been relatively few studies exploring their potential for solving bilevel optimization problems. MAs essentially attempt to combine advantages of global and local search strategies to identify optimum solutions with low computational cost (function evaluations). The approach introduced in this article is a nested Bilevel Memetic Algorithm (BLMA). At both upper and lower levels, either a global or a local search method is used during different phases of the search. The performance of BLMA is presented on twenty-five standard test problems and two real-life applications. The results are compared with other established algorithms to demonstrate the efficacy of the proposed approach.

1  Introduction and Background

A generic optimization problem involves maximization or minimization of one or more objective(s), often subject to certain constraint(s). In the presence of only one objective function, the problem is referred to as a Single-Objective (SO) optimization problem, whereas if it has multiple conflicting objectives, it is termed as a Multi-Objective (MO) optimization problem. A bilevel optimization problem is an extended version of a traditional SO or MO problem, where the optimization needs to be performed at two levels—upper and lower. In practice, such problems reflect cases where the performance of the upper-level system or authority is realizable only if the lower-level objective is optimum. For example, in engineering design, a solution optimized for performance (upper level) must satisfy certain physical conditions such as stability or equilibrium (lower level). There are several practical applications in the field of engineering (Kirjner-Neto et al., 1998), logistics (Sun et al., 2008), economics  (Sinha, Malo, Frantsev, and Deb 2013b), and transportation (Migdalas, 1995) that have inherent nested structure that needs to be modeled as a bilevel optimization problem. Given the increasing complexity of the problems emerging from these domains, there has been significant interest to develop efficient algorithms for bilevel optimization. The complexity could result from various factors such as large number of system variables, multiple objectives, and highly nonlinear objective/constraint functions. These factors combined with the nested nature of the problem pose severe challenges to the optimization methods. Furthermore, some of these problems need to be solved in real time, which calls for the need for computationally efficient algorithms. In this article, we present a memetic algorithm to deal with this challenge, one that attempts to deliver competitive solutions with relatively fewer function evaluations compared to existing techniques. The following subsection presents the mathematical details of bilevel optimization and a review of existing approaches. The proposed framework is described in Section 2. The metrics used for comparison of performance between different algorithms are presented in Section 3. The numerical experiments are presented on 25 benchmark test problems in Section 4 and two practical applications in Section 5. The conclusions are discussed in Section 6.

1.1  Mathematical Formulation

In this article, we limit our discussion to single-objective optimization problems involving only continuous variables. A single-objective optimization problem can be mathematically defined as shown in Equation 1.
formula
1

Here, is the objective function, is the set of inequality constraints, and is the set of equality constraints. is a vector of design variables in the domain  (defined using upper and lower bounds of each variable), where .

A single-objective bilevel optimization problem is an extended version of the problem described above, in which the optimization has to be performed at two levels, upper and lower. Each level has its own variables, objectives, and constraints. The subscript will be henceforth used to denote attributes of the upper-level problem, whereas subscript will be used for the lower-level problem. For a given upper-level vector , the evaluation of upper-level function is valid/feasible only if the for the corresponding lower-level problem (with held constant) is optimum. This critical relation leads to the nested nature of the bilevel problem. Mathematically, a bilevel problem can be represented as shown in Equation 2.
formula
formula
2

In Equation 2, the upper-level objective function is and the lower-level objective function is . The vectors and denote the upper-level variables (in domain ) and lower-level variables (in domain ) respectively. and are the sets representing inequality and equality constraints for the upper-level problem. Similarly, there are inequality constraints in and equality constraints in for the lower-level problem.

For the lower-level problem, acts as a fixed parameter. Corresponding to each , there is a lower-level optimum (=), which is used to evaluate the upper-level objective function . The overall goal is to optimize the upper-level (“leader”) objective function subject to the optimality of the lower-level (“follower”) problem. Bilevel problems can be difficult to solve and, interestingly, if the lower-level optimization task is non-convex, then it cannot be guaranteed that an algorithm finds even a single solution which is feasible in the context of bilevel optimization.

A number of the methods reported in the field of bilevel optimization operate using a nested mode. By a nested approach, it is meant that to evaluate an upper-level objective for a given set of upper-level variables ; the lower-level problem is solved using a suitable optimizer. If the lower-level optimizer is able to obtain an optimal (feasible) solution of the lower-level problem, then the set {} is used to evaluate the upper-level objective and constraints. This forms an individual evaluation for an optimizer solving the upper-level problem. If the lower-level optimization doesn’t return a feasible solution to the lower-level problem, the upper-level solution is discarded. There have also been several efforts toward developing non-nested techniques (especially in classical literature), which effectively convert the bilevel problem to a single-level problem. A brief review of notable approaches in the literature are discussed in the following subsections.

1.2  Classical Approaches

The term classical approach in the current context refers to the methods that either use analytical/gradient information via calculus, or may rely on derivative-free schemes such as branch and bound. They are usually straightforward to implement and require fairly low computational resources. However, their application is limited as they require assumptions on the nature of the objective and/or constraint functions, such as continuity, differentiability, linear/quadratic form, etc. Additionally, except for some cases (e.g., linear and convex problems), they identify a local optimum as opposed to global search methods, and hence are not suitable for multimodal problems.

For linear bilevel optimization problems, a number of classical approaches have been suggested in the literature. Some of these include simplex method (Dempe, 1987; Önal, 1993), branch and bound (Bard and Falk, 1982), descent methods (Vicente et al., 1994), adaptive tabu search (Gendreau et al., 1996), and penalty-based approach (White and Anandalingam, 1993). The choice of the algorithm is primarily dictated by the characteristics of the problem, such as the form of the objective function (convex/nonconvex), number of decision variables at upper and lower levels, and the type of constraints. The method presented in Dempe (1987) used a simplex method with additional rules to handle slack variables. It operates on a “full description” of the feasible set, obtained using the concept of active constraints and subgradients. However, it needs an interior feasible starting point. In contrast, the branch and bound method proposed by Bard and Falk (1982) can be applied to solve nonconvex problems and nonlinear problems. However, its complexity grows exponentially with the number of variables, as it is an exhaustive search technique. Steepest descent direction by Savard and Gauvin (1994) and the derivative-based approach by Kolstad and Lasdon (1990) are computationally faster, but may often converge to a local optimum. Additionally, for nondifferentiable or nonsmooth functions, numerical gradients are not accurate and can often mislead the search. In Colson et al. (2005), an approximation of nonlinear bilevel mathematical program was used in conjunction with trust-region–based methods. A derivation of generalized conditions for a local optimum in the context of a bilevel programming problem appears in Falk and Liu (1995). Two algorithms were proposed to solve the problems under these assumptions, namely trust-region bundle method with the subgradient information of the lower-level problems, and hybrid algorithm combining bundle methods and quasi-Newton methods. The approach was tested on two nonlinear problems. A tabu-search–based algorithm was proposed by  Rajesh et al. (2003), and its performance was demonstrated on a set of nine linear/quadratic test problems.

1.3  Complementary Constraint Approaches

A number of classical techniques use complementary constraints to reformulate the original bilevel problem into a single-level problem. The Karush-Kuhn-Tucker (KKT) condition is one such analytical condition that can be used to ensure local optimality at a given point and thus has been extensively used as a complementary constraint.

For linear bilevel programming, the complementary constraint approach using KKT has been used with the branch and bound method (Hansen et al., 1992), penalty approach (Lv et al., 2007), and mixed integer programming (Saharidis and Ierapetritou, 2009). Problems with up to 250 upper-level variables, 150 linear constraints, and up to 150 lower-level variables were solved using the technique in  Hansen et al. (1992). For nonlinear problems, KKT was used with exact and inexact penalty function approaches in Marcotte and Zhu (1996). A combination of duality theory and KKT optimality conditions was used in  Jenabi et al. (2013) to convert a bilevel problem (concerning electric power networks) into mixed integer linear and nonlinear programming problems. In  Herskovits et al. (2000), a contact shape optimization problem was formulated as a bilevel programming problem. They formulated and used alternative optimality conditions (in lieu of KKT) as the complementary constraints for this problem. In Sahin and Ciric (1998), a simulated-annealing (SA)–based algorithm was proposed to solve bilevel optimization problems. Their key idea was to stochastically relax or “fuzzify” the constraints that represent the lower-level problem during the search. This relaxation is controlled by a parameter (), whereas another parameter () is used for the annealing at the upper level; thus resulting in a “dual-temperature” SA. The approach was applied to finding a safe process layout for a chemical plant. Carrión et al. (2009) formulated a medium-term decision-making problem by a power retailer as a bilevel programming problem. To solve the problem, they transformed the bilevel problem to a single-level mixed integer problem using KKT optimality conditions, and converted the nonlinearity in the resulting problem to linear equivalents using well-known integer algebra results. Ye (2006) extended well-known constraint qualifications (Abadie, Zungwill, Kuhn-Tucker) used for nonlinear programming to bilevel programming problems. Based on these constraint qualifications, they proposed KKT-like optimality conditions which could be used as constraints for the reformulated problem. Subsequently, in Ye and Zhu (2010), classical and value function approaches were combined to derive new necessary optimality conditions for bilevel programming problems.

1.4  Population-Based Approaches

In recent decades, population-based metaheuristic optimization methods have gained popularity for a number of reasons. They are global search methods and do not require specific conditions on the mathematical form of the objective functions and are inherently suited to solve multi-objective problems. Such algorithms can also be readily implemented in parallel, which is especially advantageous for cases where evaluation of objective/constraint functions is significantly expensive compared to the cost of the operators. On the other hand, it is important to highlight that they cannot guarantee optimality. In the bilevel model, if the true optimum of the lower level is not achieved, then the solution is considered infeasible at the upper level. This is one of the key limitations of using metaheuristic methods for lower-level optimization problems, but several studies in literature have demonstrated their ability to achieve near optimal solutions (for both single and bilevel optimization problems). Besides, in some cases analytical expressions or their properties involved in the objective/constraints are unknown, for example, when their evaluation is the result of a numerical analysis such as Finite Element Analysis (FEA), Computational Fluid Dynamics (CFD), or a physical experiment. In such cases the classical techniques may not be applicable, and therefore heuristic/metaheuristic methods may be the only alternatives. However, the number of function evaluations used for optimization can be one of the concerns while using metaheuristics. Population-based methods in particular tend to evaluate a large number of solutions during the course of the search. This may not be affordable if the function evaluation is computationally expensive, as it would result in significant runtime for overall optimization exercise. Parallel computing can compensate for this issue where possible, and often a trade-off has to be made to get the best out of the available resources. Some of the population-based methods used in the context of bilevel optimization are outlined below.

A genetic-algorithm–based method called GABBA was introduced in  Mathieu et al. (1994). The approach utilized a mutation strategy specific to bilevel problems, and demonstrated improvements in accuracy compared to grid search (Bard, 1983).  Wang et al. (2007) solved a linear bilevel problem using an adaptive genetic programming approach. The approach incorporated intensification of search in the feasible space through evolutionary operators ensuring generation of feasible chromosomes. At the same time, probabilities of crossover and mutation were modified adaptively during the search. Yet another approach for linear problems was suggested in Hejazi et al. (2002), wherein a genetic algorithm was used to solve a transformed problem using KKT conditions. Oduguwa and Roy (2002) proposed a Bilevel Genetic Algorithm (BiGA), with a “co-evolutionary operator” to preserve the interactive nature of the decision making at two levels. An extension of this method was proposed in Legillon et al. (2012), referred to as the Coevolutionary Bi-level method using Repeated Algorithms (CoBRA). CoBRA incrementally improves two different subpopulations, each corresponding to one level, periodically exchanging information with each other. It was able to deal with linear bilevel optimization problems involving complex large-sized search spaces. Wang and Du (2014) used a genetic algorithm with a self-adaptive penalty function to solve constrained bilevel optimization problems. The use of Particle Swarm Optimization (PSO) has also been reported to solve bilevel linear programming problems by Kuo and Huang (2009) and Halter and Mostaghim (2006).

More recent methods have focused on nonlinear bilevel optimization problems, which are less tractable due to relatively complex landscapes. Some of the notable contributions in this area are by Sinha, Malo, and Deb (2013) and Sinha et al. (2014b). They introduced a scalable set of bilevel test problems (SMD series) incorporating different types of challenges such as nonlinearity, multimodality, and noncooperation between upper and lower levels. Problems of the SMD series include eight unconstrained and four constrained single-objective problems. Results using basic Nested Bilevel Evolutionary Algorithm (NBLEA) were presented for SMD1–SMD12 in Sinha et al. (2014b). The algorithm offered good convergence for unconstrained problems, but its performance on constrained problems was relatively inferior. In subsequent works (Sinha et al., 2014a; Sinha, Malo, and Deb, 2013), they proposed an improved evolutionary algorithm named Bilevel Evolutionary Algorithm on Quadratic Approximation (BLEAQ). In BLEAQ, computational cost was reduced by estimating the approximate values of lower-level optimum variables in lieu of exhaustive lower-level search. Consequently the accuracy of the final solutions obtained depends on the ability of the quadratic model to capture the relation between the lower-level optimum and upper-level variables accurately. Angelo et al. (2013) proposed a nested Bilevel Differential Evolution (BIDE)–based approach for single-objective problems. Exhaustive lower-level search was performed for each upper-level evaluation, resulting in competitive results, but at the expense of high numbers of function evaluations.

For single-level optimization problems, often the global search methods (such as evolutionary algorithms) are used in conjunction with local search methods (such as gradient search) to search for optimum solution(s) efficiently. This hybrid approach is referred to as memetic algorithm (Moscato, 1989). For a review on memetic algorithms, the readers are referred to Moscato et al. (2004). The advantages of using memetic algorithms have been reported in a number of prior studies including Singh et al. (2011) for conventional (single-level) optimization problems. When a bilevel problem is solved, parts of the strategy resemble single-level optimization, and a good technique for single-level optimization thus should have some potential of improving the performance for bilevel problems. However, there have been very few reports so far on application of memetic algorithms for the solution of bilevel optimization problems. Among these works, the approach has been to use a global optimizer at the upper level and a local optimizer at the lower level. Zhu et al. (2006) introduced Differential Evolution (DE) to solve leaders problems, while interior point method was used for solving follower problems. Similarly, in Koh (2013), DE was used to solve an upper-level problem and a gradient search for the lower-level optimization. Although these schemes aid in reducing the number of function evaluations, there is a considerable likelihood of obtaining a local optimum at the lower level, which affects the upper-level search.

To summarize the observations from the existing literature, it is apparent that in the context of bilevel optimization:

  • Classical approaches are relatively fast, but are severely limited in terms of the functions they can handle. Complementary constraints can be used to convert the problem to single-level, but they also need assumptions on the underlying functions and define only necessary conditions for optimality, often converging to local optimum for nonconvex problems.

  • Evolutionary approaches do not require assumptions on the objective functions, but tend to use a prohibitive number of function evaluations. Although the optimum cannot be guaranteed at either level, the performance is competitive if a sufficiently large number of evaluations is allowed.

  • The performance and applicability of hybrid approaches on bilevel optimization problems has been scarcely studied despite their success in solving conventional (single-level) single-objective optimization problems.

The principal motivation of the research presented in this article stems from the limitations of the works described above. The aim is to be able to solve generic unconstrained or constrained, linear or nonlinear problems, while using as few evaluations as possible. We also treat the objective/constraint functions as a “black-box,”that is, information about function characteristics such as linearity, convexity, and derivatives are not used as an input to the algorithm. The algorithm operates on input and output values only and wherever required, the gradients are calculated numerically. The intent is to make the algorithm suitable for a wider range of problems instead of being applicable to problems with very specific properties.

Toward the above goal, a Bilevel Memetic Algorithm (BLMA) is proposed in this article. The key features of the algorithm are:

  • BLMA uses a nested strategy, which is a combination of global (Differential Evolution) and local (Interior Point Algorithm) search methods. Unlike existing hybrid strategies the search is not prescribed to be of global nature at the upper level and local at the lower level. Instead, both levels use global or local search during different phases of the search. This is done to reduce the possibility of getting stuck at the local optimum at the lower level, as well as to make the search at the upper level efficient.

  • Local search at the lower level is enhanced by selecting a promising starting solution instead of choosing it randomly. To do this, an archive of all evaluated solutions () is maintained during the search. For a lower-level search corresponding to a new , the solution closest to it in the archive is identified as . The stored value corresponding to this solution is then selected as the starting point for the local search (more details in the next section).

  • In order to further ensure that the solutions are near-optimal at lower level, a re-evaluation is done of the best upper-level solution of the population in each generation. Re-evaluation in this context entails running a global optimizer for a large number of generations for the best in the population. All such selectively re-evaluated solutions are stored in a separate archive, from which the final solution is selected at the end of the run.

An extensive set of 25 benchmark problems has been used to evaluate the performance of the proposed algorithm. The performance is compared with other nested strategies (with fixed local or global search at both levels), and established bilevel optimization algorithms in the literature. In addition, results on two real-life applications are also presented.

2  Proposed Approach

In this section, the details of the proposed algorithm are discussed, starting with the overall nested framework and then discussing the specific underlying search strategies.

2.1  Overall Framework

The overall approach of BLMA is presented in Algorithm 1. The global search method used in this study is Differential Evolution (DE) and the local search method adopted is the Interior Point Algorithm (IP), discussed in the next subsection. The search at the upper and the lower levels proceeds as follows:

  • Upper level: At the upper level, the algorithm begins by initializing a set of solutions in space, using uniform random distribution. This forms the initial population for the global search (DE) at the upper level. For evaluation of each upper-level individual , an optimizer is run at the lower level to obtain an optimum lower-level solution and corresponding . The resulting () is then used to evaluate the upper-level function , which is the objective to be minimized by the upper-level algorithm.After the initialization, DE is run at the upper level until a prescribed number of generations . Of these, for the first number of generations, DE is used at the lower level to identify , whereas thereafter, as the upper-level solution approaches a potentially promising region, the lower-level optimization is switched to local search (IP) for faster convergence. is calculated as the nearest integer to , where is a user-defined parameter.After generations, the search at the upper level also switches from global to local. For the lower-level optimization, IP continues to operate. At the end of the run, the best obtained and the corresponding , along with the variables (, ) are reported as the final solution delivered by the algorithm. The pseudocode of the upper-level optimization process (as part of overall BLMA) is shown in Algorithm 1.

    formula
  • Lower level: As mentioned above, DE is used at the lower level for optimization for a prescribed number of upper-level generations (). The DE at the lower level uses a population size of , which is evolved for generations. The DE operators for lower-level optimization are the same as those used for upper-level DE, but they evolve the lower-level variables keeping fixed.Once the upper-level generation counter () exceeds , the search at the lower level switches to local (IP). For the local search, a starting solution needs to be selected. If the starting solution is chosen randomly, there is a higher chance that the final solution obtained would be a local optimum instead of a global one for multimodal functions. In order to improve the likelihood of the search converging to a global optimum, the starting solution needs to be selected more carefully. To do this, an archive of all evaluated solutions ( and the corresponding optimum ) encountered during the search is maintained and continually updated. The starting point is then identified as follows. The Euclidean distance of current is calculated from all other upper-level solutions in the archive. The solution in the archive with the least Euclidean distance from current is identified. The corresponding is chosen as the starting point for the local search at the lower level, with a prospect that solutions close in space are likely to have close values. Thus the available information about the landscape is exploited in order to provide promising starting solutions to the IP. The pseudocode of the lower-level optimization is given in Algorithm 2.

    formula
  • Re-evaluation: After each generation at the upper level, the best ranked upper-level solution is re-evaluated in order to ensure that its corresponding lower-level problem is optimized. In order to do this, a longer run of global optimization is performed for this , which uses same population size as used for the lower level, but five times the generations used for the other solutions (e.g., if the number of lower-level DE generations is set to 20, then re-evaluation uses 100). If the resulting is better than the original solution, it implies that the existing was not the optimum. In such a case, the solution is replaced with the re-evaluated one; otherwise, the original solution is kept. In case in a given generation the top-ranked solution has already been re-evaluated previously, then the second solution is chosen for re-evaluation instead (and so on). The upper-level population is re-sorted before moving to the next generation. A separate archive () of such re-evaluated solutions is maintained, and the final solution reported is the best solution (based on the upper-level objective) available from this archive ().

Thus, it can be seen that for the first few upper-level generations (), BLMA operates with global search at both the upper level and lower level; then for , it operates using global search at upper level and local search at the lower level, and thereafter, it uses local search at both levels. In this way, it is intended to avoid pitfalls of the algorithms that operate solely on global or local search at both levels. In addition, choosing an already evaluated solution (with the lower level optimized) from the provides a potentially good starting point for the local search at the lower level, improving the overall results, whereas the re-evaluation helps in further ensuring that the lower level is near-optimal for the best solutions. A preliminary version of the algorithm with no re-evaluation was previously studied by the authors in Islam et al. (2015).

2.2  Underlying Search Algorithms

As mentioned before, the global search used in the study is based on DE, whereas the local search is done using an IP. These algorithms are briefly discussed below.

2.2.1  Global Search

For global search, a -constrained differential evolution (DE) algorithm, proposed in Takahama and Sakai (2006), is used. The initialization of the population is done using a uniform random distribution in the search space.

For evolution, mutation and crossover operators are used. For each population member , mutation is used to generate a trial individual. This is done by choosing three unique individuals  from the population and generating the trial individual as:
formula
3
where is a parameter known as the scaling factor. This trial individual then undergoes crossover with the parent individual . In this algorithm, a binomial crossover with DE/1/rand/bin strategy is used for exploration, which has shown competitive results in the literature (Gong et al., 2011). A child individual is created using binomial crossover between the parent individual and the trial individual , by using the following operation for each variable:
formula
4

Here is a user-defined parameter known as crossover probability and is the number of variables. The term “randint[1,” denotes a random integer number in the range [1,] sampled using uniform distribution, whereas “rand [0,1]” denotes a random real number sampled using a uniform distribution in the range [0,1]. The child solution is accepted if it is better than the parent in terms of the objective value. Constraint handling is done based on level comparisons, as discussed in Takahama and Sakai (2006).

The global search using the DE is run for a prescribed number of generations. The best solution obtained is provided as an input to the local search phase, described next.

2.2.2  Local Search

In this study, Interior Point Algorithm (IP) based on Byrd et al. (1999) is chosen to complement the global search. In this approach, a barrier function is formulated and each of the subproblems are solved using two powerful techniques: Sequential Quadratic Programming (SQP) and Trust-Region (TR) methods. SQP is utilized since it can efficiently handle nonlinear constraints; whereas trust-region strategy allows for handling nonconvexity in the problem, in effect “globalizing” the SQP iteration. The strategies used in the algorithm are chosen to satisfy explicit conditions for global convergence stated in Byrd et al. (2000). The IP algorithm is used here for doing quick local improvements based on its ability to handle constraints and nonconvexity unlike other available local search methods. The implementation of this IP algorithm as available in the package within MATLAB 2013b is used in this study. Default parameters for the algorithm (as set in MATLAB) are used during optimization. When used for local search at the lower level, the obtained by IP for the given is used to evaluate the solution obtained at the upper level. When used for the upper level (during the last phase of BLMA), the IP outputs the best solution for the upper level, which is reported as the final solution delivered by BLMA (along with the corresponding , , ).

3  Comparison Metrics

In order to evaluate and compare the performance of the proposed approach, the following metrics are used. Due to the nested nature, the comparison of performance is not as straightforward for bilevel problems as it is in single-level problems. These challenges are also highlighted below.

  1. Objective values: For the standard (single-level) optimization problems, usually the performance could be judged by the statistics of the solutions obtained by each algorithm across multiple independent runs. For minimization, the algorithm which shows lower objective values statistically is considered to have better performance. However, this method of comparison does not hold for bilevel problems as an inaccurate optimum at the lower level may result in better objective value for the upper level than the true optimum. Thus, in order to unambiguously compare the performance, the availability of the true optimum is essential. This is the case for the Set 1 (SMD) problems considered above. The true optima for most problems in Set 2 are also known. For such problems, the performance is compared using an error metric, that is, for upper level , where is the best value obtained by the algorithm at the end of the run, and is the true optimum. Similarly, for the lower level, the . Lower values of and represent better accuracy. However, for the cases for which the true optimum is not explicitly available (BLTP1, BLTP3, BLTP12), the comparison is done based on the best-known objective values available in the literature.

  2. Computational effort: The next challenge lies in comparing the computational effort, as the upper and lower levels have a different number of function evaluations. This instigates a question about which of them should be compared or limited to a specified number for fair comparison. Due to the nested nature, and with different strategies of handling the function evaluations, it is almost always unattainable to terminate the algorithm at the same number of both upper and lower function evaluations. Hence in the numerical experiments that follow, some preliminary experiments are done to determine the algorithm settings required to obtain roughly similar number of function evaluations as obtained using the other algorithms. However, as expected, they often end up to be very different, depending on the search paths during different runs.It is to be noted that for real-life examples, this comparison should also take into account the relative computational complexity involved in evaluating a function at upper and lower levels. If upper-level evaluation is computationally intensive, then the total time taken will depend heavily on the upper-level function evaluations. The same goes for lower-level comparisons.

  3. Performance profile: Apart from reporting the median accuracy of the algorithms on individual problems, performance profile (Dolan and Moré, 2002; Barbosa et al., 2010) is also used for a visual comparison. The performance profile is a statistical tool for comparing the performance of multiple approaches against a given performance metric for a large set of problems. In a performance profile plot, the horizontal axis () signifies the performance values relative to the best performing algorithm for a problem. The vertical axis () signifies the cumulative distribution of the performance. For a given value of , it indicates what proportion of problems an algorithm was able to solve within a factor of the best performer. The algorithms could thus be compared based on a given level of performance . Additionally, overall performance could also be quantified using the area under the curve () of the profile. The larger value of is indicative of a better performance.

4  Numerical Experiments

In this section, we describe the numerical experiments conducted using BLMA on 25 benchmark problems. The performance of BLMA is reported on these problems and compared with four different baseline strategies and five other established bilevel optimization algorithms. This section first introduces the test problems and the algorithms used for comparisons, followed by the details of the experiments and discussion of results.

4.1  Test Problems

The test problems used in this study are divided into two sets: Set 1 and Set 2. Set 1 consists of the SMD problems proposed in  Sinha et al. (2014b). Set 2 is a collection of diverse problems from various other references in literature. The complete mathematical formulations are omitted for the sake of brevity, and could be found in the cited references.

4.1.1  Set 1: SMD Problems

SMD is a scalable set of test problems which cover a wide range of difficulties associated with bilevel optimization, including nonlinearity, multimodality, and conflict in upper-and lower-level objectives. The suite contains 12 problems of which SMD1–SMD8 are unconstrained, whereas the rest are constrained problems. At each level, the problem contains three components, and the overall objective is calculated as the summation of these three components. The generic form of the objective functions is as shown in Equation 5.
formula
5

Here , , and are the subcomponents of the upper-level objective function. Similarly, , , and are treated as subcomponents of the lower-level task. The dimension of the variables are set to , , , . For SMD6, . The value , , is set for all problems except for SMD6 for which , , , . These settings correspond to all problems having five variables, two at the upper level and three at the lower level for the studies presented. The true optimum values for each problem are reported in Table 1.

Table 1:
Bilevel test problems.
Set-1 (Sinha et al. (2014b))Set 2
ProblemProblem
NameNameSource
SMD1 (0,0) BLTP1 Aiyoshi and Shimizu (1984(0,200) 
   (Problem 2) 
SMD2 (0,0) BLTP2 Bard (1988) (Problem 1) (17,1) 
SMD3 (0,0) BLTP3 Angelo et al. (2013) (Problem 14) (1,0) 
SMD4 (0,0) BLTP4 Oduguwa and Roy (2002) (F1) (1000,1) 
SMD5 (0,0) BLTP5 Savard and Gauvin (1994(−1.4074,7.6172) 
SMD6 (0,0) BLTP6 Bard (1988) (Problem 3) (100,0) 
SMD7 (0,0) BLTP7 Anandalingam and White (1990(49,−17) 
SMD8 (0,0) BLTP8 Falk and Liu (1995(−1,0) 
SMD9 (0,0) BLTP9 Rajesh et al. (2003) (Problem 4) (85.0909,−50.1818) 
SMD10 (4,3) BLTP10 Rajesh et al. (2003) (Problem 2) (5,4) 
SMD11 (−1,1) BLTP11 Rajesh et al. (2003) (Problem 3) (9,0) 
SMD12 (3,4) BLTP12 Bard and Falk (1982) (Problem 2) (3.25,4) 
– BLTP13 Candler and Townsley (1982(29.2,−3.2) 
Set-1 (Sinha et al. (2014b))Set 2
ProblemProblem
NameNameSource
SMD1 (0,0) BLTP1 Aiyoshi and Shimizu (1984(0,200) 
   (Problem 2) 
SMD2 (0,0) BLTP2 Bard (1988) (Problem 1) (17,1) 
SMD3 (0,0) BLTP3 Angelo et al. (2013) (Problem 14) (1,0) 
SMD4 (0,0) BLTP4 Oduguwa and Roy (2002) (F1) (1000,1) 
SMD5 (0,0) BLTP5 Savard and Gauvin (1994(−1.4074,7.6172) 
SMD6 (0,0) BLTP6 Bard (1988) (Problem 3) (100,0) 
SMD7 (0,0) BLTP7 Anandalingam and White (1990(49,−17) 
SMD8 (0,0) BLTP8 Falk and Liu (1995(−1,0) 
SMD9 (0,0) BLTP9 Rajesh et al. (2003) (Problem 4) (85.0909,−50.1818) 
SMD10 (4,3) BLTP10 Rajesh et al. (2003) (Problem 2) (5,4) 
SMD11 (−1,1) BLTP11 Rajesh et al. (2003) (Problem 3) (9,0) 
SMD12 (3,4) BLTP12 Bard and Falk (1982) (Problem 2) (3.25,4) 
– BLTP13 Candler and Townsley (1982(29.2,−3.2) 

4.1.2  Set 2: BLTP Problems

The second set of 13 test problems considered here are collected from various references in the literature. The problems are well known in the bilevel optimization domain and have often been used to evaluate the efficacy of different bilevel optimization algorithms. Since a number of these problems appear in classical literature, they mostly involve linear and/or quadratic objective and constraint functions. For some of the problems (BLTP1, BLTP3, BLTP12), the true optimum values are not known, and hence the best considered in the literature are listed. For ease of reference, the problems have been labeled here as BLTP1–BLTP13 (BLTP= Bilevel Test Problems). The source reference for each problem and the optimum values are listed in Table 1.

4.2  Algorithms Used for Comparisons

As discussed in Section 2, the proposed BLMA operates using a combination of global and local search strategies to solve upper- and lower-level problems. Therefore, the performance of BLMA is first compared with a set of four basic nested algorithms, which are constructed by fixing the upper- and lower-level search methods to either global (DE) or local (IP). Thus, there are four combinations possible, which are briefly described below along with the parameters used.

  1. DE-DE: This nested approach uses the DE discussed previously for solving both upper- and lower-level problems. At both levels, the population size used is 50. The number of generations used is 28 for all problems, except the constrained problems SMD9–SMD12 for which 60 generations are used.1 The crossover probability and scaling factor are set to 0.9 and 0.7, respectively.

  2. DE-IP: This nested approach uses DE to optimize the upper level, and IP to optimize the lower level. The population size and generations used for DE are the same as in (1) above, whereas as for IP the standard parameters provided within MATLAB function are used. Further, for IP, multiple runs from 20 different starting solutions (identified using Latin Hypercube Sampling) are used and the best one is chosen as the lower-level optimum. The maximum number of evaluations for each IP run is set to 250.

  3. IP-IP: This nested approach uses IP to optimize both upper and lower levels. At the lower level, similar to (2) above, a maximum of 250 function evaluations is set with 20 restarts. At upper level, IP is run from random point(s) until a maximum of 1400 upper-level evaluations are reached for all problems except the constrained problems SMD9–SMD12 for which the limit is set as 3000 evaluations. The best solution obtained out of the upper-level IP is tagged as the overall best solution.

  4. IP-DE: In this case, IP is used to optimize the upper-level function, and DE for the lower-level function. The settings for IP and DE are the same as in (3) and (1) above, respectively.

Apart from the baseline strategies discussed above, the following three established algorithms in the field of bilevel optimization have been used to objectively assess the performance of BLMA:

  1. NBLEA (Sinha et al., 2014b): Nested Bilevel Evolutionary Algorithm (NBLEA) is an algorithm similar to DE-DE presented above, but uses a different Evolutionary Algorithm instead at both the levels. The implementation available for download at http://bilevel.org is used in this study, with population size 50 and default parameter settings.2

  2. BLEAQ (Sinha, Malo, and Deb, 2013): Bilevel Evolutionary Algorithm with Quadratic Approximation (BLEAQ) attempts to reduce the computational cost of a nested EA by estimating the values of the lower-level optimum variables instead of explicit search. Based on the truly evaluated solutions during the search, the lower-level optimum variables are approximated as a quadratic function of the upper-level variables. BLEAQ also shares some similarities with BLMA, for example, use of both local and global search methods and selection method of starting point for local search. The implementation available for download at http://bilevel.org is used in this study, with population size 50 and default parameter settings.

  3. BIDE (Angelo et al., 2013): This nested approach also uses a differential evolution-based algorithm to optimize the problems at both upper and lower levels. Exhaustive search is done at both levels to locate the optimum (in the study presented in Angelo et al., 2013), resulting in very high numbers of function evaluations. The code for the algorithm is not available; hence the results reported in Angelo et al. (2013) have been used for comparison.

  4. BiGA (Oduguwa and Roy, 2002): Bilevel Genetic Algorithm (BiGA) solves two optimization problems iteratively—one with upper level and a subset of lower variables, and one with lower variables only. Two separate populations are coevolved to solve each problem. The presented algorithm’s performance is compared with BiGA on a subset of BLTP problems, for which results are available from their article.

  5. NEA (Wang et al., 2005): This is a non-nested approach. In this method, the bilevel problem is transformed into a single-level problem with two objectives. A specific crossover and constraint handling strategy is also proposed to improve the performance. We compare the proposed algorithm for some of the BLTP problems studied in their article.

The experiments presented in the following subsections are divided into two parts. The first part deals with Set 1 (SMD problems), while the second deals with Set 2 (BLTP problems).

4.3  Numerical Experiments on Set 1 Problems (SMD)

For the proposed algorithm BLMA, a population size of 50 is used. The upper-level DE is run for a maximum of = 5 generations for all problems except SMD9–SMD12 for which 10 generations are used (since they are relatively harder problems to solve based on existing literature). For the lower level, DE uses a population of 50 and the number of generations . The value of is set to 0.8, which results in and 8 generations for the two different settings, respectively. A parametric study for different values of is presented later in Subsection 4.6. For IP, a maximum of 250 function evaluations is used for both unconstrained and constrained problems. The precision beyond 1e-6 is considered insignificant here, and hence any values below 1e-6 are simply treated as 1e-6.

For all the strategies studied here, 29 independent runs are conducted. We start by observing the performance of the four baseline strategies. The results obtained using them are shown in Tables 2 (median accuracy values) and 3 (median function evaluations). From Table 2, one can observe that DE-DE is able to reach close to the global optimum very consistently, but is unable to fully converge to the solution within given function evaluations. This behavior of population-based algorithms is not unexpected, as they depend on evolution operators even during final stages of the run to get to the optimum, which slows down convergence (as opposed to gradient-based or some other deterministic methods). At the same time, the numbers of function evaluations at the lower level are also high for this strategy owing to evolutionary search for each solution at the upper level. DE-IP shows similar behavior in terms of the function evaluations; in fact, the lower-level evaluations are higher for this case because of a high number of restarts at the lower level. However, the many restarts also increase the chance of getting to the true optimum of the lower-level problem, which is evident from the relatively low error values. One can observe that for most of the problems, at least the lower-level optimization obtains near-optimum results, with the exception of SMD4 and SMD9. Correspondingly, the upper-level accuracy is also good for most of the problems, albeit at high overall computational expense. It is seen that IP-IP shows fast convergence for most of the test problems. With similar computational expense, it is able to achieve better results in terms of accuracy for 8 out of 12 problems (SMD1-SMD3, SMD5-SMD8, SMD10) compared to the other three strategies. However, poor performance could be seen for some of the problems (e.g., SMD4, SMD9, SMD11-SMD12), which emphasizes the fact that IP can converge to a local optimum even with multiple starts. IP-DE has the worst performance among these four strategies, seemingly due to inaccurate lower-level optima which in turn misguide the upper-level search.

Table 2:
Median accuracy () using four baseline methods for SMD problems.
DE-DEDE-IPIP-IPIP-DE
Prb NameUL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()
SMD1 3.79E-04 2.89E-04 2.10E-05 9.00E-06 1.00E-06 1.00E-06 4.38E-02 4.64E-03 
SMD2 5.33E-03 6.46E-03 1.20E-05 4.00E-06 1.00E-06 1.00E-06 1.04E-01 4.47E-02 
SMD3 1.72E-03 1.02E-02 4.10E-05 1.20E-05 1.00E-06 1.00E-06 8.90E-02 1.10E-01 
SMD4 6.10E-02 9.37E-02 9.01E-01 9.53E-01 2.27E-01 9.51E-01 1.22E-02 1.17E-02 
SMD5 4.45E-02 5.22E-02 2.10E-05 9.00E-06 1.00E-06 1.00E-06 3.95E-02 1.80E-02 
SMD6 1.70E-05 2.00E-06 2.10E-05 9.00E-06 1.00E-06 1.00E-06 2.82E-04 1.28E-04 
SMD7 3.94E-03 4.80E-03 2.50E-05 1.00E-06 1.00E-06 1.00E-06 6.55E-02 6.49E-03 
SMD8 1.51E-02 3.02E-02 1.98E-03 3.03E-04 2.50E-05 1.00E-06 2.85E-01 3.50E-02 
SMD9 2.04E-04 2.54E-04 4.96E+01 5.41E+01 2.10E+01 1.88E+01 1.71E+01 1.28E+01 
SMD10 2.03E-01 2.50E-01 1.47E-02 4.69E-03 1.63E-03 3.90E-05 1.84E-01 7.76E-01 
SMD11 4.37E-01 5.12E-01 4.00E-06 4.00E-06 2.86E+00 2.76E+00 6.11E+00 5.11E+00 
SMD12 1.20E+00 2.62E+00 7.31E-02 1.00E-06 5.00E+00 1.62E+01 
DE-DEDE-IPIP-IPIP-DE
Prb NameUL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()
SMD1 3.79E-04 2.89E-04 2.10E-05 9.00E-06 1.00E-06 1.00E-06 4.38E-02 4.64E-03 
SMD2 5.33E-03 6.46E-03 1.20E-05 4.00E-06 1.00E-06 1.00E-06 1.04E-01 4.47E-02 
SMD3 1.72E-03 1.02E-02 4.10E-05 1.20E-05 1.00E-06 1.00E-06 8.90E-02 1.10E-01 
SMD4 6.10E-02 9.37E-02 9.01E-01 9.53E-01 2.27E-01 9.51E-01 1.22E-02 1.17E-02 
SMD5 4.45E-02 5.22E-02 2.10E-05 9.00E-06 1.00E-06 1.00E-06 3.95E-02 1.80E-02 
SMD6 1.70E-05 2.00E-06 2.10E-05 9.00E-06 1.00E-06 1.00E-06 2.82E-04 1.28E-04 
SMD7 3.94E-03 4.80E-03 2.50E-05 1.00E-06 1.00E-06 1.00E-06 6.55E-02 6.49E-03 
SMD8 1.51E-02 3.02E-02 1.98E-03 3.03E-04 2.50E-05 1.00E-06 2.85E-01 3.50E-02 
SMD9 2.04E-04 2.54E-04 4.96E+01 5.41E+01 2.10E+01 1.88E+01 1.71E+01 1.28E+01 
SMD10 2.03E-01 2.50E-01 1.47E-02 4.69E-03 1.63E-03 3.90E-05 1.84E-01 7.76E-01 
SMD11 4.37E-01 5.12E-01 4.00E-06 4.00E-06 2.86E+00 2.76E+00 6.11E+00 5.11E+00 
SMD12 1.20E+00 2.62E+00 7.31E-02 1.00E-06 5.00E+00 1.62E+01 
Table 3:
Median function evaluations using four baseline methods for SMD problems.
DE-DEDE-IPIP-IPIP-DE
Prb NameUL FELL FEUL FELL FEUL FELL FEUL FELL FE
SMD1 1.40E+03 1.96E+06 1.40E+03 1.69E+06 1.41E+03 1.61E+06 1.45E+03 2.01E+06 
SMD2 1.40E+03 1.96E+06 1.40E+03 1.69E+06 1.41E+03 1.56E+06 1.44E+03 2.00E+06 
SMD3 1.40E+03 1.96E+06 1.40E+03 2.67E+06 1.45E+03 2.77E+06 1.45E+03 2.01E+06 
SMD4 1.40E+03 1.96E+06 1.40E+03 2.54E+06 1.43E+03 2.64E+06 1.45E+03 2.00E+06 
SMD5 1.40E+03 1.96E+06 1.40E+03 2.85E+06 1.44E+03 3.79E+06 1.45E+03 2.01E+06 
SMD6 1.40E+03 1.96E+06 1.40E+03 1.30E+06 1.42E+03 1.25E+06 1.45E+03 2.01E+06 
SMD7 1.40E+03 1.96E+06 1.40E+03 1.71E+06 1.43E+03 1.70E+06 1.45E+03 2.00E+06 
SMD8 1.40E+03 1.96E+06 1.40E+03 3.69E+06 1.44E+03 4.34E+06 1.45E+03 2.00E+06 
SMD9 3.00E+03 9.00E+06 3.00E+03 2.46E+06 3.04E+03 3.24E+06 3.04E+03 8.99E+06 
SMD10 3.00E+03 9.00E+06 3.00E+03 4.96E+06 3.07E+03 4.96E+06 3.06E+03 9.10E+06 
SMD11 3.00E+03 9.00E+06 3.00E+03 9.26E+06 3.04E+03 4.31E+06 3.03E+03 8.95E+06 
SMD12 3.00E+03 9.00E+06 3.00E+03 3.92E+06 3.03E+03 3.92E+06 3.07E+03 9.11E+06 
DE-DEDE-IPIP-IPIP-DE
Prb NameUL FELL FEUL FELL FEUL FELL FEUL FELL FE
SMD1 1.40E+03 1.96E+06 1.40E+03 1.69E+06 1.41E+03 1.61E+06 1.45E+03 2.01E+06 
SMD2 1.40E+03 1.96E+06 1.40E+03 1.69E+06 1.41E+03 1.56E+06 1.44E+03 2.00E+06 
SMD3 1.40E+03 1.96E+06 1.40E+03 2.67E+06 1.45E+03 2.77E+06 1.45E+03 2.01E+06 
SMD4 1.40E+03 1.96E+06 1.40E+03 2.54E+06 1.43E+03 2.64E+06 1.45E+03 2.00E+06 
SMD5 1.40E+03 1.96E+06 1.40E+03 2.85E+06 1.44E+03 3.79E+06 1.45E+03 2.01E+06 
SMD6 1.40E+03 1.96E+06 1.40E+03 1.30E+06 1.42E+03 1.25E+06 1.45E+03 2.01E+06 
SMD7 1.40E+03 1.96E+06 1.40E+03 1.71E+06 1.43E+03 1.70E+06 1.45E+03 2.00E+06 
SMD8 1.40E+03 1.96E+06 1.40E+03 3.69E+06 1.44E+03 4.34E+06 1.45E+03 2.00E+06 
SMD9 3.00E+03 9.00E+06 3.00E+03 2.46E+06 3.04E+03 3.24E+06 3.04E+03 8.99E+06 
SMD10 3.00E+03 9.00E+06 3.00E+03 4.96E+06 3.07E+03 4.96E+06 3.06E+03 9.10E+06 
SMD11 3.00E+03 9.00E+06 3.00E+03 9.26E+06 3.04E+03 4.31E+06 3.03E+03 8.95E+06 
SMD12 3.00E+03 9.00E+06 3.00E+03 3.92E+06 3.03E+03 3.92E+06 3.07E+03 9.11E+06 

Next, the comparison between BLMA, NBLEA, BLEAQ and BIDE is presented in Tables 4 (in terms of accuracy) and 5 (in terms of function evaluations used3). The values shown are again median results from 29 runs. It is seen that for most of the problems, BLMA obtains very competitive results. For SMD2, SMD5, SMD7, and SMD10 it shows better/equal accuracy for both levels with fewer function evaluations compared to NBLEA, BLEAQ and BIDE as well as all four baseline strategies. The best (among different algorithms) median accuracy value obtained for each problem is shown in bold in Table 4.

Table 4:
Median accuracy () obtained using NBLEA, BLEAQ, BIDE, and BLMA on Set 1 (SMD).
BLMANBLEABLEAQBIDE
PrbUL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()
SMD1 1.00E-06 1.00E-06 5.03E-06(+) 2.37E-06(+) 1.00E-06(+) 1.00E-06(+) 3.41E-06 1.94E-06 
SMD2 1.00E-06 1.00E-06 3.17E-06(+) 3.43E-06(+) 5.44E-06(+) 5.50E-06(+) 1.29E-06 1.00E-06 
SMD3 1.00E-06 1.00E-06 1.37E-05(+) 8.95E-06(+) 7.55E-06(+) 5.50E-06(+) 4.10E-06 2.92E-06 
SMD4 1.00E-06 1.00E-06 9.29E-06(+) 6.19E-06(+) 1.00E-06(+) 1.86E-06(+) 2.30E-05 5.14E-05 
SMD5 1.00E-06 1.00E-06 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 1.58E-06 1.38E-06 
SMD6 1.00E-06 1.00E-06 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 3.47E-06 2.07E-06 
SMD7 1.00E-06 1.00E-06 3.69E-06(+) 1.00E-06(+) 5.81E-06(+) 9.23E-06(+) 
SMD8 4.55E-03 8.15E-04 4.74E-05(-) 1.54E-05(-) 2.21E-04(-) 5.53E-05(-) 
SMD9 2.00E-06 2.00E-06 4.91E-01(+) 7.55E-01(+) 4.22E-06(=) 1.16E-05(+) 
SMD10 4.00E-06 5.00E-06 2.47E-02(=) 4.08E-02(=) 1.02E-03(=) 8.55E-04(=) 
SMD11 2.89E-02 7.33E-03 5.99E-03(-) 8.78E-03(=) 1.28E-03(-) 2.08E-03(-) 
SMD12 4.00E-06 1.24E-02 3.82E-02(=) 7.81E-02(=) 4.56E-02(=) 2.00E-02(=) 
BLMANBLEABLEAQBIDE
PrbUL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()
SMD1 1.00E-06 1.00E-06 5.03E-06(+) 2.37E-06(+) 1.00E-06(+) 1.00E-06(+) 3.41E-06 1.94E-06 
SMD2 1.00E-06 1.00E-06 3.17E-06(+) 3.43E-06(+) 5.44E-06(+) 5.50E-06(+) 1.29E-06 1.00E-06 
SMD3 1.00E-06 1.00E-06 1.37E-05(+) 8.95E-06(+) 7.55E-06(+) 5.50E-06(+) 4.10E-06 2.92E-06 
SMD4 1.00E-06 1.00E-06 9.29E-06(+) 6.19E-06(+) 1.00E-06(+) 1.86E-06(+) 2.30E-05 5.14E-05 
SMD5 1.00E-06 1.00E-06 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 1.58E-06 1.38E-06 
SMD6 1.00E-06 1.00E-06 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 3.47E-06 2.07E-06 
SMD7 1.00E-06 1.00E-06 3.69E-06(+) 1.00E-06(+) 5.81E-06(+) 9.23E-06(+) 
SMD8 4.55E-03 8.15E-04 4.74E-05(-) 1.54E-05(-) 2.21E-04(-) 5.53E-05(-) 
SMD9 2.00E-06 2.00E-06 4.91E-01(+) 7.55E-01(+) 4.22E-06(=) 1.16E-05(+) 
SMD10 4.00E-06 5.00E-06 2.47E-02(=) 4.08E-02(=) 1.02E-03(=) 8.55E-04(=) 
SMD11 2.89E-02 7.33E-03 5.99E-03(-) 8.78E-03(=) 1.28E-03(-) 2.08E-03(-) 
SMD12 4.00E-06 1.24E-02 3.82E-02(=) 7.81E-02(=) 4.56E-02(=) 2.00E-02(=) 
Table 5:
Median function evaluations using NBLEA, BLEAQ, BIDE, and BLMA on Set 1 (SMD).
BLMANBLEABLEAQBIDE
Prb NameUL FELL FEUL FELL FEUL FELL FEUL FELL FE
SMD1 4.12E+02 3.05E+05 1.52E+03 9.52E+05 1.19E+03 2.37E+05 6.00E+03 1.80E+07 
SMD2 4.24E+02 3.01E+05 1.56E+03 9.63E+05 1.20E+03 4.08E+05 6.00E+03 1.80E+07 
SMD3 4.12E+02 3.09E+05 1.56E+03 1.04E+06 1.29E+03 3.02E+05 6.00E+03 1.80E+07 
SMD4 5.52E+02 3.29E+05 1.53E+03 8.33E+05 1.31E+03 3.07E+05 6.00E+03 1.80E+07 
SMD5 5.52E+02 3.28E+05 3.40E+03 2.22E+06 2.06E+03 8.42E+05 6.00E+03 1.80E+07 
SMD6 4.88E+02 3.26E+05 4.06E+03 1.11E+05 4.08E+03 1.98E+04 6.00E+03 1.80E+07 
SMD7 4.24E+02 3.01E+05 1.58E+03 9.56E+05 1.27E+03 3.45E+05 
SMD8 5.52E+02 3.25E+05 4.09E+03 2.85E+06 3.54E+03 1.48E+06 
SMD9 6.78E+02 5.30E+05 4.24E+03 3.08E+06 1.26E+03 4.35E+05 
SMD10 5.85E+02 5.13E+05 2.96E+03 2.93E+06 1.92E+03 5.25E+05 
SMD11 6.06E+02 5.12E+05 3.35E+03 3.10E+06 2.39E+03 2.26E+06 
SMD12 5.84E+02 5.10E+05 2.99E+03 3.25E+06 1.50E+03 4.15E+05 
BLMANBLEABLEAQBIDE
Prb NameUL FELL FEUL FELL FEUL FELL FEUL FELL FE
SMD1 4.12E+02 3.05E+05 1.52E+03 9.52E+05 1.19E+03 2.37E+05 6.00E+03 1.80E+07 
SMD2 4.24E+02 3.01E+05 1.56E+03 9.63E+05 1.20E+03 4.08E+05 6.00E+03 1.80E+07 
SMD3 4.12E+02 3.09E+05 1.56E+03 1.04E+06 1.29E+03 3.02E+05 6.00E+03 1.80E+07 
SMD4 5.52E+02 3.29E+05 1.53E+03 8.33E+05 1.31E+03 3.07E+05 6.00E+03 1.80E+07 
SMD5 5.52E+02 3.28E+05 3.40E+03 2.22E+06 2.06E+03 8.42E+05 6.00E+03 1.80E+07 
SMD6 4.88E+02 3.26E+05 4.06E+03 1.11E+05 4.08E+03 1.98E+04 6.00E+03 1.80E+07 
SMD7 4.24E+02 3.01E+05 1.58E+03 9.56E+05 1.27E+03 3.45E+05 
SMD8 5.52E+02 3.25E+05 4.09E+03 2.85E+06 3.54E+03 1.48E+06 
SMD9 6.78E+02 5.30E+05 4.24E+03 3.08E+06 1.26E+03 4.35E+05 
SMD10 5.85E+02 5.13E+05 2.96E+03 2.93E+06 1.92E+03 5.25E+05 
SMD11 6.06E+02 5.12E+05 3.35E+03 3.10E+06 2.39E+03 2.26E+06 
SMD12 5.84E+02 5.10E+05 2.99E+03 3.25E+06 1.50E+03 4.15E+05 

In comparison to the four baseline strategies, BLMA obtains better results for almost all cases. In comparison to BLEAQ, for six problems (SMD1, SMD3, SMD4, SMD6, SMD9, SMD12), the accuracy obtained by BLMA is better/comparable with marginally higher function evaluations, whereas BLMA is better for all others. In comparison with NBLEA, BLMA shows better/equal performance for all problems in both accuracy and function evaluations, except for SMD6, SMD8, and SMD11. For SMD6, BLMA obtains the same results but ends up using more function evaluations at the lower level compared to NBLEA, whereas for SMD8 it is marginally worse in accuracy at both levels but uses fewer function evaluations. Another notable observation from the results is the consistently good performance of BLMA for constrained problems (SMD9–SMD12). For the case of BIDE, the results for SMD1–SMD6 are taken from Angelo et al. (2013). The results using BLMA are better for these problems, even though it uses much fewer evaluations compared to BIDE. A statistical comparison of the results using NBLEA and BLEAQ against BLMA is done using the Wilcoxon signed rank test with 95% confidence level. The results are indicated in Table 4 using (=) if the values are statistically equivalent, (+) if they are higher, and (-) if they are lower. For graphical visualization of the results, the comparison between these four algorithms is also presented in Figure 1 (upper-and lower-level accuracies), and Figure 2 (upper-and lower-level function evaluations).

Figure 1:

Comparison of the median accuracy() obtained using NBLEA, BLEAQ, BIDE, and BLMA on Set 1 (SMD). The data points are joined with lines for easier visualization.

Figure 1:

Comparison of the median accuracy() obtained using NBLEA, BLEAQ, BIDE, and BLMA on Set 1 (SMD). The data points are joined with lines for easier visualization.

Figure 2:

Comparison of function evaluations used by NBLEA, BLEAQ, BIDE, and BLMA on Set 1 (SMD).

Figure 2:

Comparison of function evaluations used by NBLEA, BLEAQ, BIDE, and BLMA on Set 1 (SMD).

4.4  Numerical Experiments on Set 2 Problems (BLTP)

In this subsection, the numerical experiments on the BLTP problems (Set 2) are presented. The experiments are similar to those presented for the SMD problems, starting with baseline methods. The parameter settings for all algorithms are kept the same as those used for SMD experiments. The results shown for BIDE, NEA, and BiGA algorithm are taken from Angelo et al. (2013), Wang et al. (2005), and Oduguwa and Roy (2002), respectively.

The results obtained using the four baseline methods are shown in Tables 6 (median accuracy) and 7 (median function evaluations). From the results, it seems evident that all four strategies show better results on this set of problems compared to SMD problems. This can be attributed to the fact that all problems in this set are at most quadratic, in contrast with the SMD problems which incorporate relatively more nonlinear functions. DE-DE is able to get close to the optimum for all problems. DE-IP obtains reasonably good results, except for BLTP4. IP-IP performs well except for problems BLTP1 and BLTP4. The worst results are obtained by IP-DE, similar to what was observed for the SMD problems, because the upper-level local search is misguided by inaccurate lower-level optima.

Table 6:
Median accuracy () using four baseline methods for BLTP.
DE-DEDE-IPIP-IPIP-DE
PrbUL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()
BLTP1 1.00E-06 1.00E-06 1.00E-06 5.00E-06 1.12E+01 1.14E+02 1.28E+01 1.04E+02 
BLTP2 2.98E-02 1.00E-06 3.30E-05 8.70E-05 1.00E-06 1.00E-06 8.33E-03 9.45E-03 
BLTP3 1.00E-06 1.00E-06 1.00E-06 1.50E-05 4.00E-06 5.00E-06 1.00E-06 1.00E-06 
BLTP4 1.00E-06 1.00E-06 5.00E+02 1.00E-06 5.18E+02 1.00E-06 9.00E+02 1.00E-06 
BLTP5 5.60E-05 1.37E-04 3.50E-05 5.00E-05 7.00E-06 8.40E-05 1.10E-02 2.31E-02 
BLTP6 7.46E-04 1.00E-06 2.17E-03 1.00E-06 1.84E-01 8.40E-05 5.00E+01 4.29E+00 
BLTP7 4.28E-04 3.02E-04 5.57E-04 3.98E-04 6.00E-06 4.00E-06 2.62E-03 1.86E-03 
BLTP8 7.54E-03 6.40E-05 1.27E-04 1.00E-06 1.57E-03 1.00E-06 8.25E-03 1.00E-06 
BLTP9 1.00E-03 3.50E-04 2.15E-03 7.39E-04 9.00E-06 1.80E-05 3.68E-01 1.29E-01 
BLTP10 2.27E-03 2.34E-02 1.00E-06 3.60E-05 6.00E-06 8.60E-03 3.48E-04 3.09E-02 
BLTP11 7.44E-04 1.00E-06 1.00E-06 1.00E-06 2.83E-01 1.00E-06 4.24E-01 1.00E-06 
BLTP12 3.70E-04 2.96E-03 1.62E-04 4.31E-04 1.79E-04 9.20E-05 2.36E-01 3.91E-01 
BLTP13 2.61E-02 4.59E-03 1.23E-03 2.39E-04 3.67E-02 2.25E-03 9.39E-01 6.90E-02 
DE-DEDE-IPIP-IPIP-DE
PrbUL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()
BLTP1 1.00E-06 1.00E-06 1.00E-06 5.00E-06 1.12E+01 1.14E+02 1.28E+01 1.04E+02 
BLTP2 2.98E-02 1.00E-06 3.30E-05 8.70E-05 1.00E-06 1.00E-06 8.33E-03 9.45E-03 
BLTP3 1.00E-06 1.00E-06 1.00E-06 1.50E-05 4.00E-06 5.00E-06 1.00E-06 1.00E-06 
BLTP4 1.00E-06 1.00E-06 5.00E+02 1.00E-06 5.18E+02 1.00E-06 9.00E+02 1.00E-06 
BLTP5 5.60E-05 1.37E-04 3.50E-05 5.00E-05 7.00E-06 8.40E-05 1.10E-02 2.31E-02 
BLTP6 7.46E-04 1.00E-06 2.17E-03 1.00E-06 1.84E-01 8.40E-05 5.00E+01 4.29E+00 
BLTP7 4.28E-04 3.02E-04 5.57E-04 3.98E-04 6.00E-06 4.00E-06 2.62E-03 1.86E-03 
BLTP8 7.54E-03 6.40E-05 1.27E-04 1.00E-06 1.57E-03 1.00E-06 8.25E-03 1.00E-06 
BLTP9 1.00E-03 3.50E-04 2.15E-03 7.39E-04 9.00E-06 1.80E-05 3.68E-01 1.29E-01 
BLTP10 2.27E-03 2.34E-02 1.00E-06 3.60E-05 6.00E-06 8.60E-03 3.48E-04 3.09E-02 
BLTP11 7.44E-04 1.00E-06 1.00E-06 1.00E-06 2.83E-01 1.00E-06 4.24E-01 1.00E-06 
BLTP12 3.70E-04 2.96E-03 1.62E-04 4.31E-04 1.79E-04 9.20E-05 2.36E-01 3.91E-01 
BLTP13 2.61E-02 4.59E-03 1.23E-03 2.39E-04 3.67E-02 2.25E-03 9.39E-01 6.90E-02 
Table 7:
Median function evaluations using four baseline methods for BLTP.
DE-DEDE-IPIP-IPIP-DE
Prb NameUL FELL FEUL FELL FEUL FELL FEUL FELL FE
BLTP1 1.40E+03 1.96E+06 1.40E+03 6.94E+05 1.43E+03 7.38E+05 1.43E+03 1.97E+06 
BLTP2 1.40E+03 1.96E+06 1.40E+03 5.75E+05 1.43E+03 6.52E+05 1.41E+03 1.83E+06 
BLTP3 1.40E+03 1.96E+06 1.40E+03 2.91E+05 1.43E+03 2.93E+05 1.47E+03 2.04E+06 
BLTP4 1.40E+03 1.96E+06 1.40E+03 7.73E+05 1.40E+03 5.32E+05 1.41E+03 1.87E+06 
BLTP5 1.40E+03 1.96E+06 1.40E+03 1.13E+06 1.43E+03 1.10E+06 1.42E+03 1.92E+06 
BLTP6 1.40E+03 1.96E+06 1.40E+03 6.17E+05 1.43E+03 4.89E+05 1.42E+03 1.94E+06 
BLTP7 1.40E+03 1.96E+06 1.40E+03 6.01E+05 1.49E+03 6.07E+05 1.43E+03 1.96E+06 
BLTP8 1.40E+03 1.96E+06 1.40E+03 3.40E+06 1.45E+03 3.32E+06 1.45E+03 2.01E+06 
BLTP9 1.40E+03 1.96E+06 1.40E+03 1.32E+06 1.44E+03 2.05E+06 1.43E+03 1.96E+06 
BLTP10 1.40E+03 1.96E+06 1.40E+03 4.97E+05 1.43E+03 5.36E+05 1.42E+03 1.94E+06 
BLTP11 1.40E+03 1.96E+06 1.40E+03 2.11E+06 1.44E+03 2.13E+06 1.42E+03 1.96E+06 
BLTP12 1.40E+03 1.96E+06 1.40E+03 1.29E+06 1.46E+03 1.82E+06 1.43E+03 1.93E+06 
BLTP13 1.40E+03 1.96E+06 1.40E+03 1.79E+06 1.45E+03 2.33E+06 1.43E+03 1.97E+06 
DE-DEDE-IPIP-IPIP-DE
Prb NameUL FELL FEUL FELL FEUL FELL FEUL FELL FE
BLTP1 1.40E+03 1.96E+06 1.40E+03 6.94E+05 1.43E+03 7.38E+05 1.43E+03 1.97E+06 
BLTP2 1.40E+03 1.96E+06 1.40E+03 5.75E+05 1.43E+03 6.52E+05 1.41E+03 1.83E+06 
BLTP3 1.40E+03 1.96E+06 1.40E+03 2.91E+05 1.43E+03 2.93E+05 1.47E+03 2.04E+06 
BLTP4 1.40E+03 1.96E+06 1.40E+03 7.73E+05 1.40E+03 5.32E+05 1.41E+03 1.87E+06 
BLTP5 1.40E+03 1.96E+06 1.40E+03 1.13E+06 1.43E+03 1.10E+06 1.42E+03 1.92E+06 
BLTP6 1.40E+03 1.96E+06 1.40E+03 6.17E+05 1.43E+03 4.89E+05 1.42E+03 1.94E+06 
BLTP7 1.40E+03 1.96E+06 1.40E+03 6.01E+05 1.49E+03 6.07E+05 1.43E+03 1.96E+06 
BLTP8 1.40E+03 1.96E+06 1.40E+03 3.40E+06 1.45E+03 3.32E+06 1.45E+03 2.01E+06 
BLTP9 1.40E+03 1.96E+06 1.40E+03 1.32E+06 1.44E+03 2.05E+06 1.43E+03 1.96E+06 
BLTP10 1.40E+03 1.96E+06 1.40E+03 4.97E+05 1.43E+03 5.36E+05 1.42E+03 1.94E+06 
BLTP11 1.40E+03 1.96E+06 1.40E+03 2.11E+06 1.44E+03 2.13E+06 1.42E+03 1.96E+06 
BLTP12 1.40E+03 1.96E+06 1.40E+03 1.29E+06 1.46E+03 1.82E+06 1.43E+03 1.93E+06 
BLTP13 1.40E+03 1.96E+06 1.40E+03 1.79E+06 1.45E+03 2.33E+06 1.43E+03 1.97E+06 

Next, the median results obtained using BLMA, NBLEA, BLEAQ, and BIDE are shown in Table 8, and the corresponding function evaluations used by each algorithm are shown in Table 9. The best median objective value obtained for each problem is shown in bold in Table 8. It can be seen that BIDE gives the best results for most of the problems in this case, but at the cost of relatively much higher (nearly 100–1000 times) computational expense than other algorithms. BLMA is seen to converge to the global optimum (both and are 1e-6) for 6 problems out of 13. It also converges very close (1e-6 1e-3) to the optimum for the remaining ones, except for BLTP6 where it gets at the upper level. The results obtained by NBLEA and BLEAQ for BLTP6 are better than that of BLMA for this case, but with higher function evaluations at the upper level. For BLTP4, which seems to pose significant challenge to most algorithms studied (including baseline), the results obtained using BLMA are the best. A statistical comparison of the results using NBLEA and BLEAQ against BLMA is done using the Wilcoxon signed rank test with 95% confidence level, indicated as (=), (+), or (-) for equivalent, higher, or lower values.

Table 8:
Median accuracy () using NBLEA, BLEAQ, BIDE, NEA, BiGA, and BLMA on Set 2 (BLTP).
BLMANBLEABLEAQBIDENEABiGA
PrbUL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()
BLTP1 1.00E-06 1.00E-06 1.00E-06(+) 1.00E-06(+) 7.40E-06(+) 1.00E+02(+) 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
BLTP2 1.00E-06 1.00E-06 6.26E-04(+) 1.64E-03(+) 2.64E-03(+) 6.89E-03(+) 1.00E-06 1.00E-06 
BLTP3 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06 1.00E-06 1.11 6.07E-01 
BLTP4 1.00E-06 1.00E-06 3.78E+02(+) 1.00E-06(+) 3.87E+02(+) 1.00E-06(+) 6.34E+00 1.00E-06 2.10E-03 1.00E-06 1.00E-06 1.00E-06 
BLTP5 8.00E-06 8.80E-05 1.65E-03(+) 3.41E-03(+) 3.63E-03(+) 7.60E-03(+) 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
BLTP6 7.48E-01 1.39E-03 4.98E-03(-) 1.00E-06(-) 1.05E-02(-) 1.00E-06(-) 1.00E-06 1.00E-06 1.98E-01 1.55E-03 8.80E-01 1.00E-03 
BLTP7 1.20E-05 9.00E-06 1.00E-06(-) 1.00E-06(-) 1.00E-06(-) 1.00E-06(-) 1.00E-06 1.00E-06 7.05E-03 1.00E-06 
BLTP8 3.06E-03 1.00E-06 1.79E-05(-) 1.00E-06(+) 2.46E-05(-) 1.00E-06(+) 1.00E-06 1.00E-06 
BLTP9 5.00E-06 1.60E-05 3.75E-03(+) 1.30E-03(+) 2.47E-02(+) 8.63E-03(+) 6.11E-01 2.14E-01 
BLTP10 1.00E-06 1.02E-04 1.00E-06(+) 1.22E-03(+) 1.00E-06(+) 9.05E-04(+) 1.00E-06 1.00E-06 
BLTP11 1.00E-06 1.00E-06 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 1.00E-06 1.00E-06 
BLTP12 1.00E-06 1.00E-06 3.64E-01 3.02E+00 
BLTP13 7.65E-03 3.10E-04 1.03E-02(+) 2.00E-03(+) 4.30E-02(+) 8.32E-03(+) 1.00E-06 1.00E-06 1.00E-06 1.44E+00 
BLMANBLEABLEAQBIDENEABiGA
PrbUL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()
BLTP1 1.00E-06 1.00E-06 1.00E-06(+) 1.00E-06(+) 7.40E-06(+) 1.00E+02(+) 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
BLTP2 1.00E-06 1.00E-06 6.26E-04(+) 1.64E-03(+) 2.64E-03(+) 6.89E-03(+) 1.00E-06 1.00E-06 
BLTP3 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06 1.00E-06 1.11 6.07E-01 
BLTP4 1.00E-06 1.00E-06 3.78E+02(+) 1.00E-06(+) 3.87E+02(+) 1.00E-06(+) 6.34E+00 1.00E-06 2.10E-03 1.00E-06 1.00E-06 1.00E-06 
BLTP5 8.00E-06 8.80E-05 1.65E-03(+) 3.41E-03(+) 3.63E-03(+) 7.60E-03(+) 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
BLTP6 7.48E-01 1.39E-03 4.98E-03(-) 1.00E-06(-) 1.05E-02(-) 1.00E-06(-) 1.00E-06 1.00E-06 1.98E-01 1.55E-03 8.80E-01 1.00E-03 
BLTP7 1.20E-05 9.00E-06 1.00E-06(-) 1.00E-06(-) 1.00E-06(-) 1.00E-06(-) 1.00E-06 1.00E-06 7.05E-03 1.00E-06 
BLTP8 3.06E-03 1.00E-06 1.79E-05(-) 1.00E-06(+) 2.46E-05(-) 1.00E-06(+) 1.00E-06 1.00E-06 
BLTP9 5.00E-06 1.60E-05 3.75E-03(+) 1.30E-03(+) 2.47E-02(+) 8.63E-03(+) 6.11E-01 2.14E-01 
BLTP10 1.00E-06 1.02E-04 1.00E-06(+) 1.22E-03(+) 1.00E-06(+) 9.05E-04(+) 1.00E-06 1.00E-06 
BLTP11 1.00E-06 1.00E-06 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 1.00E-06(+) 1.00E-06 1.00E-06 
BLTP12 1.00E-06 1.00E-06 3.64E-01 3.02E+00 
BLTP13 7.65E-03 3.10E-04 1.03E-02(+) 2.00E-03(+) 4.30E-02(+) 8.32E-03(+) 1.00E-06 1.00E-06 1.00E-06 1.44E+00 
Table 9:
Median function evaluation using NBLEA, BLEAQ, BIDE, and BLMA on Set 2 (BLTP).
BLMANBLEABLEAQBIDENEABIGA
Prb NameUL FELL FEUL FELL FEUL FELL FEUL FELL FETotal FETotal FE
BLTP1 3.04E+02 2.81E+05 7.33E+02 1.29E+04 8.89E+02 4.62E+03 6.00E+03 1.80E+07 2.56E+05 
BLTP2 3.93E+02 2.84E+05 2.55E+03 2.41E+04 1.58E+03 4.76E+03 6.00E+03 1.80E+07 
BLTP3 3.03E+02 2.80E+05 2.96E+02 2.79E+03 3.05E+02 7.69E+02 6.00E+03 1.80E+07 5.01E+04 
BLTP4 3.39E+02 2.87E+05 9.00E+02 1.61E+04 9.93E+02 1.85E+03 6.00E+03 1.80E+07 2.85E+04 5.01E+04 
BLTP5 3.87E+02 2.85E+05 2.17E+03 3.86E+04 4.28E+02 3.92E+03 6.00E+03 1.80E+07 4.28E+03 
BLTP6 3.09E+02 2.81E+05 1.82E+03 1.71E+04 1.61E+03 3.24E+03 6.00E+03 1.80E+07 1.64E+05 5.01E+04 
BLTP7 3.38E+02 2.82E+05 3.61E+02 3.32E+03 3.60E+02 7.24E+02 6.00E+03 1.80E+07 4.34E+03 
BLTP8 3.70E+02 2.88E+05 1.70E+03 3.04E+04 1.37E+03 6.11E+03 6.00E+03 1.80E+07 
BLTP9 3.88E+02 2.92E+05 2.36E+03 2.23E+04 4.21E+02 1.92E+03 6.00E+03 1.80E+07 
BLTP10 4.24E+02 2.82E+05 2.09E+03 1.97E+04 1.08E+03 9.28E+02 6.00E+03 1.80E+07 
BLTP11 3.11E+02 2.83E+05 1.19E+03 1.12E+04 1.06E+03 9.02E+02 6.00E+03 1.80E+07 
BLTP12 3.18E+02 2.84E+05 6.00E+03 1.80E+07 
BLTP13 4.51E+02 2.91E+05 1.61E+03 4.36E+04 3.37E+02 4.49E+03 6.00E+03 1.80E+07 2.92E+05 
BLMANBLEABLEAQBIDENEABIGA
Prb NameUL FELL FEUL FELL FEUL FELL FEUL FELL FETotal FETotal FE
BLTP1 3.04E+02 2.81E+05 7.33E+02 1.29E+04 8.89E+02 4.62E+03 6.00E+03 1.80E+07 2.56E+05 
BLTP2 3.93E+02 2.84E+05 2.55E+03 2.41E+04 1.58E+03 4.76E+03 6.00E+03 1.80E+07 
BLTP3 3.03E+02 2.80E+05 2.96E+02 2.79E+03 3.05E+02 7.69E+02 6.00E+03 1.80E+07 5.01E+04 
BLTP4 3.39E+02 2.87E+05 9.00E+02 1.61E+04 9.93E+02 1.85E+03 6.00E+03 1.80E+07 2.85E+04 5.01E+04 
BLTP5 3.87E+02 2.85E+05 2.17E+03 3.86E+04 4.28E+02 3.92E+03 6.00E+03 1.80E+07 4.28E+03 
BLTP6 3.09E+02 2.81E+05 1.82E+03 1.71E+04 1.61E+03 3.24E+03 6.00E+03 1.80E+07 1.64E+05 5.01E+04 
BLTP7 3.38E+02 2.82E+05 3.61E+02 3.32E+03 3.60E+02 7.24E+02 6.00E+03 1.80E+07 4.34E+03 
BLTP8 3.70E+02 2.88E+05 1.70E+03 3.04E+04 1.37E+03 6.11E+03 6.00E+03 1.80E+07 
BLTP9 3.88E+02 2.92E+05 2.36E+03 2.23E+04 4.21E+02 1.92E+03 6.00E+03 1.80E+07 
BLTP10 4.24E+02 2.82E+05 2.09E+03 1.97E+04 1.08E+03 9.28E+02 6.00E+03 1.80E+07 
BLTP11 3.11E+02 2.83E+05 1.19E+03 1.12E+04 1.06E+03 9.02E+02 6.00E+03 1.80E+07 
BLTP12 3.18E+02 2.84E+05 6.00E+03 1.80E+07 
BLTP13 4.51E+02 2.91E+05 1.61E+03 4.36E+04 3.37E+02 4.49E+03 6.00E+03 1.80E+07 2.92E+05 

A notable observation from Table 9 is that the numbers of lower-level evaluations used by NBLEA are very small. For example, for BLTP3, the total number of upper-level function evaluations is 296, whereas the lower-level function evaluations total 2790, less than 10 per upper evaluation, even though the lower-level (default) population size is 50. This is because of a modification done in NBLEA code further to that reported earlier in Sinha et al. (2014b). The current version of the code (available at http://www.bilevel.org/) first checks if the problem is solvable using Quadratic Programming (QP). If so, instead of using the evolutionary algorithm, QP is used to solve the lower-level problem. For the BLTP problems discussed here, this criterion is often satisfied as the constraints are linear and objective functions are linear or quadratic. Thus, the resulting function evaluations are from QP which are much lower compared to an EA. On the other hand, such low numbers were not observed for SMD problems, indicating they could not be solved using QP. In the proposed BLMA, this knowledge about the test problems was not considered. Consequently, the lower-level function evaluations are large, but at the same time, it can be seen that BLMA is able to achieve close to global optimum solution for all problems. This includes BLTP4 where NBLEA and BLEAQ obtain poor upper-level solutions. For the case of both NBLEA and BLEAQ, no feasible solution was found for BLTP12 using the default algorithm settings, and hence the data are omitted. The comparisons with NEA (Wang et al., 2005) and BiGA (Oduguwa and Roy, 2002) are done for only a few problems for which the results were available. The median accuracy for these problems achieved using BLMA is also better than these algorithms. A graphical comparison of accuracy obtained and function evaluations used by the algorithms is given in Figures 3 and 4, respectively.

Figure 3:

Comparison of the median accuracy() obtained using BLMA, NBLEA, BLEAQ, and BIDE for Set 2 (BLTP). The data points are joined with lines for easier visualization.

Figure 3:

Comparison of the median accuracy() obtained using BLMA, NBLEA, BLEAQ, and BIDE for Set 2 (BLTP). The data points are joined with lines for easier visualization.

Figure 4:

Comparison of function evaluations used by NBLEA, BLEAQ, BIDE, and BLMA for Set 2 (BLTP).

Figure 4:

Comparison of function evaluations used by NBLEA, BLEAQ, BIDE, and BLMA for Set 2 (BLTP).

The interquartile range is also observed in order to check the robustness of the performance obtained by different algorithms. Interquartile range is the difference between the median values of the upper and lower quartile, representing a statistical measure of dispersion. The values are shown in Table 10. It is seen that the spread is low for most of the problems, representing that the algorithms converged in the vicinity of the global optimum consistently, with some exceptions (e.g., BLMA for SMD10, NBLEA for BLTP4, and BLEAQ for BLTP1).

Table 10:
Comparison of interquartile range among BLMA, NBLEA, and BLEAQ.
BLMANBLEABLEAQ
Prb NameULLLULLLULLL
SMD1 1.00E-06 1.00E-06 6.15E-06 4.29E-06 1.00E-06 1.00E-06 
SMD2 1.00E-06 1.00E-06 6.07E-06 1.26E-05 1.59E-05 2.50E-05 
SMD3 1.00E-06 1.00E-06 1.31E-05 1.10E-05 1.12E-05 7.60E-06 
SMD4 1.00E-06 1.00E-06 1.48E-05 9.20E-06 1.00E-06 4.37E-06 
SMD5 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
SMD6 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
SMD7 1.00E-06 1.00E-06 1.45E-05 1.69E-06 9.75E-02 2.43E+02 
SMD8 4.97E-03 1.25E-03 1.31E-04 4.01E-05 4.80E-04 1.19E-04 
SMD9 1.00E-06 2.00E-06 6.43E-01 4.99E-01 1.35E-05 1.95E-05 
SMD10 1.57E+01 6.89E-01 3.86E-02 5.15E-02 8.01E-03 3.77E-03 
SMD11 7.98E-02 1.40E-01 6.95E-03 1.46E-02 2.24E-03 3.84E-03 
SMD12 1.99E-02 4.32E+00 6.43E-02 8.05E-02 5.93E-02 3.49E-02 
BLTP1 1.00E-06 1.00E-06 1.00E-06 1.00E-06 6.77E-04 1.00E+02 
BLTP2 1.00E-06 5.25E-06 1.97E-03 5.15E-03 4.70E-03 1.21E-02 
BLTP3 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
BLTP4 6.20E+01 1.00E-06 2.69E+01 1.00E-06 3.44E+01 1.00E-06 
BLTP5 1.00E-06 6.25E-06 2.68E-03 5.66E-03 3.83E-03 8.09E-03 
BLTP6 8.36E-01 3.11E-03 5.84E-03 1.00E-06 2.06E-02 1.54E-06 
BLTP7 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
BLTP8 1.31E-02 1.00E-06 5.24E-05 1.00E-06 1.78E-04 1.00E-06 
BLTP9 3.25E-06 3.25E-06 5.43E-03 1.90E-03 3.84E-02 1.34E-02 
BLTP10 1.00E-06 7.90E-05 1.00E-06 2.13E-03 1.00E-06 1.84E-03 
BLTP11 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
BLTP12 1.00E-06 1.00E-06 
BLTP13 1.87E-02 1.75E-03 2.89E-02 5.59E-03 8.11E-02 1.39E-02 
BLMANBLEABLEAQ
Prb NameULLLULLLULLL
SMD1 1.00E-06 1.00E-06 6.15E-06 4.29E-06 1.00E-06 1.00E-06 
SMD2 1.00E-06 1.00E-06 6.07E-06 1.26E-05 1.59E-05 2.50E-05 
SMD3 1.00E-06 1.00E-06 1.31E-05 1.10E-05 1.12E-05 7.60E-06 
SMD4 1.00E-06 1.00E-06 1.48E-05 9.20E-06 1.00E-06 4.37E-06 
SMD5 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
SMD6 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
SMD7 1.00E-06 1.00E-06 1.45E-05 1.69E-06 9.75E-02 2.43E+02 
SMD8 4.97E-03 1.25E-03 1.31E-04 4.01E-05 4.80E-04 1.19E-04 
SMD9 1.00E-06 2.00E-06 6.43E-01 4.99E-01 1.35E-05 1.95E-05 
SMD10 1.57E+01 6.89E-01 3.86E-02 5.15E-02 8.01E-03 3.77E-03 
SMD11 7.98E-02 1.40E-01 6.95E-03 1.46E-02 2.24E-03 3.84E-03 
SMD12 1.99E-02 4.32E+00 6.43E-02 8.05E-02 5.93E-02 3.49E-02 
BLTP1 1.00E-06 1.00E-06 1.00E-06 1.00E-06 6.77E-04 1.00E+02 
BLTP2 1.00E-06 5.25E-06 1.97E-03 5.15E-03 4.70E-03 1.21E-02 
BLTP3 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
BLTP4 6.20E+01 1.00E-06 2.69E+01 1.00E-06 3.44E+01 1.00E-06 
BLTP5 1.00E-06 6.25E-06 2.68E-03 5.66E-03 3.83E-03 8.09E-03 
BLTP6 8.36E-01 3.11E-03 5.84E-03 1.00E-06 2.06E-02 1.54E-06 
BLTP7 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
BLTP8 1.31E-02 1.00E-06 5.24E-05 1.00E-06 1.78E-04 1.00E-06 
BLTP9 3.25E-06 3.25E-06 5.43E-03 1.90E-03 3.84E-02 1.34E-02 
BLTP10 1.00E-06 7.90E-05 1.00E-06 2.13E-03 1.00E-06 1.84E-03 
BLTP11 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 1.00E-06 
BLTP12 1.00E-06 1.00E-06 
BLTP13 1.87E-02 1.75E-03 2.89E-02 5.59E-03 8.11E-02 1.39E-02 

4.5  Performance Profile

Apart from reporting the performance on individual problems above, we also present a statistical comparison using performance profiles in Figure 5. The function is the cumulative distribution function for performance measure . The value effectively represents the factor (ratio) of the performance compared to the best performance. Thus the minimum value for is 1 (or 0 if plotted in log scale). In our study, the median accuracy ( and ) is used as the performance measure and is accordingly the ratio of median accuracy obtained by an algorithm to the best accuracy obtained for that problem, plotted in log scale. The resulting performance profile of the algorithms BLMA, NBLEA, BLEAQ, and BIDE are shown in Figure 5 considering the 25 problems (SMD and BLTP) for both levels. It is to be noted that for the case of SMD problems, BIDE is compared for only the first six problems (SMD1–SMD6) for which results are available. Therefore the values do not reach 1 for BIDE.

Figure 5:

Performance profile for NBLEA, BLEAQ, BIDE, and BLMA on the set of 25 benchmark problems. The value of is plotted on log scale.

Figure 5:

Performance profile for NBLEA, BLEAQ, BIDE, and BLMA on the set of 25 benchmark problems. The value of is plotted on log scale.

From Figure 5, it can be seen that the profile of BLMA always lies to the left of other profiles for any given value of . Thus it can be seen that the cumulative area under the curve (AUC) is the highest for it, for both upper and lower levels, establishing it as the best performing algorithm among the four considered. As far as the upper-level median accuracy is concerned, for 72% of the problems, BLMA is the best performing algorithm (). BIDE is the next best algorithm up to , beyond which the performance of BLEAQ is comparable to it. Beyond , both NBLEA and BLEAQ perform better than BIDE. The trends for the lower-level objective are quite similar. BLMA performs the best for 72% of the problems. NBLEA and BLEAQ perform worse than BIDE up to , after which their relative performance improves.

4.6  Sensitivity of BLMA Performance to

Most of the parameters used in the proposed approach are conventional DE or IP parameters which have standard interpretations in the literature. For example, an increase in population size and generations is likely to improve diversity and convergence at the cost of more computational effort (function evaluations). The settings used for upper-and lower-level populations and generations have been selected so that a comparison could be made with similar function evaluations as reported in the previous studies, although as discussed in Section 3, this is difficult to achieve. However, the reliability of the proposed approach is reflected in the fact that the parameters are not varied across problems, which means that they do not need significant tuning.

The key parameter that is additionally introduced in BLMA is , which determines the number of generations () after which the lower-level search switches from global to local. The sensitivity of the proposed algorithm to this parameter is of interest to establish reliability of the algorithm. The default value of has been set to 0.8 for the experiments presented so far. In this section, we present the results obtained by varying the value of  (while keeping other parameters unchanged). The results using are shown in Tables 11 (accuracy) and 12 (function evaluations). It is seen that the variations in the results are very marginal, indicating that for this set of problems the performance of BLMA is not too sensitive to . The Wilcoxon signed rank test with 95% confidence is also conducted for each of the cases compared to , and results are indicated using (+), (-) or (=) for higher, lower, or equivalent values respectively.

Table 11:
Median accuracy () obtained using different values of .
PrbUL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()
SMD1 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD2 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD3 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD4 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD5 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD6 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD7 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD8 4.55E-03 8.15E-04 3.57E-03(=) 6.35E-04(=) 3.82E-03(=) 8.15E-04(=) 3.70E-03(=) 7.42E-04(=) 
SMD9 2.00E-06 2.00E-06 3.00E-06(+) 3.00E-06(+) 2.00E-06(=) 2.00E-06(=) 2.00E-06(=) 2.00E-06(=) 
SMD10 4.00E-06 5.00E-06 4.00E-06(=) 4.00E-06(=) 4.00E-06(=) 5.00E-06(=) 1.60E-05(=) 4.00E-06(=) 
SMD11 2.89E-02 7.33E-03 1.12E-02(-) 1.00E-06(-) 5.56E-03(-) 1.00E-06(-) 5.56E-03(-) 1.00E-06(-) 
SMD12 4.00E-06 1.24E-02 4.00E-06(=) 3.71E-03(=) 4.00E-06(=) 1.19E-02(=) 4.00E-06(=) 1.29E-02(=) 
BLTP1 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
BLTP2 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 2.00E-06(=) 1.00E-06(=) 3.00E-06(-) 
BLTP3 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
BLTP4 1.00E-06 1.00E-06 1.00E-06(-) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
BLTP5 8.00E-06 8.80E-05 8.00E-06(=) 8.80E-05(=) 7.00E-06(=) 8.40E-05(-) 8.00E-06(=) 8.80E-05(=) 
BLTP6 7.48E-01 1.39E-03 4.66E-01(=) 5.39E-04(=) 9.45E-01(=) 2.21E-03(=) 3.95E-01(-) 3.89E-04(-) 
BLTP7 1.20E-05 9.00E-06 1.20E-05(=) 9.00E-06(=) 1.20E-05(=) 9.00E-06(=) 1.20E-05(=) 9.00E-06(=) 
BLTP8 3.06E-03 1.00E-06 4.17E-03(=) 1.00E-06(=) 4.11E-03(=) 1.00E-06(=) 3.22E-03(=) 1.00E-06(=) 
BLTP9 5.00E-06 1.60E-05 5.00E-06(=) 1.60E-05(=) 6.00E-06(=) 1.60E-05(=) 7.00E-06(=) 1.60E-05(=) 
BLTP10 1.00E-06 1.02E-04 1.00E-06(=) 1.04E-04(=) 1.00E-06(=) 1.03E-04(=) 1.00E-06(=) 9.70E-05(=) 
BLTP11 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
BLTP12 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(+) 1.00E-06(+) 2.00E-06(+) 6.00E-06(+) 
BLTP13 7.65E-03 3.10E-04 8.58E-03(=) 1.25E-03(=) 8.58E-03(+) 1.25E-03(+) 8.56E-03(=) 1.18E-03(=) 
PrbUL ()LL ()UL ()LL ()UL ()LL ()UL ()LL ()
SMD1 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD2 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD3 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD4 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD5 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD6 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD7 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
SMD8 4.55E-03 8.15E-04 3.57E-03(=) 6.35E-04(=) 3.82E-03(=) 8.15E-04(=) 3.70E-03(=) 7.42E-04(=) 
SMD9 2.00E-06 2.00E-06 3.00E-06(+) 3.00E-06(+) 2.00E-06(=) 2.00E-06(=) 2.00E-06(=) 2.00E-06(=) 
SMD10 4.00E-06 5.00E-06 4.00E-06(=) 4.00E-06(=) 4.00E-06(=) 5.00E-06(=) 1.60E-05(=) 4.00E-06(=) 
SMD11 2.89E-02 7.33E-03 1.12E-02(-) 1.00E-06(-) 5.56E-03(-) 1.00E-06(-) 5.56E-03(-) 1.00E-06(-) 
SMD12 4.00E-06 1.24E-02 4.00E-06(=) 3.71E-03(=) 4.00E-06(=) 1.19E-02(=) 4.00E-06(=) 1.29E-02(=) 
BLTP1 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
BLTP2 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 2.00E-06(=) 1.00E-06(=) 3.00E-06(-) 
BLTP3 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
BLTP4 1.00E-06 1.00E-06 1.00E-06(-) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
BLTP5 8.00E-06 8.80E-05 8.00E-06(=) 8.80E-05(=) 7.00E-06(=) 8.40E-05(-) 8.00E-06(=) 8.80E-05(=) 
BLTP6 7.48E-01 1.39E-03 4.66E-01(=) 5.39E-04(=) 9.45E-01(=) 2.21E-03(=) 3.95E-01(-) 3.89E-04(-) 
BLTP7 1.20E-05 9.00E-06 1.20E-05(=) 9.00E-06(=) 1.20E-05(=) 9.00E-06(=) 1.20E-05(=) 9.00E-06(=) 
BLTP8 3.06E-03 1.00E-06 4.17E-03(=) 1.00E-06(=) 4.11E-03(=) 1.00E-06(=) 3.22E-03(=) 1.00E-06(=) 
BLTP9 5.00E-06 1.60E-05 5.00E-06(=) 1.60E-05(=) 6.00E-06(=) 1.60E-05(=) 7.00E-06(=) 1.60E-05(=) 
BLTP10 1.00E-06 1.02E-04 1.00E-06(=) 1.04E-04(=) 1.00E-06(=) 1.03E-04(=) 1.00E-06(=) 9.70E-05(=) 
BLTP11 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 1.00E-06(=) 
BLTP12 1.00E-06 1.00E-06 1.00E-06(=) 1.00E-06(=) 1.00E-06(+) 1.00E-06(+) 2.00E-06(+) 6.00E-06(+) 
BLTP13 7.65E-03 3.10E-04 8.58E-03(=) 1.25E-03(=) 8.58E-03(+) 1.25E-03(+) 8.56E-03(=) 1.18E-03(=) 
Table 12:
Median function evaluations (UL FE, LL FE) obtained using different values of .
Prb.UL FELL FEUL FELL FEUl FELL FEUl FELL FE
SMD1 4.12E+02 3.05E+05 4.12E+02 2.63E+05 4.12E+02 2.21E+05 4.12E+02 1.79E+05 
SMD2 4.24E+02 3.01E+05 4.24E+02 2.60E+05 4.24E+02 2.17E+05 4.24E+02 1.76E+05 
SMD3 4.12E+02 3.09E+05 4.12E+02 2.69E+05 4.12E+02 2.29E+05 4.12E+02 1.90E+05 
SMD4 5.52E+02 3.29E+05 5.52E+02 2.86E+05 5.52E+02 2.43E+05 5.52E+02 1.99E+05 
SMD5 5.52E+02 3.28E+05 5.52E+02 2.86E+05 5.52E+02 2.44E+05 5.52E+02 1.99E+05 
SMD6 4.88E+02 3.26E+05 5.08E+02 2.89E+05 5.20E+02 2.48E+05 4.92E+02 2.05E+05 
SMD7 4.24E+02 3.01E+05 4.24E+02 2.58E+05 4.24E+02 2.16E+05 4.24E+02 1.74E+05 
SMD8 5.52E+02 3.25E+05 5.52E+02 2.84E+05 5.51E+02 2.39E+05 5.51E+02 1.95E+05 
SMD9 6.78E+02 5.30E+05 6.61E+02 4.38E+05 6.71E+02 3.47E+05 6.76E+02 2.57E+05 
SMD10 5.85E+02 5.13E+05 5.87E+02 4.19E+05 5.89E+02 3.24E+05 5.87E+02 2.28E+05 
SMD11 6.06E+02 5.12E+05 6.02E+02 4.17E+05 6.07E+02 3.21E+05 5.97E+02 2.25E+05 
SMD12 5.84E+02 5.10E+05 5.81E+02 4.14E+05 5.85E+02 3.17E+05 5.79E+02 2.19E+05 
BLTP1 3.04E+02 2.81E+05 3.04E+02 2.32E+05 3.04E+02 1.82E+05 3.04E+02 1.33E+05 
BLTP2 3.93E+02 2.84E+05 3.87E+02 2.34E+05 3.90E+02 1.85E+05 3.84E+02 1.35E+05 
BLTP3 3.03E+02 2.80E+05 3.03E+02 2.30E+05 3.03E+02 1.80E+05 3.03E+02 1.30E+05 
BLTP4 3.39E+02 2.87E+05 3.39E+02 2.39E+05 3.39E+02 1.91E+05 3.39E+02 1.43E+05 
BLTP5 3.87E+02 2.85E+05 3.86E+02 2.36E+05 3.84E+02 1.88E+05 3.87E+02 1.40E+05 
BLTP6 3.09E+02 2.81E+05 3.09E+02 2.31E+05 3.09E+02 1.82E+05 3.09E+02 1.32E+05 
BLTP7 3.38E+02 2.82E+05 3.38E+02 2.32E+05 3.38E+02 1.83E+05 3.37E+02 1.33E+05 
BLTP8 3.70E+02 2.88E+05 3.68E+02 2.41E+05 3.60E+02 1.94E+05 3.60E+02 1.49E+05 
BLTP9 3.88E+02 2.92E+05 3.94E+02 2.44E+05 3.89E+02 1.96E+05 3.87E+02 1.49E+05 
BLTP10 4.24E+02 2.82E+05 4.27E+02 2.32E+05 4.36E+02 1.83E+05 4.40E+02 1.34E+05 
BLTP11 3.11E+02 2.83E+05 3.11E+02 2.35E+05 3.11E+02 1.86E+05 3.11E+02 1.38E+05 
BLTP12 3.18E+02 2.84E+05 3.18E+02 2.37E+05 3.18E+02 1.89E+05 3.59E+02 1.43E+05 
BLTP13 4.51E+02 2.91E+05 4.72E+02 2.43E+05 4.53E+02 1.96E+05 4.37E+02 1.47E+05 
Prb.UL FELL FEUL FELL FEUl FELL FEUl FELL FE
SMD1 4.12E+02 3.05E+05 4.12E+02 2.63E+05 4.12E+02 2.21E+05 4.12E+02 1.79E+05 
SMD2 4.24E+02 3.01E+05 4.24E+02 2.60E+05 4.24E+02 2.17E+05 4.24E+02 1.76E+05 
SMD3 4.12E+02 3.09E+05 4.12E+02 2.69E+05 4.12E+02 2.29E+05 4.12E+02 1.90E+05 
SMD4 5.52E+02 3.29E+05 5.52E+02 2.86E+05 5.52E+02 2.43E+05 5.52E+02 1.99E+05 
SMD5 5.52E+02 3.28E+05 5.52E+02 2.86E+05 5.52E+02 2.44E+05 5.52E+02 1.99E+05 
SMD6 4.88E+02 3.26E+05 5.08E+02 2.89E+05 5.20E+02 2.48E+05 4.92E+02 2.05E+05 
SMD7 4.24E+02 3.01E+05 4.24E+02 2.58E+05 4.24E+02 2.16E+05 4.24E+02 1.74E+05 
SMD8 5.52E+02 3.25E+05 5.52E+02 2.84E+05 5.51E+02 2.39E+05 5.51E+02 1.95E+05 
SMD9 6.78E+02 5.30E+05 6.61E+02 4.38E+05 6.71E+02 3.47E+05 6.76E+02 2.57E+05 
SMD10 5.85E+02 5.13E+05 5.87E+02 4.19E+05 5.89E+02 3.24E+05 5.87E+02 2.28E+05 
SMD11 6.06E+02 5.12E+05 6.02E+02 4.17E+05 6.07E+02 3.21E+05 5.97E+02 2.25E+05 
SMD12 5.84E+02 5.10E+05 5.81E+02 4.14E+05 5.85E+02 3.17E+05 5.79E+02 2.19E+05 
BLTP1 3.04E+02 2.81E+05 3.04E+02 2.32E+05 3.04E+02 1.82E+05 3.04E+02 1.33E+05 
BLTP2 3.93E+02 2.84E+05 3.87E+02 2.34E+05 3.90E+02 1.85E+05 3.84E+02 1.35E+05 
BLTP3 3.03E+02 2.80E+05 3.03E+02 2.30E+05 3.03E+02 1.80E+05 3.03E+02 1.30E+05 
BLTP4 3.39E+02 2.87E+05 3.39E+02 2.39E+05 3.39E+02 1.91E+05 3.39E+02 1.43E+05 
BLTP5 3.87E+02 2.85E+05 3.86E+02 2.36E+05 3.84E+02 1.88E+05 3.87E+02 1.40E+05 
BLTP6 3.09E+02 2.81E+05 3.09E+02 2.31E+05 3.09E+02 1.82E+05 3.09E+02 1.32E+05 
BLTP7 3.38E+02 2.82E+05 3.38E+02 2.32E+05 3.38E+02 1.83E+05 3.37E+02 1.33E+05 
BLTP8 3.70E+02 2.88E+05 3.68E+02 2.41E+05 3.60E+02 1.94E+05 3.60E+02 1.49E+05 
BLTP9 3.88E+02 2.92E+05 3.94E+02 2.44E+05 3.89E+02 1.96E+05 3.87E+02 1.49E+05 
BLTP10 4.24E+02 2.82E+05 4.27E+02 2.32E+05 4.36E+02 1.83E+05 4.40E+02 1.34E+05 
BLTP11 3.11E+02 2.83E+05 3.11E+02 2.35E+05 3.11E+02 1.86E+05 3.11E+02 1.38E+05 
BLTP12 3.18E+02 2.84E+05 3.18E+02 2.37E+05 3.18E+02 1.89E+05 3.59E+02 1.43E+05 
BLTP13 4.51E+02 2.91E+05 4.72E+02 2.43E+05 4.53E+02 1.96E+05 4.37E+02 1.47E+05 

5  Practical Applications

Having established the efficacy of BLMA in solving the benchmark problems, in this section we demonstrate its performance on two practical application problems.

Both problems considered for this study are continuous bilevel problems with two objectives at each level. Since the presented BLMA is designed for single-objective problems only, a scalarization technique is required to convert them into corresponding single-objective bilevel problems. For scalarizing the lower-level problem, the adaptive scalarization technique discussed in Gupta and Ong (2015) is used. The process involves introducing a set of weights {} as additional variables at the upper level ( = number of lower-level objectives). Since the sum of weights should equal 1, the actual number of variables introduced is . Subsequently, for any given set of weights, the problem at the lower level is scalarized as using weighted sum (other type of scalarizations could also be used). For the bi-objective cases considered here, only one weight parameter () is required, and the lower-level objectives are thus modified as shown in Equation 6.
formula
6

Note that now acts as an upper-level variable and is evolved along with the . For the corresponding lower-level problem, and remain constant. This type of interpretation has also been proposed earlier in Dempe et al. (2013).

The next step is to scalarize the upper-level problem. This is done by using a fixed set of weights {}, where is the number of upper-level objectives. Once again, the sum of weights must be 1, so the number of actual parameters introduced are . Thus, the bi-objective upper-level problem is reformulated as shown in Equation 7. It is to be noted that these weights are not evolved during the algorithm run. For the studies presented here, is fixed to 0.5.
formula
7

By using the above scalarization technique, the problems at both levels are converted to single-objective problems, and the BLMA algorithm should ideally deliver one solution on the Pareto Optimal Front (POF) of the problem. The experiments on these scalarized problems are presented in the following subsections.

For solving the problem using BLMA, a population size of 20 is used at both levels. The number of generations used at the lower-level is 50 for both problems. The number of generations used for the upper-level is 5 and 50 for applications 1 and 2, respectively. Correspondingly, the population size used for NBLEA and BLEAQ is also set as 20 at both levels. The remaining parameters are the same as those used in the previous section. Twenty-nine independent runs were conducted using each algorithm.

5.1  Application 1: Gold Mining in Kuusamo

This application deals with an environmental problem arising out of the Kuusamo region of Finland, as described in Sinha et al. (2016). The region has a rich amount of gold deposits which could bring significant economic benefits, but could at the same time cause environmental damage. Therefore, the government needs to decide whether to allow mining and to what extent. Thus, the government has two objectives: maximize revenue and minimize the resulting effect on the environment. On the other hand, the mining company’s objectives are to maximize its profits, and at the same time its reputation (as excessive profits at the cost of environment may lead to poor public image). Considering the factors that influence the above objectives, the problem statement is mathematically represented as follows (Sinha et al., 2016):
formula
8
Here denotes per unit tax imposed on the mine and represents the amount of metal extracted by the mine.

The median solutions obtained using different algorithms are shown in Table 13 and Figure 6. Since the two objectives are of different order of magnitudes, the overall upper-level objective using is predominantly dictated by . As a result, the target solution (Maximum ) lies close to the extremity of the Pareto front where is maximum, as shown in Figure 6(a). It is seen that for this problem all of the algorithms are able to converge very close to this target solution. Further, the lower-level solution obtained also lies on the corresponding lower-level Pareto front as shown in Figure 6(b). This confirms that the solutions obtained by each of the algorithms satisfy the lower-level optimality. The Pareto optimal fronts plotted in Figure 6 are merely for the verification of optimality and are obtained using a separate, multi-objective heuristic, the discussion of which is not considered within scope of this study. The fronts are consistent with those presented in the original reference from which the application problems have been taken (Sinha et al., 2016).

Figure 6:

Median solutions obtained for application 1 (Pareto fronts are approximate).

Figure 6:

Median solutions obtained for application 1 (Pareto fronts are approximate).

Table 13:
Median results for Application 1.
Func EvalUL VarLL VarUL ObjLL Obj
Algo.ULLL
BLMA 1.93E+02 1.37E+05 5.05E+01 9.75E-01 1.21E+01 6.12E+02 −1.21E+01 3.00E+02 2.94E+02 −1.21E+01 2.87E+02 
NBLEA 5.07E+02 4.78E+03 5.00E+01 1.00E+00 1.23E+01 6.13E+02 −1.23E+01 3.00E+02 3.00E+02 −1.23E+01 3.00E+02 
BLEAQ 7.36E+02 1.77E+03 5.00E+01 1.00E+00 1.22E+01 6.12E+02 −1.22E+01 3.00E+02 3.00E+02 −1.22E+01 3.00E+02 
DE-DE 1.40E+03 1.96E+06 5.00E+01 1.00E+00 1.22E+01 6.12E+02 −1.22E+01 3.00E+02 3.00E+02 −1.22E+01 3.00E+02 
DE-IP 1.40E+03 5.02E+05 5.00E+01 1.00E+00 1.22E+01 6.12E+02 −1.22E+01 3.00E+02 3.00E+02 −1.22E+01 3.00E+02 
IP-IP 1.51E+03 5.23E+05 5.00E+01 9.97E-01 1.22E+01 6.12E+02 −1.22E+01 3.00E+02 3.00E+02 −1.22E+01 2.99E+02 
IP-DE 1.43E+03 1.97E+06 5.05E+01 9.99E-01 1.21E+01 6.12E+02 −1.21E+01 3.00E+02 2.94E+02 −1.21E+01 2.94E+02 
Func EvalUL VarLL VarUL ObjLL Obj
Algo.ULLL
BLMA 1.93E+02 1.37E+05 5.05E+01 9.75E-01 1.21E+01 6.12E+02 −1.21E+01 3.00E+02 2.94E+02 −1.21E+01 2.87E+02 
NBLEA 5.07E+02 4.78E+03 5.00E+01 1.00E+00 1.23E+01 6.13E+02 −1.23E+01 3.00E+02 3.00E+02 −1.23E+01 3.00E+02 
BLEAQ 7.36E+02 1.77E+03 5.00E+01 1.00E+00 1.22E+01 6.12E+02 −1.22E+01 3.00E+02 3.00E+02 −1.22E+01 3.00E+02 
DE-DE 1.40E+03 1.96E+06 5.00E+01 1.00E+00 1.22E+01 6.12E+02 −1.22E+01 3.00E+02 3.00E+02 −1.22E+01 3.00E+02 
DE-IP 1.40E+03 5.02E+05 5.00E+01 1.00E+00 1.22E+01 6.12E+02 −1.22E+01 3.00E+02 3.00E+02 −1.22E+01 3.00E+02 
IP-IP 1.51E+03 5.23E+05 5.00E+01 9.97E-01 1.22E+01 6.12E+02 −1.22E+01 3.00E+02 3.00E+02 −1.22E+01 2.99E+02 
IP-DE 1.43E+03 1.97E+06 5.05E+01 9.99E-01 1.21E+01 6.12E+02 −1.21E+01 3.00E+02 2.94E+02 −1.21E+01 2.94E+02 

5.2  Application 2: Decision Making in a Company

This problem represents a scenario of a company operating with different branches under it (Sinha et al., 2016). The objectives of the chief executive officer (CEO) of a company are to maximize the profit and the quality of products. On the other hand, for the branch heads, the objectives are to maximize the branch profits and their workers’ satisfaction. Consequently, this is a bi-objective bilevel optimization problem. There are four constraints at the upper level and five at the lower level. The mathematical formulation of the problem is given in Equation 9.
formula
9

The median solutions using each of the algorithms are shown in Table 14 and Figure 7. The upper-level Pareto front for the problem contains a “knee” which is at the intersection of two linear segments. For , this knee point has the maximum and hence is the target solution, shown in Figure 7(a). Three of the algorithms, BLMA, DE-DE, and DE-IP were able to converge close to this target solution. The solutions obtained by NBLEA, IP-IP, and IP-DE were relatively farther from the upper-level front. Each of these solutions seems valid, since the lower-level solution lies on (or very close to) the corresponding Pareto optimal front as seen in Figure 7(b). BLEAQ was not able to achieve any feasible solution, and was hence excluded from the plot.

Table 14:
Median results for Application 2.
Func EvalUL VarLL VarUL ObjLL Obj
Algo.ULLL
BLMA 1.13E+03 1.08E+06 1.24E+02 3.81E+01 6.33E-01 1.90E-01 1.20E+01 1.48E+01 7.64E+02 1.10E+03 9.32E+02 8.15E+02 9.66E+02 8.70E+02 
NBLEA 5.73E+02 1.54E+04 1.10E+02 1.45E+01 8.58E-01 2.93E+01 5.74E-09 1.03E+01 7.87E+02 9.69E+02 8.78E+02 4.65E+02 9.77E+02 5.38E+02 
BLEAQ 
DE-DE 1.40E+03 1.96E+06 1.24E+02 3.80E+01 1.00E+00 0.00E+00 1.24E+01 1.41E+01 7.57E+02 1.10E+03 9.27E+02 8.08E+02 9.61E+02 8.08E+02 
DE-IP 1.40E+03 2.04E+06 1.23E+02 3.69E+01 9.90E-01 2.03E+00 1.08E+01 1.50E+01 7.70E+02 1.09E+03 9.30E+02 7.98E+02 9.65E+02 8.00E+02 
IP-IP 1.44E+03 2.67E+06 1.10E+02 8.69E-01 1.00E+00 4.18E+01 1.76E+00 6.81E-01 7.36E+02 9.42E+02 8.39E+02 2.34E+02 1.00E+03 2.34E+02 
IP-DE 1.44E+03 1.98E+06 1.06E+02 9.57E-01 4.89E-01 4.41E+01 0.00E+00 3.33E+00 7.59E+02 9.26E+02 8.43E+02 2.29E+02 1.00E+03 6.25E+02 
Func EvalUL VarLL VarUL ObjLL Obj
Algo.ULLL
BLMA 1.13E+03 1.08E+06 1.24E+02 3.81E+01 6.33E-01 1.90E-01 1.20E+01 1.48E+01 7.64E+02 1.10E+03 9.32E+02 8.15E+02 9.66E+02 8.70E+02 
NBLEA 5.73E+02 1.54E+04 1.10E+02 1.45E+01 8.58E-01 2.93E+01 5.74E-09 1.03E+01 7.87E+02 9.69E+02 8.78E+02 4.65E+02 9.77E+02 5.38E+02 
BLEAQ 
DE-DE 1.40E+03 1.96E+06 1.24E+02 3.80E+01 1.00E+00 0.00E+00 1.24E+01 1.41E+01 7.57E+02 1.10E+03 9.27E+02 8.08E+02 9.61E+02 8.08E+02 
DE-IP 1.40E+03 2.04E+06 1.23E+02 3.69E+01 9.90E-01 2.03E+00 1.08E+01 1.50E+01 7.70E+02 1.09E+03 9.30E+02 7.98E+02 9.65E+02 8.00E+02 
IP-IP 1.44E+03 2.67E+06 1.10E+02 8.69E-01 1.00E+00 4.18E+01 1.76E+00 6.81E-01 7.36E+02 9.42E+02 8.39E+02 2.34E+02 1.00E+03 2.34E+02 
IP-DE 1.44E+03 1.98E+06 1.06E+02 9.57E-01 4.89E-01 4.41E+01 0.00E+00 3.33E+00 7.59E+02 9.26E+02 8.43E+02 2.29E+02 1.00E+03 6.25E+02 
Figure 7:

Median solutions obtained for application 2 (Pareto fronts are approximate).

Figure 7:

Median solutions obtained for application 2 (Pareto fronts are approximate).

6  Conclusions

Bilevel optimization is a challenging problem with practical applications in several fields. In this article, major approaches used for solving bilevel optimization problems are reviewed. Thereafter, a Bilevel Memetic Algorithm (BLMA) is proposed to solve single-objective, continuous nonlinear problems, which attempts to extend the advantages of hybridization of global and local search to bilevel optimization problems. In the proposed BLMA, both upper and lower levels use either a global or a local search, depending on which phase the search is in. Further, a “re-evaluation” step is performed on selected promising solutions in order to improve confidence on lower-level optimality. These provisions enable the algorithm to solve for global optimum at both levels reliably, without requiring an excessive number of function evaluations. The algorithm is tested on an extensive set of benchmark problems, and compared with a set of baseline nested strategies and other established algorithms, such as NBLEA, BLEAQ, and BIDE. The efficiency of BLMA is demonstrated through competitive performance on the benchmark problems, within a reasonable number of function evaluations. Further, BLMA is also used to solve a scalarized version of two practical bilevel optimization problems. The numerical experiments highlight the significant potential of BLMA in solving single-objective bilevel optimization problems.

References

Aiyoshi
,
E.
, and
Shimizu
,
K
. (
1984
).
A solution method for the static constrained Stackelberg problem via penalty method
.
IEEE Transactions on Automatic Control
,
29
(
12
):
1111
1114
.
Anandalingam
,
G.
, and
White
,
D
. (
1990
).
A solution method for the linear static Stackelberg problem using penalty functions
.
IEEE Transactions on Automatic Control
,
35
(
10
):
1170
1173
.
Angelo
,
J. S.
,
Krempser
,
E.
, and
Barbosa
,
H. J.
(
2013
).
Differential evolution for bilevel programming
.
IEEE Congress on Evolutionary Computation
,
470
477
.
Barbosa
,
H. J.
,
Bernardino
,
H. S.
, and
Barreto
,
A.
(
2010
).
Using performance profiles to analyze the results of the 2006 CEC constrained optimization competition
.
IEEE Congress on Evolutionary Computation
,
1
8
.
Bard
,
J. F
. (
1983
).
An efficient point algorithm for a linear two-stage optimization problem
.
Operations Research
,
31
(
4
):
670
684
.
Bard
,
J. F
. (
1988
).
Convex two-level optimization
.
Mathematical Programming
,
40
(
1–3
):
15
27
.
Bard
,
J. F.
, and
Falk
,
J. E
. (
1982
).
An explicit solution to the multi-level programming problem
.
Computers & Operations Research
,
9
(
1
):
77
100
.
Byrd
,
R. H.
,
Gilbert
,
J. C.
, and
Nocedal
,
J
. (
2000
).
A trust region method based on interior point techniques for nonlinear programming
.
Mathematical Programming
,
89
(
1
):
149
185
.
Byrd
,
R. H.
,
Hribar
,
M. E.
, and
Nocedal
,
J
. (
1999
).
An interior point algorithm for large-scale nonlinear programming
.
SIAM Journal on Optimization
,
9
(
4
):
877
900
.
Candler
,
W.
, and
Townsley
,
R
. (
1982
).
A linear two-level programming problem
.
Computers & Operations Research
,
9
(
1
):
59
76
.
Carrión
,
M.
,
Arroyo
,
J. M.
, and
Conejo
,
A. J
. (
2009
).
A bilevel stochastic programming approach for retailer futures market trading
.
IEEE Transactions on Power Systems
,
24
(
3
):
1446
1456
.
Colson
,
B.
,
Marcotte
,
P.
, and
Savard
,
G
. (
2005
).
A trust-region method for nonlinear bilevel programming: Algorithm and computational experience
.
Computational Optimization and Applications
,
30
(
3
):
211
227
.
Dempe
,
S
. (
1987
).
A simple algorithm for the-linear bilevel programming problem
.
Optimization
,
18
(
3
):
373
385
.
Dempe
,
S.
,
Gadhi
,
N.
, and
Zemkoho
,
A. B
. (
2013
).
New optimality conditions for the semivectorial bilevel optimization problem
.
Journal of Optimization Theory and Applications
,
157
(
1
):
54
74
.
Dolan
,
E. D.
, and
Moré
,
J. J
. (
2002
).
Benchmarking optimization software with performance profiles
.
Mathematical Programming
,
91
(
2
):
201
213
.
Falk
,
J. E.
, and
Liu
,
J
. (
1995
).
On bilevel programming, Part i: General nonlinear cases
.
Mathematical Programming
,
70
(
1–3
):
47
72
.
Gendreau
,
M.
,
Marcotte
,
P.
, and
Savard
,
G
. (
1996
).
A hybrid tabu-ascent algorithm for the linear bilevel programming problem
.
Journal of Global Optimization
,
8
(
3
):
217
233
.
Gong
,
W.
,
Cai
,
Z.
,
Ling
,
C. X.
, and
Li
,
C
. (
2011
).
Enhanced differential evolution with adaptive strategies for numerical optimization
.
IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics
,
41
(
2
):
397
413
.
Gupta
,
A.
, and
Ong
,
Y.-S.
(
2015
).
An evolutionary algorithm with adaptive scalarization for multiobjective bilevel programs
.
IEEE Congress on Evolutionary Computation
,
1636
1642
.
Halter
,
W.
, and
Mostaghim
,
S.
(
2006
).
Bilevel optimization of multi-component chemical systems using particle swarm optimization
.
IEEE Congress on Evolutionary Computation
,
1240
1247
.
Hansen
,
P.
,
Jaumard
,
B.
, and
Savard
,
G
. (
1992
).
New branch-and-bound rules for linear bilevel programming
.
SIAM Journal on Scientific and Statistical Computing
,
13
(
5
):
1194
1217
.
Hejazi
,
S. R.
,
Memariani
,
A.
,
Jahanshahloo
,
G.
, and
Sepehri
,
M. M
. (
2002
).
Linear bilevel programming solution by genetic algorithm
.
Computers & Operations Research
,
29
(
13
):
1913
1925
.
Herskovits
,
J.
,
Leontiev
,
A.
,
Dias
,
G.
, and
Santos
,
G
. (
2000
).
Contact shape optimization: A bilevel programming approach
.
Structural and Multidisciplinary Optimization
,
20
(
3
):
214
221
.
Islam
,
M. M.
,
Singh
,
H. K.
, and
Ray
,
T.
(
2015
).
A memetic algorithm for solving single objective bilevel optimization problems
.
2015 IEEE Congress on Evolutionary Computation
,
1643
1650
.
Jenabi
,
M.
,
Fatemi Ghomi
,
S.
, and
Smeers
,
Y
. (
2013
).
Bi-level game approaches for coordination of generation and transmission expansion planning within a market environment
.
IEEE Transactions on Power Systems
,
28
(
3
):
2639
2650
.
Kirjner-Neto
,
C.
,
Polak
,
E.
, and
Der Kiureghian
,
A
. (
1998
).
An outer approximations approach to reliability-based optimal design of structures
.
Journal of Optimization Theory and Applications
,
98
(
1
):
1
16
.
Koh
,
A.
(
2013
).
A metaheuristic framework for bi-level programming problems with multi-disciplinary applications
.
Metaheuristics for Bi-level Optimization
,
153
187
.
Kolstad
,
C. D.
, and
Lasdon
,
L. S
. (
1990
).
Derivative evaluation and computational experience with large bilevel mathematical programs
.
Journal of Optimization Theory and Applications
,
65
(
3
):
485
499
.
Kuo
,
R.
, and
Huang
,
C
. (
2009
).
Application of particle swarm optimization algorithm for solving bi-level linear programming problems
.
Computers & Mathematics with Applications
,
58
(
4
):
678
685
.
Legillon
,
F.
,
Liefooghe
,
A.
, and
Talbi
,
E.
(
2012
).
Cobra: A cooperative coevolutionary algorithm for bi-level optimization
.
IEEE Congress on Evolutionary Computation
,
1
8
.
Lv
,
Y.
,
Hu
,
T.
,
Wang
,
G.
, and
Wan
,
Z
. (
2007
).
A penalty function method based on Kuhn–Tucker condition for solving linear bilevel programming
.
Applied Mathematics and Computation
,
188
(
1
):
808
813
.
Marcotte
,
P.
, and
Zhu
,
D. L
. (
1996
).
Exact and inexact penalty methods for the generalized bilevel programming problem
.
Mathematical Programming
,
74
(
2
):
141
157
.
Mathieu
,
R.
,
Pittard
,
L.
, and
Anandalingam
,
G
. (
1994
).
Genetic algorithm based approach to bi-level linear programming
.
RAIRO. Recherche Opérationnelle
,
28
(
1
):
1
21
.
Migdalas
,
A
. (
1995
).
Bilevel programming in traffic planning: Models, methods and challenge
.
Journal of Global Optimization
,
7
(
4
):
381
405
.
Moscato
,
P.
(
1989
).
On evolution, search, optimization, genetic algorithms and martial arts: Towards memetic algorithms
.
Technical report, Caltech concurrent computation program, C3P Report 826
.
Moscato
,
P.
,
Cotta
,
C.
, and
Mendes
,
A.
(
2004
).
Memetic algorithms
.
New Optimization Techniques in Engineering
,
53
85
.
Oduguwa
,
V.
, and
Roy
,
R.
(
2002
).
Bi-level optimisation using genetic algorithm
.
IEEE International Conference on Artificial Intelligence Systems
,
322
327
.
Önal
,
H
. (
1993
).
A modified simplex approach for solving bilevel linear programming problems
.
European Journal of Operational Research
,
67
(
1
):
126
135
.
Rajesh
,
J.
,
Gupta
,
K.
,
Kusumakar
,
H. S.
,
Jayaraman
,
V. K.
, and
Kulkarni
,
B. D
. (
2003
).
A tabu search based approach for solving a class of bilevel programming problems in chemical engineering
.
Journal of Heuristics
,
9
(
4
):
307
319
.
Saharidis
,
G. K.
, and
Ierapetritou
,
M. G
. (
2009
).
Resolution method for mixed integer bi-level linear problems based on decomposition technique
.
Journal of Global Optimization
,
44
(
1
):
29
51
.
Sahin
,
K. H.
, and
Ciric
,
A. R
. (
1998
).
A dual temperature simulated annealing approach for solving bilevel programming problems
.
Computers & Chemical Engineering
,
23
(
1
):
11
25
.
Savard
,
G.
, and
Gauvin
,
J
. (
1994
).
The steepest descent direction for the nonlinear bilevel programming problem
.
Operations Research Letters
,
15
(
5
):
265
272
.
Singh
,
H. K.
,
Ray
,
T.
, and
Smith
,
W.
(
2011
).
Performance of infeasibility empowered memetic algorithm (IEMA) on engineering design problems
.
Advances in Artificial Intelligence
,
425
434
.
Sinha
,
A.
,
Malo
,
P.
, and
Deb
,
K.
(
2013a
).
Efficient evolutionary algorithm for single-objective bilevel optimization
. Retrieved from
arXiv preprint arXiv:1303.3901
.
Sinha
,
A.
,
Malo
,
P.
, and
Deb
,
K.
(
2014a
).
An improved bilevel evolutionary algorithm based on quadratic approximations
.
IEEE Congress on Evolutionary Computation
,
1870
1877
.
Sinha
,
A.
,
Malo
,
P.
, and
Deb
,
K
. (
2014b
).
Test problem construction for single-objective bilevel optimization
.
Evolutionary Computation
,
22
(
3
):
439
477
.
Sinha
,
A.
,
Malo
,
P.
,
Deb
,
K.
,
Korhonen
,
P.
, and
Wallenius
,
J.
(
2016
).
Solving bilevel multi-criterion optimization problems with lower level decision uncertainty
.
IEEE Transactions on Evolutionary Computation
,
20
(
2
):
199
217
.
Sinha
,
A.
,
Malo
,
P.
,
Frantsev
,
A.
, and
Deb
,
K.
(
2013
).
Multi-objective Stackelberg game between a regulating authority and a mining company: A case study in environmental economics
. In
IEEE Congress on Evolutionary Computation (CEC)
, pages
478
485
.
Sun
,
H.
,
Gao
,
Z.
, and
Wu
,
J
. (
2008
).
A bi-level programming model and solution algorithm for the location of logistics distribution centers
.
Applied Mathematical Modelling
,
32
(
4
):
610
616
.
Takahama
,
T.
, and
Sakai
,
S.
(
2006
).
Constrained optimization by the constrained differential evolution with gradient-based mutation and feasible elites
.
IEEE Congress on Evolutionary Computation
,
1
8
.
Vicente
,
L.
,
Savard
,
G.
, and
Júdice
,
J
. (
1994
).
Descent approaches for quadratic bilevel programming
.
Journal of Optimization Theory and Applications
,
81
(
2
):
379
399
.
Wang
,
D.
, and
Du
,
G
. (
2014
).
A self adaptive penalty function based genetic algorithm for value-bilevel programming problem
.
International Journal on Computer Science and Engineering
,
3
(
2
):
136
146
.
Wang
,
G.-m.
,
Wang
,
X.-j.
,
Wan
,
Z.-p.
, and
Jia
,
S.-h.
(
2007
).
An adaptive genetic algorithm for solving bilevel linear programming problems
.
Applied Mathematics and Mechanics
,
28
:
1605
1612
.
Wang
,
Y.
,
Jiao
,
Y.-C.
, and
Li
,
H
. (
2005
).
An evolutionary algorithm for solving nonlinear bilevel programming based on a new constraint-handling scheme
.
IEEE Transactions on Systems, Man, and Cybernetics, Part C: Applications and Reviews
,
35
(
2
):
221
232
.
White
,
D. J.
, and
Anandalingam
,
G
. (
1993
).
A penalty function approach for solving bi-level linear programs
.
Journal of Global Optimization
,
3
(
4
):
397
419
.
Ye
,
J. J
. (
2006
).
Constraint qualifications and kkt conditions for bilevel programming problems
.
Mathematics of Operations Research
,
31
(
4
):
811
824
.
Ye
,
J. J.
, and
Zhu
,
D
. (
2010
).
New necessary optimality conditions for bilevel programs by combining the mpec and value function approaches
.
SIAM Journal on Optimization
,
20
(
4
):
1885
1905
.
Zhu
,
X.
,
Yu
,
Q.
, and
Wang
,
X.
(
2006
).
A hybrid differential evolution algorithm for solving nonlinear bilevel programming with linear constraints
.
Proceedings of 5th IEEE International Conference on Cognitive Informatics
,
126
131
.

Notes

1 

Please note that based on prior studies, SMD9–12 are relatively complex nonlinear constrained problems and hence are assigned a higher number of evaluations compared to other problems, for all algorithms in this study.

2 

Please note that the current version of the code has certain modifications from the one used for the study in  Sinha et al. (2014b). This is discussed further in Section 4.4.

3 

Please note that population size used in Sinha et al. (2014b) and dimensionality of problems studied in Sinha, Malo, and Deb (2013) are different from those in this article, and hence there are differences in the reported numbers of function evaluations.