TY - JOUR AU - Ma, Mei AB - Introduction Problem models BLPP is a typical representative of multilevel hierarchical optimisation problems. In contrast to multi-objective optimisation, the decision makers in BLPP are at two different levels. This hierarchical structure often leads to a corresponding problem that is neither convex nor differentiable, which is also strongly non-deterministic polynomial hard as well. The BLPP model is expressed as follows: (1) Here, x = (x1, …, xn) and y = (y1, …, ym) are the leader’s and follower’s variables. F: Rn+m → R and f: Rn+m → R are the leader’s and follower’s objective function. G: Rn+m → R and g: Rn+m → R are the leader’s and follower’s constraints, respectively. It can be seen from (1) that BLPP is an interaction optimisation model between two decision makers with their own objectives. The decision-making procedure is executed as follows: The leader, located at upper level, makes a decision by selecting a variable value. Then the follower observes the follower’s selection and responds to the leader’s decision by optimizing his/her objective and providing an optimal solution y. The point pair (x, y) is called a feasible point of BLPP. The bilevel optimisation arms to select a value x that optimizes the leader’s objective among all feasible points. During the optimisation procedure, the leader influences the follower by providing a value x fixed as a parametric value. Subsequently, the follower alters the values of the objective and constraints at the leader’s problems by reacting a value y to address the leader’s problem. When the optimal solution (x, y) is obtained, the optimisation process stops. The nested nature of BLPP poses a number of additional challenges compared with conventional single-level optimization problems. In particular: The leader’s behaviour of a bilevel problem may be nonlinear, even if the problems at both levels are linear. It has been proved that BLPP are non-deterministic polynomial hard problem. Theoretically, a leader’s solution is considered valid/feasible only if the corresponding follower’s variables are the true global optimum of the follower problem. Global optimality can only be assured in very limited cases, such as convex and linear problems. However, for most nonlinear and black-box problems, it is not possible to ensure global optimality. In the deceptive case, incorrect follower’s optimal values may cause the objective value to be better than the leader’s true optimal value, which poses a severe challenge to the ranking strategy used in evolutionary optimization technology. Due to each variable of the leader needs to solve the follower programming problem to obtain a feasible solution of the BLPP, the cost of calculation is significantly high. In particularly, if the problem is mathematically not suited to be solved using exact techniques, thus, evolutionary/hybrid techniques are used. Related work With the advancements and developments in human science and technology, various increasingly complex optimisation models have emerged, particularly in the fields of equipment installation, task scheduling, production planning, line optimisation, software design, and tariff setting. Therefore, how these models can be effectively solved has emerged as an important main research topic in the optimisation field. Over the past few decades, many optimisation research results have mainly focused on single-objective and multi-objective optimisation models, whereas relatively few studies have been conducted on hierarchical optimisation models. Leveraging the one-time decision theory to produce multiple short-life cycle products, Zhu and Guo [1]studied BLPP applications with follower problem operators for manufacturers and used classical optimisation methods to solve them. Nasrolahpour et al. [2] developed an energy storage system for merchant pricing based on a two-tier complementary model that can determine the most favourable transaction behaviour. Under a multi-objective framework, Ahmad et al. [3] proposed a simple multi-objective bilevel linear programming method, which considered reservoir managers and multiple water-use departments as a hierarchical structure to optimally allocate limited water resources. More applications of optimisation models can be found in the literature [4–13]. Many practical applications have promoted theoretical research on bilevel programming problems, such as developing efficient algorithms and obtaining optimality conditions. However, owing to the computational complexity of the BLPP itself, it is often very challenging to adopt traditional optimisation methods based on gradient information to address such problems. Currently, only special bilevel programming approachs, such as linear and convex quadratic programs, can obtain the optimal solution to these problem through optimality conditions. For other types of bilevel programming problems research has focused on designing swarm intelligence and hybrid algorithms, which are currently more effective algorithm frameworks and have achieved better solutions when tested against certain problems. Existing methods for solving the BLPP can be grouped into the following categories: A). Classical approaches Some classic methods for BLPP include the simplex [14], branch-and-bound method [15], descending gradient method [16], and penalty function methods [17–19]. Classical methods typically apply the optimality conditions of the follower to convert the BLPP into a single-level problem. Dempe [14]used the simplex method to propose an algorithm for solving linear BLPP. The algorithm introduces slack variables to find a basis as a feasible solution of the BLPP, through a simplex and iterative method. Although this method is more effective for solving small-scale linear BLPP it cannot be directly extended to nonlinear problems. Susanne et al. [16] replaced the the follower problem a with Karush-Kuhn-Tucker condition and applied the optimal value method to transform the original problem into a single-level problems. This approach could effectively satisfy the complementary conditions in the mathematical programming problem in Banach space and introduce M stationarity. Using the dual gap as a penalty function, the literature of White [17] proposed an effective algorithm for solving the BLPP. It adopts a newly designed precise penalty function method to obtain the global optimal solution of the linear situation and conducts a novel theoretical analysis. In the literature [18, 19], a weak linear BLPP dealt with the penalty function and Karush-Kuhn-Tucker conditions. B). Evolutionary approaches The evolutionary algorithms, as representatives of the swarm intelligence optimisation technology, have been widely used to solve various BLPP over the past few decades. An early evolutionary algorithm was proposed by Mathieu et al. [20], who applied linear programming methods to handle followers’ problems and used a genetic algorithm to explore the search space of leaders’ problems. Focused on BLPPs whose follower’s problem is convex programming, Wang [21] proposed an evolutionary algorithm with an embedded constraint processing method. This method applies Karush-Kuhn-Tucker conditions to transform the followers’ problems, turning the original model into a single-level programming model. In addition, to obtain sufficient feasible solutions, a constraint-processing method was designed. Based on similar optimality conditions, Li [22] presented a genetic algorithm that solves nonlinear/linear fractional BLPP. In another study on the presence of the so-called pseudo-feasible solutions in evolutionary bilevel optimisation(QBCA-2), the focus was on determining how pseudo-feasible solutions can affect the performance of an evolutionary algorithm. Moreover, a novel and scalable set of test problems with characterised pseudo-feasible solutions was introduced in [23]. In the literature [24], a co-evolutionary algorithm was proposed to solve BLPP. In the evolution process, the follower problem was solved in two stages. Aboelnaga et al. [25] proposed an improved genetic algorithm and chaotic search method. Goshu et al. [26] proposed a metaheuristic algorithm for random BLPP. An evolutionary algorithm for solving nonlinear bilevel programming problems has been presented in the literature [27]. The algorithm is designed by reflecting the optimal solution of the follower problem to the leader problem. To ensure the quality of each iteration, the algorithm adaptive changes the population size during the evolution process and generates individuals using the tabu search method. C). Hybrid approaches The hybrid algorithm [28] is a common method used to solve BLPP. Abo-Elnaga et al. [29] proposed a multi-sine-cosine algorithm to solve nonlinear BLPP. The sine-cosine algorithms based on three different populations were presented. The first population was used to deal with the leader programming problem, while the second one addressed the follower programming problem. In addition, the Karush-Kuhn-Tucker condition was applied to transform the initial problem into a constrained optimisation problem, which was solved using the third population. If the objective function value is equal to zero, then the solution obtained by solving the leader and follower problems is deemed feasible. Wang [30] proposed a particle swarm distribution estimation algorithm with an embedded Estimation of Distribution Algorithm to solve nonlinear BLPP. Before executing the speed and location update rules in Partical Swarm Optimization, a Gaussian distribution was applied to generate offspring to replace some inferior individuals (particles) in the current population. D). EA based on approximate methods To avoid lengthy calculation processes, some evolutionary algorithms utilise approximation techniques to improve the efficiency of the intermediate calculation process. A bilayer covariance matrix adaptive evolution strategy(BLCMAES) is proposed in [31]. The method designed a sharing mechanism so that prior knowledge of the follower problem can be extracted from the leader optimizer, reduced the number of evaluations of the follower problem. Furthermore, an optimization-based elite retention mechanism is proposed to keep track of the elite and avoid incorrect solutions. Sinha et al. [32] proposed an evolutionary algorithm that uses approximate functions to adress BLPPs. For offspring generated by evolutionary operators, the approximation method can reduce the number of evaluations of the follower objective function. Using approximate Karush-Kuhn-Tucker conditions, Sinha et al. [33, 34] transformed the BLPP into a single-level optimisation problem; then leveraged an evolutionary algorithm embedded with the idea of neighbourhood measurements for the transformed model. Sinha et al. [35] presented an evolutionary optimisation algorithm (BLEAQ-2) that fits an extreme value mapping using the relationship between the leader and follower variables. Islam et al. [36] introduced an evolutionary algorithm that employs three types of surrogate models to approximate the optimal solution of the follower problem. Based on the linear programming optimality conditions, Li [37] proposed a genetic algorithm for solving linear BLPP. In this method, the follower is adopted as the search object of the evolutionary algorithm, and the lower variables in the leader problems are replaced by possible solution functions. This process locally optimises the leader variables. Research motivation In the field of engineering and economic management, optimization models involving different decision-making levels often appear, such problems are called multi-level optimization problems. Since the decision-making process of the problem is hierarchical, it is also called a hierarchical optimization problem. BLPP is a typical representative of multi-level hierarchical optimization problems, and has become an important research field for optimization problems due to its extensive practical application background and algorithmic challenges. Different from multi-objective optimization, the decision makers of BLPP are at two different levels, and this hierarchical structure often leads to problems that are non-convex and non-differentiable. Based on the above characteristics, it is often difficult for traditional optimization algorithms based on gradients to find the global optimal solution to the BLPP. Evolutionary algorithm be increasingly used to solve BLPP because it have the characteristics of global convergence, and there is no restriction on functions that are convex and differentiable. In this paper, driven the correlation coefficient and surrogate model, an evolutionary algorithm (TCEA) is proposed to solve complex nonlinear BLPP. The algorithm has the following characteristics: First, the isodata clustering method [38] is used to group the initial populations, and then determine the correlation coefficients of the leader and follower objective functions in each group based on the rank of the leader (follower) objective function, after that, some points are selected and updated in the offsprings based on the correlation coefficient value in each group. Second, for the offspring individuals produced by the crossover and mutation operators, the surrogate model is used to approximate the solution of the follower programming problem, and so does to reduce the number of evaluations of follower problems. The remainder of this paper is organized as follows. The basic concepts of BLPP are described in the next Section. The correlation coefficients and surrogate models are presented in Section 3. A new evolutionary algorithm based on the above methods are stated in Section 4. Experimental results and comparisons are provided in Section 5. We conclude our approach in Section 6. Problem models BLPP is a typical representative of multilevel hierarchical optimisation problems. In contrast to multi-objective optimisation, the decision makers in BLPP are at two different levels. This hierarchical structure often leads to a corresponding problem that is neither convex nor differentiable, which is also strongly non-deterministic polynomial hard as well. The BLPP model is expressed as follows: (1) Here, x = (x1, …, xn) and y = (y1, …, ym) are the leader’s and follower’s variables. F: Rn+m → R and f: Rn+m → R are the leader’s and follower’s objective function. G: Rn+m → R and g: Rn+m → R are the leader’s and follower’s constraints, respectively. It can be seen from (1) that BLPP is an interaction optimisation model between two decision makers with their own objectives. The decision-making procedure is executed as follows: The leader, located at upper level, makes a decision by selecting a variable value. Then the follower observes the follower’s selection and responds to the leader’s decision by optimizing his/her objective and providing an optimal solution y. The point pair (x, y) is called a feasible point of BLPP. The bilevel optimisation arms to select a value x that optimizes the leader’s objective among all feasible points. During the optimisation procedure, the leader influences the follower by providing a value x fixed as a parametric value. Subsequently, the follower alters the values of the objective and constraints at the leader’s problems by reacting a value y to address the leader’s problem. When the optimal solution (x, y) is obtained, the optimisation process stops. The nested nature of BLPP poses a number of additional challenges compared with conventional single-level optimization problems. In particular: The leader’s behaviour of a bilevel problem may be nonlinear, even if the problems at both levels are linear. It has been proved that BLPP are non-deterministic polynomial hard problem. Theoretically, a leader’s solution is considered valid/feasible only if the corresponding follower’s variables are the true global optimum of the follower problem. Global optimality can only be assured in very limited cases, such as convex and linear problems. However, for most nonlinear and black-box problems, it is not possible to ensure global optimality. In the deceptive case, incorrect follower’s optimal values may cause the objective value to be better than the leader’s true optimal value, which poses a severe challenge to the ranking strategy used in evolutionary optimization technology. Due to each variable of the leader needs to solve the follower programming problem to obtain a feasible solution of the BLPP, the cost of calculation is significantly high. In particularly, if the problem is mathematically not suited to be solved using exact techniques, thus, evolutionary/hybrid techniques are used. Related work With the advancements and developments in human science and technology, various increasingly complex optimisation models have emerged, particularly in the fields of equipment installation, task scheduling, production planning, line optimisation, software design, and tariff setting. Therefore, how these models can be effectively solved has emerged as an important main research topic in the optimisation field. Over the past few decades, many optimisation research results have mainly focused on single-objective and multi-objective optimisation models, whereas relatively few studies have been conducted on hierarchical optimisation models. Leveraging the one-time decision theory to produce multiple short-life cycle products, Zhu and Guo [1]studied BLPP applications with follower problem operators for manufacturers and used classical optimisation methods to solve them. Nasrolahpour et al. [2] developed an energy storage system for merchant pricing based on a two-tier complementary model that can determine the most favourable transaction behaviour. Under a multi-objective framework, Ahmad et al. [3] proposed a simple multi-objective bilevel linear programming method, which considered reservoir managers and multiple water-use departments as a hierarchical structure to optimally allocate limited water resources. More applications of optimisation models can be found in the literature [4–13]. Many practical applications have promoted theoretical research on bilevel programming problems, such as developing efficient algorithms and obtaining optimality conditions. However, owing to the computational complexity of the BLPP itself, it is often very challenging to adopt traditional optimisation methods based on gradient information to address such problems. Currently, only special bilevel programming approachs, such as linear and convex quadratic programs, can obtain the optimal solution to these problem through optimality conditions. For other types of bilevel programming problems research has focused on designing swarm intelligence and hybrid algorithms, which are currently more effective algorithm frameworks and have achieved better solutions when tested against certain problems. Existing methods for solving the BLPP can be grouped into the following categories: A). Classical approaches Some classic methods for BLPP include the simplex [14], branch-and-bound method [15], descending gradient method [16], and penalty function methods [17–19]. Classical methods typically apply the optimality conditions of the follower to convert the BLPP into a single-level problem. Dempe [14]used the simplex method to propose an algorithm for solving linear BLPP. The algorithm introduces slack variables to find a basis as a feasible solution of the BLPP, through a simplex and iterative method. Although this method is more effective for solving small-scale linear BLPP it cannot be directly extended to nonlinear problems. Susanne et al. [16] replaced the the follower problem a with Karush-Kuhn-Tucker condition and applied the optimal value method to transform the original problem into a single-level problems. This approach could effectively satisfy the complementary conditions in the mathematical programming problem in Banach space and introduce M stationarity. Using the dual gap as a penalty function, the literature of White [17] proposed an effective algorithm for solving the BLPP. It adopts a newly designed precise penalty function method to obtain the global optimal solution of the linear situation and conducts a novel theoretical analysis. In the literature [18, 19], a weak linear BLPP dealt with the penalty function and Karush-Kuhn-Tucker conditions. B). Evolutionary approaches The evolutionary algorithms, as representatives of the swarm intelligence optimisation technology, have been widely used to solve various BLPP over the past few decades. An early evolutionary algorithm was proposed by Mathieu et al. [20], who applied linear programming methods to handle followers’ problems and used a genetic algorithm to explore the search space of leaders’ problems. Focused on BLPPs whose follower’s problem is convex programming, Wang [21] proposed an evolutionary algorithm with an embedded constraint processing method. This method applies Karush-Kuhn-Tucker conditions to transform the followers’ problems, turning the original model into a single-level programming model. In addition, to obtain sufficient feasible solutions, a constraint-processing method was designed. Based on similar optimality conditions, Li [22] presented a genetic algorithm that solves nonlinear/linear fractional BLPP. In another study on the presence of the so-called pseudo-feasible solutions in evolutionary bilevel optimisation(QBCA-2), the focus was on determining how pseudo-feasible solutions can affect the performance of an evolutionary algorithm. Moreover, a novel and scalable set of test problems with characterised pseudo-feasible solutions was introduced in [23]. In the literature [24], a co-evolutionary algorithm was proposed to solve BLPP. In the evolution process, the follower problem was solved in two stages. Aboelnaga et al. [25] proposed an improved genetic algorithm and chaotic search method. Goshu et al. [26] proposed a metaheuristic algorithm for random BLPP. An evolutionary algorithm for solving nonlinear bilevel programming problems has been presented in the literature [27]. The algorithm is designed by reflecting the optimal solution of the follower problem to the leader problem. To ensure the quality of each iteration, the algorithm adaptive changes the population size during the evolution process and generates individuals using the tabu search method. C). Hybrid approaches The hybrid algorithm [28] is a common method used to solve BLPP. Abo-Elnaga et al. [29] proposed a multi-sine-cosine algorithm to solve nonlinear BLPP. The sine-cosine algorithms based on three different populations were presented. The first population was used to deal with the leader programming problem, while the second one addressed the follower programming problem. In addition, the Karush-Kuhn-Tucker condition was applied to transform the initial problem into a constrained optimisation problem, which was solved using the third population. If the objective function value is equal to zero, then the solution obtained by solving the leader and follower problems is deemed feasible. Wang [30] proposed a particle swarm distribution estimation algorithm with an embedded Estimation of Distribution Algorithm to solve nonlinear BLPP. Before executing the speed and location update rules in Partical Swarm Optimization, a Gaussian distribution was applied to generate offspring to replace some inferior individuals (particles) in the current population. D). EA based on approximate methods To avoid lengthy calculation processes, some evolutionary algorithms utilise approximation techniques to improve the efficiency of the intermediate calculation process. A bilayer covariance matrix adaptive evolution strategy(BLCMAES) is proposed in [31]. The method designed a sharing mechanism so that prior knowledge of the follower problem can be extracted from the leader optimizer, reduced the number of evaluations of the follower problem. Furthermore, an optimization-based elite retention mechanism is proposed to keep track of the elite and avoid incorrect solutions. Sinha et al. [32] proposed an evolutionary algorithm that uses approximate functions to adress BLPPs. For offspring generated by evolutionary operators, the approximation method can reduce the number of evaluations of the follower objective function. Using approximate Karush-Kuhn-Tucker conditions, Sinha et al. [33, 34] transformed the BLPP into a single-level optimisation problem; then leveraged an evolutionary algorithm embedded with the idea of neighbourhood measurements for the transformed model. Sinha et al. [35] presented an evolutionary optimisation algorithm (BLEAQ-2) that fits an extreme value mapping using the relationship between the leader and follower variables. Islam et al. [36] introduced an evolutionary algorithm that employs three types of surrogate models to approximate the optimal solution of the follower problem. Based on the linear programming optimality conditions, Li [37] proposed a genetic algorithm for solving linear BLPP. In this method, the follower is adopted as the search object of the evolutionary algorithm, and the lower variables in the leader problems are replaced by possible solution functions. This process locally optimises the leader variables. Research motivation In the field of engineering and economic management, optimization models involving different decision-making levels often appear, such problems are called multi-level optimization problems. Since the decision-making process of the problem is hierarchical, it is also called a hierarchical optimization problem. BLPP is a typical representative of multi-level hierarchical optimization problems, and has become an important research field for optimization problems due to its extensive practical application background and algorithmic challenges. Different from multi-objective optimization, the decision makers of BLPP are at two different levels, and this hierarchical structure often leads to problems that are non-convex and non-differentiable. Based on the above characteristics, it is often difficult for traditional optimization algorithms based on gradients to find the global optimal solution to the BLPP. Evolutionary algorithm be increasingly used to solve BLPP because it have the characteristics of global convergence, and there is no restriction on functions that are convex and differentiable. In this paper, driven the correlation coefficient and surrogate model, an evolutionary algorithm (TCEA) is proposed to solve complex nonlinear BLPP. The algorithm has the following characteristics: First, the isodata clustering method [38] is used to group the initial populations, and then determine the correlation coefficients of the leader and follower objective functions in each group based on the rank of the leader (follower) objective function, after that, some points are selected and updated in the offsprings based on the correlation coefficient value in each group. Second, for the offspring individuals produced by the crossover and mutation operators, the surrogate model is used to approximate the solution of the follower programming problem, and so does to reduce the number of evaluations of follower problems. The remainder of this paper is organized as follows. The basic concepts of BLPP are described in the next Section. The correlation coefficients and surrogate models are presented in Section 3. A new evolutionary algorithm based on the above methods are stated in Section 4. Experimental results and comparisons are provided in Section 5. We conclude our approach in Section 6. Basic concepts Some basic definitions for problem (1) are summarized as follows: Constraint region: Follower’s feasible set for x fixed: Projection of S onto the leader’s decision space: Follower’s rational reaction set for each x ∈ S(x) Inducible region: In terms of the aforementioned definitions, problem (1) can also be witten as In order to ensure that problem (1) is well posed, in the remainder, we always assume that (A1) S is nonempty and compact; (A2) For all decisions taken by the leader’s, the follower’s has some room to react, that is, S(x)≠ϕ; (A3) The follower’s problem has unique optimal solution for each fixed x. Main improvement schemes Correlation coefficients The challenge faced in solving BLPP is that evaluating the follower programming problem involves a large amount of calculation. Therefore, to reduce the optimization times of the follower problem, we update a part of the population points by means of the relationship between the objective functions of the leader and follower (called correlation coefficients), thereby reducing the calculation times of the follower optimization problem. The correlation coefficients are defined as follows: Given z points, we acquire the leader and the follower objective function values F(xi, yi) and f(xi, yi), i = 1, 2, …, z. Then, we sort the objective functions based on their values, and the sequence numbers are denoted and , i = 1, 2, …, z, respectively. The sequence number difference between the sorted leader and follower objectives is obtained as follows: (2) Correlation coefficient is defined as follows: (3) The larger the value of ρ is, the more similar the changing trends of the objective functions of the leader and the follower are. Conversely, the changing trends of the objective functions of the leader and the follower are different. Particularly, if ρ = 1, the objective functions of the leader and follower exhibit exactly the same trends. When ρ = −1, the objective functions of the leader and the follower exhibit opposite trends. For example, we take z = 5 and set F(x1, y1) = 2.5, F(x2, y2) = 2.3, F(x3, y3) = 1.4, F(x4, y4) = 1.6, F(x5, y5) = 5.8. After sorting rank: F(x3, y3) = 1.4, F(x4, y4) = 1.6, F(x2, y2) = 2.3, F(x1, y1) = 2.5, F(x5, y5) = 5.8. . And given f(x1, y1) = 5.5, f(x2, y2) = 2.7, f(x3, y3) = 1.6, f(x4, y4) = 3.8, f(x5, y5) = 6.6. After sorting rank: f(x2, y2) = 1.6, f(x3, y3) = 2.7, f(x1, y1) = 3.8, f(x5, y5) = 5.5, f(x4, y4) = 6.6. . Then and the correlation coefficients It turns to be that, the objective functions of the leader and the follower have the more similar the changing trends. Surrogate models In BLPP, the procedure to finding a feasible solution results in a significant amount of computation in solving BLPP, particularly when the problem is large. And the optimal solutions to the follower’s problem are always determined by the leader’s variables. This means that the optimal solution of the follower problem is a function of the leader’s variables. However, the function is often implicit and can not be obtained analytically. In the proposed approach, we take the polynomial fitting as surrogate models [39] to estimate the optimal solutions to the follower’s problems. The polynomial fitting demonstrates better performance in fitting unknown functions and can efficiently decrease the computational times of the follower problems. It is noteworthy that for these fitting points, each follower’s variable value must be optimal when the leader’s components are fixed. In the proposed algorithm, the polynomial fitting is generated as follows. First, an initial population of N points xi, i = 1, 2, …, N is gotten, the optimal solutions to the follower problem are denoted by yi, i = 1, 2, …, N. thus N point pairs of can be obtained. These point pairs are used as fitting nodes to generate an polynomial curve. (4) where (5) i.e. each of yj, j = 1, …, m, is a function of x and y(x) = (y1, y2, …, ym). Where k is the highest degree of the polynomial, a0, a1, a2, …, ak is the undetermined coefficient and calculated by: (6) According to the above-mentioned method, we can obtain the approximate optimal solutions to the follower problem. Correlation coefficients The challenge faced in solving BLPP is that evaluating the follower programming problem involves a large amount of calculation. Therefore, to reduce the optimization times of the follower problem, we update a part of the population points by means of the relationship between the objective functions of the leader and follower (called correlation coefficients), thereby reducing the calculation times of the follower optimization problem. The correlation coefficients are defined as follows: Given z points, we acquire the leader and the follower objective function values F(xi, yi) and f(xi, yi), i = 1, 2, …, z. Then, we sort the objective functions based on their values, and the sequence numbers are denoted and , i = 1, 2, …, z, respectively. The sequence number difference between the sorted leader and follower objectives is obtained as follows: (2) Correlation coefficient is defined as follows: (3) The larger the value of ρ is, the more similar the changing trends of the objective functions of the leader and the follower are. Conversely, the changing trends of the objective functions of the leader and the follower are different. Particularly, if ρ = 1, the objective functions of the leader and follower exhibit exactly the same trends. When ρ = −1, the objective functions of the leader and the follower exhibit opposite trends. For example, we take z = 5 and set F(x1, y1) = 2.5, F(x2, y2) = 2.3, F(x3, y3) = 1.4, F(x4, y4) = 1.6, F(x5, y5) = 5.8. After sorting rank: F(x3, y3) = 1.4, F(x4, y4) = 1.6, F(x2, y2) = 2.3, F(x1, y1) = 2.5, F(x5, y5) = 5.8. . And given f(x1, y1) = 5.5, f(x2, y2) = 2.7, f(x3, y3) = 1.6, f(x4, y4) = 3.8, f(x5, y5) = 6.6. After sorting rank: f(x2, y2) = 1.6, f(x3, y3) = 2.7, f(x1, y1) = 3.8, f(x5, y5) = 5.5, f(x4, y4) = 6.6. . Then and the correlation coefficients It turns to be that, the objective functions of the leader and the follower have the more similar the changing trends. Surrogate models In BLPP, the procedure to finding a feasible solution results in a significant amount of computation in solving BLPP, particularly when the problem is large. And the optimal solutions to the follower’s problem are always determined by the leader’s variables. This means that the optimal solution of the follower problem is a function of the leader’s variables. However, the function is often implicit and can not be obtained analytically. In the proposed approach, we take the polynomial fitting as surrogate models [39] to estimate the optimal solutions to the follower’s problems. The polynomial fitting demonstrates better performance in fitting unknown functions and can efficiently decrease the computational times of the follower problems. It is noteworthy that for these fitting points, each follower’s variable value must be optimal when the leader’s components are fixed. In the proposed algorithm, the polynomial fitting is generated as follows. First, an initial population of N points xi, i = 1, 2, …, N is gotten, the optimal solutions to the follower problem are denoted by yi, i = 1, 2, …, N. thus N point pairs of can be obtained. These point pairs are used as fitting nodes to generate an polynomial curve. (4) where (5) i.e. each of yj, j = 1, …, m, is a function of x and y(x) = (y1, y2, …, ym). Where k is the highest degree of the polynomial, a0, a1, a2, …, ak is the undetermined coefficient and calculated by: (6) According to the above-mentioned method, we can obtain the approximate optimal solutions to the follower problem. Proposed algorithm In this manuscript, an evolutionary algorithm based on surrogate models and correlation coefficients, denoted by TCEA, is developed to solve BLPP. Fig 1 gives the flowchart of TCEA. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. The frame diagram of TCEA. https://doi.org/10.1371/journal.pone.0273564.g001 The detailed procedure of the proposed algorithm can be described as follows: Step 1 (Initial population) The idea of uniform design [40] is adopted to produces N points xi, i = 1, …, N, resulting in an initial population pop(0) of size N. Set gen = 0, D = Φ. Step 2 (Fitness evaluation) For each xi, solve the follower problems and obtain the optimal solution yi, i = 1, …, N. These points are put into D. The value of the leader objectives are taken as F(xi, yi), i = 1, …, N. Construct the polynomial fitting (surrogate models), as in Section 3.2. Use the isodata method to divide the generated points into p groups, denoted as I1, I2, …, Ip. Then take advantage of the correlation coefficients method in Section 3.1 to acquire the value of ρ in each group, denoted by ρ1, ρ2, …, ρp. Step 3 (Crossover) For each crossover parent individual xi, take an best individual as , and perform the following crossover operator using the spherical search method: Set: here α ∈ (0, 1) is a constant, which is called the shrinkage rate of the radius, and the method is as follows: (7) ui, li is the leader and follower’s bounds of xi. Take a uniformly distributed value of θ1, θ2, …, θJ ∈ [0, 2π],β1, β2, …, βJ ∈ (−π/2, π/2) then use the spherical search method to generate a crossover operator as follows: (n is the dimension of the leader variable) (8) (9) Step 4 (Mutation) Gaussian mutation is adopted. Suppose that is an individual chosen for mutation, then the offspring of is generated as follows: (10) Step 5 (Offspring population pop′(gen)) For offspring set (xo1, xo2, …, xoλ) generated through the crossover and mutation operation. A surrogate model, the polynomial fitting, is used to obtain the approximated solutions to the follower’s program. We only update some of offspring points based on the value of ρ1, ρ2, …, ρp in each group as follows: Case 1: If in the group τ, τ = 1, 2, …, p, the value of ρτ is greater than the given threshold μ > 0 (according to the results of the experiment), that is to say, this value is near to 1, it means that the leader’s and the follower’s objectives have the same changing trend. If the leader objective function at point (xoi, yoi) satisfies (a predetermined threshold), it means that the leader’s objective maybe become minimizer when the follower’s objective is optimized. As a result, point (xoi, yoi) needs to be updated, that is to say, the follower problem is solved and the solution to the follower’s problem is updated. If it is at point (xoi, yoi), it means that the point is unpromising even if the follower’s solution is updated. As a result, point (xoi, yoi) is not updated. Case 2: If in the group τ, τ = 1, 2, …, p, the value of ρτ is less than the given threshold −μ (according to the results of the experiment), i.e. The value is near to −1, it means that the leader and the follower objective functions have opposite trends. At this case, the leader’s objective maybe become worse when the follower’s objective is minimized. However, it is expected the worsen objective values are still better than the predetermined threshold F*(a bilevel feasible objective value). This means these points with small objective values are potential to be refined. The smallest objective is denoted by Fbest(maybe infeasible), and these points with objective value in [Fbest, F*] should be further updated in a probabilistic sense. At point (xoi, yoi), if the objective satisfies that F(xoi, yoi)