Bilevel Optimization based on Iterative Approximation of Multiple Mappings

TL;DR

This paper introduces an evolutionary optimization method for bilevel problems that reduces computational costs by iteratively approximating key mappings, showing significant performance improvements over existing methods.

Contribution

It presents a novel combined theory-based and population-based approach using mappings for bilevel optimization, which has not been previously proposed.

Findings

Significant performance gains over existing algorithms

Effective approximation of lower level mappings

Validated on numerous test problems

Abstract

A large number of application problems involve two levels of optimization, where one optimization task is nested inside the other. These problems are known as bilevel optimization problems and have been studied by both classical optimization community and evolutionary optimization community. Most of the solution procedures proposed until now are either computationally very expensive or applicable to only small classes of bilevel optimization problems adhering to mathematically simplifying assumptions. In this paper, we propose an evolutionary optimization method that tries to reduce the computational expense by iteratively approximating two important mappings in bilevel optimization; namely, the lower level rational reaction mapping and the lower level optimal value function mapping. The algorithm has been tested on a large number of test problems and comparisons have been performed with other algorithms. The results show the performance gain to be quite significant. To the best knowledge of the authors, a combined theory-based and population-based solution procedure utilizing mappings has not been suggested yet for bilevel problems.

Tables

Table 1. TABLE I: Median function evaluations required at upper level (UL) and the lower level (LL) from 31 runs of Ψ Ψ \Psi -approximation algorithm, φ 𝜑 \varphi -approximation algorithm and nested algorithm.

	UL Func. Evals.			LL Func. Evals.			Savings
	$\phi$ Appx	$\Psi$ Appx	Nested	$\phi$ Appx	$\Psi$ Appx	Nested	$\frac{\text{Nested}-\phi \text{ Appx}}{\text{Nested}}$	$\frac{\text{Nested}-\Psi \text{ Appx}}{\text{Nested}}$
	Med	Med	Med	Med	Med	Med
TP1	137	138	-	1539	1948	-	Large	Large
TP2	158	187	444	1614	2820	5252	69%	47%
TP3	198	132	642	2710	1461	6530	59%	78%
TP4	309	419	1760	2976	6347	18073	83%	66%
TP5	165	244	633	2571	2763	6616	62%	59%
TP6	115	91	153	1512	1170	2169	30%	46%
TP7	169	130	196	2386	1426	2606	9%	44%
TP8	201	328	423	2443	4722	8233	69%	42%

Table 2. TABLE II: Mean of total function evaluations (UL evaluations +LL evaluations) required by different approaches.

	Mean Func. Evals. (UL+LL)
	$\phi$ Appx.	$\Psi$ Appx.	Nested	WJL [50]	WLD [51]
TP1	1,611	2,421	34,462	85,499	86,067
TP2	1,923	3,262	6,235	256,227	171,346
TP3	2,624	1,482	8,125	92,526	95,851
TP4	3,612	6,721	19,948	291,817	211,937
TP5	2,812	3,388	7,398	77,302	69,471
TP6	1,578	1,034	1,534	163,701	65,942
TP7	2,110	1,456	2,286	1,074,742	944,105
TP8	2,734	4,434	5,325	213,522	182,121

Table 3. TABLE III: Statistics for upper level function evaluations for φ 𝜑 \varphi -approximation algorithm on the modified test problems (m-TP). The φ 𝜑 \varphi -approximation algorithm and nested algorithm fail on all the modified test problems.

	$\phi$ Appx.			$\Psi$ Appx.	Nested
	Min	Med	Max	Min/Med/Max	Min/Med/Max
m-TP1	138	192	344	-	-
m-TP2	124	236	-	-	-
m-TP3	140	242	699	-	-
m-TP4	185	545	2,582	-	-
m-TP5	172	242	977	-	-
m-TP6	159	181	559	-	-
m-TP7	119	227	501	-	-
m-TP8	158	462	2,119	-	-

Table 4. TABLE IV: Statistics for lower level function evaluations from 31 runs of the φ 𝜑 \varphi -approximation algorithm on the modified test problems (m-TP). The φ 𝜑 \varphi -approximation algorithm and nested algorithm fail on all the modified test problems.

	$\phi$ Appx.			$\Psi$ Appx.	Nested
	Min	Med	Max	Min/Med/Max	Min/Med/Max
m-TP1	1,988	2,477	8,334	-	-
m-TP2	2,394	4,420	-	-	-
m-TP3	1,404	3,321	12,353	-	-
m-TP4	1,911	5,632	25,356	-	-
m-TP5	3,129	4,166	15,345	-	-
m-TP6	2,498	3,464	9,325	-	-
m-TP7	1,476	5,635	12,256	-	-
m-TP8	2,721	6,324	28,993	-	-

Table 5. TABLE V: Step-by-step procedure for BLEAQ-II

Step	Description
1	Initialization:
	Generate an initial upper level population $x_u^{\left(1\right)},\dots ,x_u^{\left(N\right)}$ randomly or by a problem-specific method (see Section VI-A).
	(a) For each $x_u^{\left(j\right)}$, find a corresponding optimal lower level solution $x_l^{\left(j\right)}\in \Psi \left(x_u^{\left(j\right)}\right)$ by solving the lower level problem. Set $𝒫=\{(x_u^{\left(j\right)},x_l^{\left(j\right)}),j=1,\dots ,N\}$ (see Section VI-E).
	(b) Tag all vectors $(x_u^{\left(j\right)},x_l^{\left(j\right)})\in 𝒫$ for which a lower level optimization has been successfully performed as 1 and store them in the archive $𝒜$.
	(c) Assign fitness to all the members based on upper level function and constraints (refer to Section VI-B).
2	Reproduction:
	(a) Parent selection: Randomly choose $2\mu$ members from the population $𝒫$, and perform a tournament selection based on the upper level fitness. This produces $\mu$ parents, denoted by ${𝒫}_{\text{par}}$.
	(b) Offspring generation: Create $\lambda$ offsprings, denoted by ${𝒫}_{\text{off}}$, from the set of parents ${𝒫}_{\text{par}}$ using genetic operators (refer to Section VI-C).
3	Offspring Update:
	For each offspring $x^{\left(j\right)}=(x_u^{\left(j\right)},x_l^{\left(j\right)})\in {𝒫}_{\text{off}}$ produced in the previous step, update the lower level decision $x_l^{\left(j\right)}$ using the following strategy:
	(a) Optimization: If the number of Tag 1 members in $𝒫$ is less than half of the size of $𝒫$, perform lower level optimization to ensure that $x_l^{\left(j\right)}\in \Psi \left(x_u^{\left(j\right)}\right)$ (as described in Step 1.(b)). If the lower level optimization is successful, tag the offsprings as 1 and add it to the archive $𝒜$.
	(b) Approximations: If the number of Tag 1 members in $𝒫$ is more than half of the size of $𝒫$, for each offspring $x^{\left(j\right)}=(x_u^{\left(j\right)},x_l^{\left(j\right)})\in {𝒫}_{\text{off}}$, use its neighboring members in the archive $𝒜$ to construct a local quadratic approximation for $\Psi$-mapping ($q_{\Psi }$) as well as $\phi$-mapping ($q_{\phi }$). Compare the mean squared error of the approximations ($e_{mse}^{\Psi }$, $e_{mse}^{\phi }$). If $e_{mse}^{\Psi }\le e_{mse}^{\phi }$ then update the lower level decision associated with the upper level $x_u^{\left(j\right)}$ by setting $x_l^{\left(j\right)}=q_{\Psi }\left(x_u^{\left(j\right)}\right)$; otherwise solve the auxiliary optimization problem in Section VI-F by fixing $x_u^{\left(j\right)}$ and varying $x_l^{\left(j\right)}$; the optimal $x_l^{\left(j\right)}$ is paired with $x_u^{\left(j\right)}$ to form the offspring.
4	Improvements:
	Identify the Tag 1 member in the current generation in $𝒫$ with the best fitness, denoted as $x_{best}^{\left(j\right)}$. Perform a local search in the vicinity of $x_{best}^{\left(j\right)}$ after every $k$ generations and update $x_{best}^{\left(j\right)}$ if there is an improvement.
	(a) Local search: Construct local quadratic approximations of both $\Psi$-mapping and $\phi$-mapping using members in the vicinity of $x_{best}^{\left(j\right)}$ in the archive $𝒜$ and record the mean squared error of the approximations ($e_{mse}^{\Psi }$, $e_{mse}^{\phi }$). Apply local search in $x_{best}^{\left(j\right)}$ vicinity using one of the two single level reduction methods described in Sections III-A and III-B (refer to Section VI-G).
5	Termination check:
	Perform a termination check. If false, proceed to the next generation (Step 2).

Table 6. TABLE VI: The lower level evolutionary algorithm is described below that takes an upper level member as input and solves the corresponding lower level problem.

Step	Description
1	Initialization: Generate an initial lower level population $x_l^{\left(1\right)},\dots ,x_l^{\left(N\right)}$ randomly and assign fitness using lower level objective and constraints.
2	Genetic Operations: Randomly choose $2\mu$ members from the population, and perform a tournament selection leading to $\mu$ parents. Create $\lambda$ offspring using the genetic operations described in Section VI-C. Assign fitness to each offspring.
3	Update: Choose $r$ members randomly from the population and pool them with $\lambda$ offspring. Sort the pool by fitness and replace the $r$ members from the population by the best $r$ members from the pool.
4	Termination check: Perform a termination check as described in Section VI-D. If false, proceed to the next generation (Step 2).

Table 7. TABLE VII: Median function evaluations on TP test suite

	UL Func. Evals.			LL Func. Evals.			BLEAQ-II Savings
	BLEAQ-II	BLEAQ	Nested	BLEAQ-II	BLEAQ	Nested	$\frac{\text{BLEAQ}-\text{BLEAQ-II}}{\text{BLEAQ}}$	$\frac{\text{Nested}-\text{BLEAQ-II}}{\text{Nested}}$
	Med	Med	Med	Med	Med	Med
TP1	136	155	-	242	867	-	63%	98%
TP2	255	185	436	440	971	5,686	40%	63%
TP3	158	155	633	224	894	6,867	64%	98%
TP4	198	357	1,755	788	1,772	19,764	54%	98%
TP5	272	243	576	967	1,108	6,558	8%	84%
TP6	161	155	144	323	687	1,984	43%	97%
TP7	112	255	193	287	987	2,870	68%	95%
TP8	241	189	403	467	913	7,996	36%	61%

Table 8. TABLE VIII: Median function evaluations on low dimension SMD test suite

	UL Func. Evals.			LL Func. Evals.			BLEAQ-II Savings
	BLEAQ-II	BLEAQ	Nested	BLEAQ-II	BLEAQ	Nested	$\frac{\text{BLEAQ}-\text{BLEAQ-II}}{\text{BLEAQ}}$	$\frac{\text{Nested}-\text{BLEAQ-II}}{\text{Nested}}$
	Med	Med	Med	Med	Med	Med
SMD1	123	98	164	8,462	13,425	104,575	37%	92%
SMD2	114	88	106	7,264	11,271	74,678	35%	90%
SMD3	264	91	136	12,452	15,197	101,044	17%	87%
SMD4	272	110	74	8,600	12,469	59,208	29%	85%
SMD5	126	80	93	14,490	19,081	73,500	24%	80%
SMD6	259	-	116	914	-	3,074	Large	63%
SMD7	180	98	67	8,242	12,580	56,056	34%	85%
SMD8	644	228	274	22,866	35,835	175,686	35%	87%
SMD9	201	125	127	10,964	16,672	101,382	34%	89%
SMD10	780	431	-	19,335	43,720	-	54%	Large
SMD11	1735	258	260	134,916	158,854	148,520	14%	8%
SMD12	203	557	-	25,388	135,737	-	81%	Large
SMD13	317	126	211	13,729	17,752	138,089	21%	90%
SMD14	1,014	-	168	12,364	-	91,197	Large	85%

Table 9. TABLE IX: Median function evaluations on high dimension SMD test suite

	UL Func. Evals.			LL Func. Evals.			BLEAQ-II Savings
	BLEAQ-II	BLEAQ	Nested	BLEAQ-II	BLEAQ	Nested	$\frac{\text{BLEAQ}-\text{BLEAQ-II}}{\text{BLEAQ}}$	$\frac{\text{Nested}-\text{BLEAQ-II}}{\text{Nested}}$
	Med	Med	Med	Med	Med	Med
SMD1	670	370	760	52,866	61,732	1,776,426	14%	97%
SMD2	510	363	652	44,219	57,074	1,478,530	22%	97%
SMD3	1369	630	820	68,395	90,390	1,255,015	23%	94%
SMD4	580	461	765	35,722	59,134	1,028,802	39%	96%
SMD5	534	464	645	65,873	92,716	1,841,569	29%	96%
SMD6	584	-	824	3,950	-	156,2003	Large	99%
SMD7	1,486	-	-	83,221	-	-	Large	Large
SMD8	6,551	-	-	231,040	-	-	Large	Large

Table 10. TABLE X: Standard test problems TP1-TP5. (Note that x = x u 𝑥 subscript 𝑥 𝑢 x=x_{u} and y = x l 𝑦 subscript 𝑥 𝑙 y=x_{l} )

Problem	Formulation	Best Known Sol.
TP1
$n=2$, $m=2$	$\begin{array}{c}\underset{(x,y)}{\text{Minimize}}F(x,y)={\left(x_1-30\right)}^2+{\left(x_2-20\right)}^2-20y_1+20y_2,\\ \text{s.t.}\\ y\in \underset{\left(y\right)}{argmin}\left\{\begin{array}{c}f(x,y)={\left(x_1-y_1\right)}^2+{\left(x_2-y_2\right)}^2\\ 0\le y_i\le 10,i=1,2\end{array}\right\},\\ x_1+2x_2\ge 30,x_1+x_2\le 25,x_2\le 15\end{array}$	$\begin{array}{c}\\ \\ \\ \\ \\ \\ \\ F=225.0\\ f=100.0\end{array}$
TP2
$n=2$, $m=2$	$\begin{array}{c}\underset{(x,y)}{\text{Minimize}}F(x,y)=2x_1+2x_2-3y_1-3y_2-60,\\ \text{s.t.}\\ y\in \underset{\left(y\right)}{argmin}\left\{\begin{array}{c}f(x,y)={\left(y_1-x_1+20\right)}^2+{\left(y_2-x_2+20\right)}^2\\ x_1-2y_1\ge 10,x_2-2y_2\ge 10\\ -10\ge y_i\ge 20,i=1,2\end{array}\right\},\\ x_1+x_2+y_1-2y_2\le 40,\\ 0\le x_i\le 50,i=1,2.\end{array}$	$\begin{array}{c}\\ \\ \\ \\ \\ \\ \\ F=0.0\\ f=100.0\end{array}$
TP3
$n=2$, $m=2$	$\begin{array}{c}\underset{(x,y)}{\text{Minimize}}F(x,y)=-{\left(x_1\right)}^2-3{\left(x_2\right)}^2-4y_1+{\left(y_2\right)}^2,\\ \text{s.t.}\\ y\in \underset{\left(y\right)}{argmin}\left\{\begin{array}{c}f(x,y)=2{\left(x_1\right)}^2+{\left(y_1\right)}^2-5y_2\\ {\left(x_1\right)}^2-2x_1+{\left(x_2\right)}^2-2y_1+y_2\ge -3\\ x_2+3y_1-4y_2\ge 4\\ 0\le y_i,i=1,2\end{array}\right\},\\ {\left(x_1\right)}^2+2x_2\le 4,\\ 0\le x_i,i=1,2\end{array}$	$\begin{array}{c}\\ \\ \\ \\ \\ \\ \\ F=-18.6787\\ f=-1.0156\end{array}$
TP4
$n=2$, $m=3$	$\begin{array}{c}\underset{(x,y)}{\text{Minimize}}F(x,y)=-8x_1-4x_2+4y_1-40y_2-4y_3,\\ \text{s.t.}\\ y\in \underset{\left(y\right)}{argmin}\left\{\begin{array}{c}f(x,y)=x_1+2x_2+y_1+y_2+2y_3\\ y_2+y_3-y_1\le 1\\ 2x_1-y_1+2y_2-0.5y_3\le 1\\ 2x_2+2y_1-y_2-0.5y_3\le 1\\ 0\le y_i,i=1,2,3\end{array}\right\},\\ 0\le x_i,i=1,2\end{array}$	$\begin{array}{c}\\ \\ \\ \\ \\ \\ \\ F=-29.2\\ f=3.2\end{array}$
TP5
$n=2$, $m=2$		$\begin{array}{c}\\ \\ \\ \\ \\ \\ \\ F=-3.6\\ f=-2.0\end{array}$

Table 11. TABLE XI: Standard test problems TP6-TP8. (Note that x = x u 𝑥 subscript 𝑥 𝑢 x=x_{u} and y = x l 𝑦 subscript 𝑥 𝑙 y=x_{l} )

Problem	Formulation	Best Known Sol.
TP6
$n=1$, $m=2$	$\begin{array}{c}\underset{(x,y)}{\text{Minimize}}F(x,y)={\left(x_1-1\right)}^2+2y_1-2x_1,\\ \text{s.t.}\\ y\in \underset{\left(y\right)}{argmin}\left\{\begin{array}{c}f(x,y)={\left(2y_1-4\right)}^2+\\ {\left(2y_2-1\right)}^2+x_1{y}_1\\ 4x_1+5y_1+4y_2\le 12\\ 4y_2-4x_1-5y_1\le -4\\ 4x_1-4y_1+5y_2\le 4\\ 4y_1-4x_1+5y_2\le 4\\ 0\le y_i,i=1,2\end{array}\right\},\\ 0\le x_1\end{array}$	$\begin{array}{c}\\ \\ \\ \\ \\ \\ \\ F=-1.2091\\ f=7.6145\end{array}$
TP7
$n=2$, $m=2$	$\begin{array}{c}\underset{(x,y)}{\text{Minimize}}F(x,y)=-\frac{\left(x_1+y_1\right)\left(x_2+y_2\right)}{1+x_1{y}_1+x_2{y}_2},\\ \text{s.t.}\\ y\in \underset{\left(y\right)}{argmin}\left\{\begin{array}{c}f(x,y)=\frac{\left(x_1+y_1\right)\left(x_2+y_2\right)}{1+x_1{y}_1+x_2{y}_2}\\ 0\le y_i\le x_i,i=1,2\end{array}\right\},\\ {\left(x_1\right)}^2+{\left(x_2\right)}^2\le 100\\ x_1-x_2\le 0\\ 0\le x_i,i=1,2\end{array}$	$\begin{array}{c}\\ \\ \\ \\ \\ \\ \\ F=-1.96\\ f=1.96\end{array}$
TP8
$n=2$, $m=2$	$\begin{array}{c}\underset{(x,y)}{\text{Minimize}}F(x,y)=\left|2x_1+2x_2-3y_1-3y_2-60\right|,\\ \text{s.t.}\\ y\in \underset{\left(y\right)}{argmin}\left\{\begin{array}{c}f(x,y)={\left(y_1-x_1+20\right)}^2+\\ {\left(y_2-x_2+20\right)}^2\\ 2y_1-x_1+10\le 0\\ 2y_2-x_2+10\le 0\\ -10\le y_i\le 20,i=1,2\end{array}\right\},\\ x_1+x_2+y_1-2y_2\le 40\\ 0\le x_i\le 50,i=1,2\end{array}$	$\begin{array}{c}\\ \\ \\ \\ \\ \\ \\ F=0.0\\ f=100.0\end{array}$

Table 12. TABLE XII: SMD Test Problems. (Note that ( x u ​ 1 , x u ​ 2 ) = ( a , b ) subscript 𝑥 𝑢 1 subscript 𝑥 𝑢 2 𝑎 𝑏 (x_{u1},x_{u2})=(a,b) and ( x l ​ 1 , x l ​ 2 ) = ( c , d ) subscript 𝑥 𝑙 1 subscript 𝑥 𝑙 2 𝑐 𝑑 (x_{l1},x_{l2})=(c,d) )

Problem	Formulation	Solution
SMD13
	$\begin{array}{c}{F}_1={\left(a_1-1\right)}^2+{\sum }_{i=1}^{p-1}\left({\left(a_i-1\right)}^2+{\left(a_{i+1}-{\left(a_i\right)}^2\right)}^2\right),\\ F_2=-{\sum }_{i=1}^q{\sum }_{j=1}^i{\left(c_j\right)}^2,\\ F_3={\sum }_{i=1}^r{\sum }_{j=1}^i{\left(b_j\right)}^2-{\sum }_{i=1}^r{\left(b_i-\log d_i\right)}^2,\\ f_1={\sum }_{i=1}^p\left(\left|a_i\right|+2\left|\sin \left(a_i\right)\right|\right),\\ f_2={\sum }_{i=1}^q{\sum }_{j=1}^i{\left(c_j\right)}^2,\\ f_3={\sum }_{i=1}^r{\left(b_i-\log d_i\right)}^2,\\ a_i\in [-5,10],\forall i\in \{1,2,\dots ,p\},\\ b_i\in [-5,e],\forall i\in \{1,2,\dots ,r\},\\ c_i\in [-5,10],\forall i\in \{1,2,\dots ,q\},\\ d_i\in (0,10],\forall i\in \{1,2,\dots ,r\}.\end{array}$	$\begin{array}{c}{a}_i=1\forall i,\\ b_i=0\forall i,\\ c_i=0\forall i,\\ d_i=1\forall i.\end{array}$
SMD14
	$\begin{array}{c}{F}_1={\left(a_1-1\right)}^2+{\sum }_{i=1}^{p-1}\left({\left(a_i-1\right)}^2+{\left(a_{i+1}-{\left(a_i\right)}^2\right)}^2\right),\\ F_2=-{\sum }_{i=1}^q{\left|c_i\right|}^{i+1}+{\sum }_{i=q+1}^{q+s}{\left(c_i\right)}^2,\\ F_3={\sum }_{i=1}^{r}i{\left(b_i\right)}^2-{\sum }_{i=1}^r\left|d_i\right|,\\ f_1={\sum }_{i=1}^p\lfloor a_i\rfloor ,\\ f_2={\sum }_{i=1}^q{\left|c_i\right|}^{i+1}+{\sum }_{i=q+1,i=i+2}^{q+s-1}{\left(c_{i+1}-c_i\right)}^2,\\ f_3={\sum }_{i=1}^r\left|{\left(b_i\right)}^2-{\left(d_i\right)}^2\right|,\\ a_i\in [-5,10],\forall i\in \{1,2,\dots ,p\},\\ b_i\in [-5,10],\forall i\in \{1,2,\dots ,r\},\\ c_i\in [-5,10],\forall i\in \{1,2,\dots ,q+s\},\\ d_i\in [-5,10],\forall i\in \{1,2,\dots ,r\}.\end{array}$	$\begin{array}{c}{a}_i=1\forall i,\\ b_i=0\forall i,\\ c_i=0\forall i,\\ d_i=0\forall i.\end{array}$

Equations

$$\underset{x_u\in X_U,x_l\in X_L}{\text{mbox}‘‘\min \text{mbox}{}^{\prime \prime }}$$

$$\mathrm{subject}\mathrm{to}$$

$$x_l\in \underset{x_l\in X_L}{\mathrm{argmin}}\{f(x_u,x_l):g_j(x_u,x_l)\le 0,j=1,\dots ,J\}$$

$$G_k(x_u,x_l)\le 0,k=1,\dots ,K$$

$$\Psi (x_u)=$$

$$\underset{x_u\in X_U,x_l\in X_L}{\min }$$

$$\mathrm{subject}\mathrm{to}$$

$$x_l\in \Psi (x_u)$$

$$G_k(x_u,x_l)\le 0,k=1,\dots ,K$$

$$\phi (x_u)=\underset{x_l\in X_L}{\min }\{f(x_u,x_l):g_j(x_u,x_l)\le 0,j=1,\dots ,J\},$$

$$\underset{x_u\in X_U,x_l\in X_L}{\min }$$

$$\mathrm{subject}\mathrm{to}$$

$$f(x_u,x_l)\le \phi (x_u)$$

$$g_j(x_u,x_l)\le 0,j=1,\dots ,J$$

$$G_k(x_u,x_l)\le 0,k=1,\dots ,K.$$

$$\hat{\psi }=\underset{h\in \mathcal{H}}{\mathrm{argmin}}\sum _{i\in \mathcal{I}}L(h(x_u^{(i)}),{\bar{x}}_l^{(i)}),$$

$$L(h(x_u^{(i)}),{\bar{x}}_l^{(i)})=|{\bar{x}}_l^{(i)}-h(x_u^{(i)}){|}^2.$$

$$\hat{\phi }=\underset{u\in \mathcal{H}}{\mathrm{argmin}}\sum _{i\in \mathcal{I}}L(u(x_u^{(i)}),{\bar{f}}^{(i)}),$$

$$L(u(x_u^{(i)}),f^{(i)})=|{\bar{f}}^{(i)}-u(x_u^{(i)}){|}^2.$$

$$F^{new}(x_u,x_l)$$

$$f^{new}(x_u,x_l)$$

$$y_p,y_q$$

$$\underset{x_u\in X_U,x_l\in X_L}{\min }$$

$$\mathrm{subject}\mathrm{to}$$

$$G_k(x_u,x_l)\le 0,k=1,\dots ,K,$$

$$g_j(x_u,x_l)\le 0,j=1,\dots ,J.$$

$$c=z^{(p)}+{\omega }_{\xi }d+{\omega }_{\eta }\frac{p^{(2)}-p^{(1)}}{2}$$

$$\begin{array}[]{l}\alpha_{u}^{T}=\frac{\sum_{i=1}^{n}\sigma^{2}(x_{i})|_{T}}{\sum_{i=1}^{n}\sigma^{2}(x_{i})|_{0}},\end{array}$$

$$\underset{x_l\in X_L}{\min }$$

$$\mathrm{subject}\mathrm{to}$$

$$\hat{f}(x_u,x_l)\le \hat{\phi }(x_u)$$

$${\hat{g}}_j(x_u,x_l)\le 0,j=1,\dots ,J$$

$${\hat{G}}_k(x_u,x_l)\le 0,k=1,\dots ,K.$$

$$\begin{array}[]{l}F(x_{u},x_{l})=F_{1}(x_{u1})+F_{2}(x_{l1})+F_{3}(x_{u2},x_{l2})\\ f(x_{u},x_{l})=f_{1}(x_{u1},x_{u2})+f_{2}(x_{l1})+f_{3}(x_{u2},x_{l2})\\ \mbox{where}\\ \quad\quad x_{u}=(x_{u1},x_{u2})\quad\mbox{and}\quad x_{l}=(x_{l1},x_{l2})\end{array}$$

Full text

I Introduction

Interest in bilevel optimization has been growing due to a number of new applications that are arising in different fields of science and engineering. Bilevel programming is quite common in the area of defense where these problems are studied as attacker-defender problems. The problem was introduced by Bracken and McGill [9] in the area of mathematical programming, where an inner optimization problem acts as a constraint to an outer optimization problem. One of the follow-up papers by Bracken and McGill [10] highlighted the applications of bilevel programming in defense. Since then a number of studies on homeland security [12, 52, 3] have been performed, where it is common to have bilevel, trilevel and even multilevel optimization models. In the area of operations research, bilevel optimization is gaining importance in the context of interdiction and protection of hub-and-spoke networks [26], as most of the critical infrastructures like transportation and communications are predominantly hub-and-spoke. In other game theoretic settings, bilevel optimization has been used in transportation [33, 15, 11], optimal tax policies [25, 42, 39], investigation of strategic behavior in deregulated markets [22], model production processes [35] and optimization of retail channel structures [55]. The applications extend to a variety of other domains, like, facility location [24, 47, 46], chemical engineering [45, 14], structural optimization [8, 13], and optimal control [34, 2] problems. While new applications that are inherently bilevel in nature are arising at a fast pace, the development of computationally efficient algorithms for such problems has not kept the pace with the applications.

A significant body of literature exists on bilevel optimization and its optimality conditions [30, 18, 17, 19, 54] in the classical optimization literature. However, on the algorithm front most attention has been given to only simple instances of bilevel optimization where the objective functions and constraints are linear [53, 7], quadratic [6, 20, 1] or convex [31]. This is not surprising given the fact that bilevel optimization is difficult to an extent that merely evaluating the bilevel optimality of a given solution is an NP-hard task [48]. Researchers have also attempted to solve these problems using computational techniques like evolutionary algorithms. Most of the bilevel algorithms relying on evolutionary framework have been nested in nature [32, 57, 28, 58]. One of the drawbacks of such an approach is that it might be able to solve small instances of bilevel problems, but as soon as the problem scales-up beyond a few variables, the computational requirements increase tremendously. However, the evolutionary algorithms still have a niche in solving these problems as it maintains a population at each iteration of the algorithm. A population of points may allow modeling various mappings in bilevel optimization to reduce the computational expense [40]. Some studies in this direction are [41, 44, 36, 37]. We believe that exploiting some of the mathematical properties of bilevel problems through modeling of various mappings in bilevel is the way forward in solving such problems.

In this paper, we focus on two important mappings in bilevel optimization borrowed from the mathematical optimization literature. The first mapping is the lower level reaction set mapping (known as $\Psi$-mapping), which provides the lower level optimal solution(s) corresponding to any given upper level vector. Considering the upper level problem as the leader’s problem and the lower level problem as the follower’s problem, the reaction set mapping represents the rational decisions of the follower corresponding to any decision taken by the leader. The second mapping is the lower level value function mapping (known as $\varphi$-mapping) that provides the optimal objective function value to the follower’s problem for any given leader’s decision. While the first mapping can be a set-valued mapping, the second mapping is always single-valued. We work with meta-modeling techniques that try to approximate these two mappings and develop a computationally efficient evolutionary algorithm for solving bilevel problems. The algorithm has been tested on a number of test problems, and the computational gain when compared with other techniques is found to be significant. In this paper, we also extend an existing test-suite of bilevel test problems [38] with a couple of additional problems to better evaluate our proposed solution procedure.

The paper is organized as follows. To begin with, we provide a brief literature survey of bilevel optimization using evolutionary algorithms. This is followed by various formulations of the bilevel optimization problem and discussion of the two mappings that we approximate in this paper. Thereafter, we provide the bilevel evolutionary optimization algorithm which is an extension of the algorithm proposed in the previous studies [44, 36, 37]. Following this, we provide the empirical results on a number of test problems. A comparative study with other approaches is also included. Finally, we end the paper with the conclusions section.

II A Survey on Evolutionary Bilevel Optimization

Most of the evolutionary algorithms for bilevel optimization are nested in nature, where one optimization algorithm is used within the other. The outer algorithm handles the upper level problem and the inner algorithm handles the lower level problem. Such a structure necessitates that the inner algorithm is called for every upper level point generated by the outer algorithm. Therefore, nested approaches can be quite computationally demanding, and can only be applied to small scale problems. One can find studies with evolutionary algorithm being used for the upper level problem and classical approach being used for the lower level problem. If the lower level problem is complex, researchers have used evolutionary algorithms at both levels. Below we provide a review of evolutionary bilevel optimization algorithms from the past.

Mathieu et al. [32] was one of the first to propose a bilevel algorithm using evolutionary algorithms. He used a genetic algorithm to handle the upper level problem and linear programming to solve the lower level problem for every upper level member generated using genetic operations. This study was followed by nesting the Frank-Wolfe algorithm (reduced gradient method) within a genetic algorithm in Yin [57]. Other authors utilized similar nested schemes in [29, 28, 58]. Studies involving evolutionary algorithms at both levels include [4, 5], where authors have used differential evolution at both levels in the first study, and differential evolution within ant colony optimization in the second study.

Replacing the lower level problem in bilevel optimization with its KKT conditions is a common approach for solving the problem both in classical as well as evolutionary computation literature. Some of the evolutionary studies that utilize this idea include [21, 50]. The approach has been popular and even recently researchers are relying on reducing the bilevel problem into single level problem using KKT and solving the reduced problem using evolutionary algorithm, for example, see [51, 23, 27, 49].

While KKT conditions can only be applied to problems where the lower level adheres to certain mathematically simplifying assumptions, the researchers are exploring techniques that can solve more general instances of bilevel optimization problems. Some of the approaches are based on meta-modeling the mappings within bilevel optimization, while others may be based on meta-modeling the entire bilevel problem itself. Studies in this direction include [44, 36, 37]. In this paper, we aim to develop an algorithm that tries to capture two important mappings in bilevel optimization; namely, the lower level reaction set mapping and the lower level value function mapping, in order to reduce the computational complexity of the problem.

III Different Bilevel Formulations

We will start this section by providing a general formulation for bilevel optimization. This is followed by various proposals that researchers have made for reducing a bilevel problem into a single-level problem. The two levels in a bilevel problem are also known as the leader’s (upper) and follower’s (lower) problems in the domain of game theory. In general, the variables, objectives and constraints are different for the two levels. The upper level variables are treated as parameters while optimizing the lower level problem. A general bilevel formulation has been provided below (for brevity, we ignore equality constraints):

Definition 1

For the upper-level objective function $F:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$ and lower-level objective function $f:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$, the bilevel optimization problem is given by

[TABLE]

where $G_{k}:X_{U}\times X_{L}\to\mathbb{R}$, $k=1,\dots,K$ denotes the upper level constraints, and $g_{j}:X_{U}\times X_{L}\to\mathbb{R}$ denotes the lower level constraints.

The bilevel problem in Definition 1 is ill-defined if there exists more than one lower level optimal solutions for some upper level variables. In such a case, the decision of the lower level remains unclear. However, in case where the lower level optimal solution is single-valued, the lone optimal solution is the rational choice. Optimizing bilevel problems from either optimistic or pessimistic position are two common approaches that researchers have utilized to handle the ambiguity arising from multiple lower level optimal solutions. In an optimistic position, it is assumed that the lower level chooses that optimal solution which is favorable at the upper level. In a pessimistic position, the upper level optimizes its problem according to the worst case scenario. In other words, the lower level may choose a solution from the optimal set that is least favorable at the upper level. In this paper, we assume an optimistic position while solving bilevel optimization problems.

In case when certain mathematically simplifying assumptions like continuities and convexities are satisfied, often the lower level optimization task in Definition 1 is replaced with its KKT conditions. However, the reduced formulation is not simple to handle, as it induces non-convexities and discreteness into the problem through the complementary slackness conditions. We do not utilize any properties of the KKT based reduction in this paper, rather we focus on two different formulations in the development of the evolutionary algorithm in this paper.

III-A Lower Level Reaction Set Mapping

The formulation provided in Definition 1 can also be stated as follows:

Definition 2

Let $\Psi:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{m}$ be the reaction set mapping,

[TABLE]

which represents the constraint defined by the lower-level optimization problem, i.e. $\Psi(x_{u})\subset X_{L}$ for every $x_{u}\in X_{U}$. Then the following gives an alternative formulation for the bilevel optimization problem:

[TABLE]

Using the above definition, a bilevel problem can be reduced to a single level constrained problem given that the $\Psi$-mapping can somehow be determined. Unfortunately this is rarely the case. Studies in the evolutionary computation literature that rely on iteratively approximation of this mapping to reduce the lower level optimization calls could be found in [44, 36, 37]. To illustrate the idea, let’s consider the Figure 2. To acquire sufficient data for constructing the $\Psi$-mapping approximation, a few lower level problems need to be optimized completely for their corresponding upper level decision vectors in the beginning. For instance, the lower level decisions for the upper level decisions $a,b,c,d,e$ and $f$ are determined by optimizing the lower level problem, which are then used to locally approximate the $\Psi$-mapping. This has been shown in Figure 2. Even though the actual $\Psi$-mapping is still unknown, the local approximation can then be substituted to identify the lower level optimal decision for every new upper level member to avoid the lower level optimization task. This procedure of approximating the mapping and utilizing it to predict the lower level optimum needs to be repeated iteratively until convergence to the bilevel optimum. The idea works well when the $\Psi$-mapping is single valued. If the lower level has multiple optimal solutions for some upper level members as shown in Figure 2, then identifying as well as approximating the mapping is not a straightforward task.

III-B Lower Level Optimal Value Function Mapping

Another formulation for the bilevel optimization problem in Definition 1 can be written using the optimal lower level value function: [56]:

Definition 3

Let $\varphi:X_{U}\to R$ be the lower level optimal value function mapping,

[TABLE]

which represents the optimal function value at the lower level for any given upper level decision vector. Using this lower level optimal value function, the bilevel optimization problem can be expressed as:

[TABLE]

As in the case of $\Psi$-mapping, if the $\varphi$-mapping can be somehow determined, a bilevel problem can be reduced to a single level problem as described in Definition 3. Along the process of an algorithm, the $\varphi$-mapping can be approximated and used to solve the reduced single level problem formulation in an iterative manner. Such an evolutionary algorithm has been recently discussed in [41]. The approximation of the optimal value function ($\varphi$) mapping is, in general, less complicated than the reaction set ($\Psi$) mapping, in the sense that, the $\varphi$-mapping is always scalar-valued regardless of the lower level variable dimension and whether or not there exist multiple lower level optimal solutions. However, the $\varphi$-mapping based reduction is not necessarily always better than the $\Psi$-mapping based reduction. Definition 3 requires the problem to be solved with respect to upper as well as lower level variables, while in Definition 2 the lower level variables are readily available from the $\Psi$-mapping. The $\Psi$-mapping based reduction also contains fewer constraints. Therefore, clearly there is a trade-off.

IV Evaluating the Performance of $\Psi$ and $\varphi$ mappings on test problems

In this section, we implement the $\Psi$ and $\varphi$ mappings separately in two different nested algorithms to evaluate the advantages and disadvantages of using the two mappings as a local search. For evaluating the two mappings, we choose a set of simple test problems that are provided in Tables X and XI. Firstly, we create a nested algorithm that utilizes an evolutionary approach for solving the upper level problem and sequential quadratic programming (SQP) for solving the lower level problem. Most of the lower level problems in the considered test cases being convex, explains the choice for using sequential quadratic programming (SQP) at the lower level. We enhance the nested approach by allowing it to approximate the $\Psi$ and $\varphi$ mappings and measure the performance gain provided by using each of the mappings separately. The implementation of the approaches has been outlined through the Figure 4. The flowchart without the overlapping box provides the steps involved in the nested approach. In case the idea involving $\Psi$ and $\varphi$ mappings has to be used then the local search (as mentioned in the overlapping box) is conducted every $k$ generations of the nested algorithm after the update step. A detailed description of the nested algorithm has been provided below.

Create a random population of size N comprising of upper level variables 2. 2.

Solve the lower level optimization problem using SQP for each upper level variable. 3. 3.

Evaluate the fitness of each population member using upper level function and constraints (refer to Section VI-B) 4. 4.

Choose $2\mu$ population members using tournament selection and apply genetic operators (refer to Section VI-C) to produce $\lambda$ offspring. 5. 5.

Solve the lower level optimization problem using SQP for each offspring. 6. 6.

Evaluate the fitness of each offspring using using upper level function and constraints 7. 7.

Form a pool consisting of $r+\lambda$ members, where $r$ members are chosen randomly from the population and $\lambda$ members are the offspring. Use the best $r$ members from this pool to replace the chosen $r$ members from the population. 8. 8.

Perform a termination check (refer to Section VI-D) and proceed to Step 5 if termination check is false, otherwise stop.

The parameters used in the implementation of the above procedure are $N=50,\mu=2,\lambda=3$ and $r=2$. The lower level SQP terminates when the improvement in the lower level function value is less than 1e-6.

IV-A Approximating the $\Psi$-mapping

Let $\mathcal{H}$ be the hypothesis space. The hypothesis space consists of all functions that can be used to generate a mapping between the upper level decision vectors and optimal lower level decision vectors. Given a sample $\mathcal{S}$, consisting of upper level points and corresponding optimal lower level points, we would like to identify a model $\hat{\Psi}\in\mathcal{H}$ that minimizes the empirical error on the sample, i.e.

[TABLE]

where $L:X_{L}\times X_{L}\to\mathbb{R}$ denotes the prediction error, $x_{u}^{(i)}$ is any given upper level vector and $\bar{x}_{l}^{(i)}$ is its corresponding optimal solution. The prediction error may be calculated as follows:

[TABLE]

We have restricted the hypothesis space $\mathcal{H}$ to consist of second-order polynomials which reduces the error minimization problem to an ordinary quadratic regression problem. The sample $\mathcal{S}$ can be created from the population members or an archive. It should be noted that this can approximate only single-valued mapping and will fail if the mapping becomes set-valued.

IV-B Approximating the $\varphi$-mapping

Once again, let $\mathcal{H}$ be the hypothesis space of functions and $\mathcal{S}$ be a sample of upper and corresponding lower level points, our aim is to identify a model $\hat{\varphi}\in\mathcal{H}$ that minimizes the empirical error on the sample, i.e.

[TABLE]

where $L:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ denotes the prediction error, $x_{u}^{(i)}$ is any given upper level vector and $\bar{f}^{(i)}$ is its corresponding optimal function value. The prediction error can once again be computed as follows:

[TABLE]

We have once again restricted the hypothesis space $\mathcal{H}$ to consist of second-order polynomials. Since the $\varphi$-mapping is always single valued, approximating it will not involve similar issues as for the $\Psi$-mapping.

V Comparison Results for $\Psi$ vs $\varphi$ Approximations

For comparing $\Psi$-approximation approach against $\varphi$-approximation approach, we use a set of 8 test problems selected from the literature given in Tables X and XI. Table I compares the median function evaluations at both level for three algorithms $\Psi$-approximation, $\varphi$-approximation and nested algorithm. The results have been produced from 31 runs of the algorithm and further details about the runs can be found in Figures 6 and 6. Both $\Psi$-approximation and $\varphi$-approximation perform equally well and outperform the nested approach in this study. The differences in the performance of $\Psi$-approximation and $\varphi$-approximation can be attributed to differences in the quality of approximation produced during the intermediate steps of the algorithm. In Table II, we provide a comparison of the meta-modeling results with other evolutionary approaches [50, 51] to provide an idea about the extent of savings that can be produced using meta-modeling techniques. The advantage is quite clear as the savings are better by multiple order of magnitudes on the set of test problems considered in this study.

It should be noted that the $\Psi$ approximation idea would fail if the $\Psi$-mapping in bilevel optimization is set valued. Next, we test this hypothesis, by modifying the 8 test problems such that each test problem necessarily has a set-valued $\Psi$-mapping. To achieve this, we add two additional lower level variables ($y_{p}$ and $y_{q}$) in each test problem. Both the upper and lower level functions are modified as shown below:

[TABLE]

The modification makes the lower level problem have infinitely many optimal solutions (for all $y_{p}=y_{q}$) for any given upper level vector. Out of the many optimal solutions the upper level prefers the solution where $y_{p}=y_{q}=0$. After this simple modification, we once again solve the test problems using $\varphi$-approximation and $\Psi$-approximation approaches. As shown in Tables III and IV, the $\varphi$-approximation algorithm still works but $\Psi$-approximation algorithm completely fails. Function evaluations for $\varphi$-approximation algorithm increases slightly than before because of additional variables in the problem.

Therefore, the $\varphi$-approximation idea clearly has an advantage over the $\Psi$-approximation idea. Moreover, $\varphi$-mapping is always a scalar valued mapping as compared to $\Psi$-mapping which is usually vector valued and can also be set-valued. However, there is a trade-off. The reduced single level problem formed using $\Psi$-mapping may usually be a little easier to handle as compared to the single level problem formed using $\varphi$-mapping. The reason being that in case of $\Psi$-mapping the lower level variables are readily available, and the reduced problem does not involve lower level constraints. For $\varphi$-mapping, the reduced problem has to be solved both with respect to upper and lower level variables, and the formulation involves both upper and lower level constraints. Given the pros and cons of using the two mappings, next, we would like to develop an evolutionary algorithm that is capable of utilizing the better of the two mappings while solving a bilevel optimization problem.

VI Bilevel Evolutionary Algorithm based on $\Psi$ and $\varphi$-mapping Approximations

In this section, we provide the bilevel evolutionary algorithm that approximates the $\Psi$ as well as the $\varphi$ mapping during the intermediate steps of the algorithm. From the previous experiments and the properties of the two mappings we infer that there can be situations when the approximation of the $\Psi$-mapping may fail, while when $\Psi$-mapping can be approximated it offers the advantage of completely ignoring the lower level functions and constraints. Acknowledging this fact, we utilize both the approximations in our algorithm. The algorithm adaptively decides to use one of the mappings based on the quality of fit obtained when approximating the two mappings. Local quadratic approximations are created for the two mappings from a sample of points in the vicinity of the point around which we want to create an approximation. Introducing local approximation, is expected to improve the quality of approximations significantly. Given a sample dataset, the steps for creating an approximation are the same as discussed in Sections IV-A and IV-B. The algorithm also maintains an archive so as to maintain a large dataset for creating and validating the approximations. Deviating from the nested algorithm, we employ the approximated $\Psi$ and $\varphi$ mappings to avoid frequent lower level optimization calls. An earlier version of the algorithm [44, 36, 37] that relied on $\Psi$-mapping approximation alone was referred as Bilevel Evolutionary Algorithm based on Quadratic Approximations (BLEAQ). We keep the same terminology and refer to the newer version of the algorithm as BLEAQ-II. The pseudocode for the algorithm has been provided in Table V.

VI-A Initialization

The initialization in the algorithm is done by creating random upper level members $x_{u}^{(1)},\dots,x_{u}^{(N)}$, and then solving the lower level optimization problem for each member to get optimal $x_{l}^{(1)},\dots,x_{l}^{(N)}$. There can be situations, where finding random feasible upper and lower level pair that satisfy both lower and upper level constraints in the problem can be difficult. In such situations, one can solve the following problem to create $(x_{u}^{(i)},x_{l}^{(i)})$ pairs that satisfies all the constraints to begin with.

[TABLE]

The above problem can be solved using any standard procedure like a greedy GA or SQP with a random starting point to arrive at a feasible solution. As soon as a feasible member is found, the method stops. Solving the above method repeatedly with random population (in case of GA) or a random starting point (in case of SQP) can provide the starting population of upper level members $x_{u}^{(1)},\dots,x_{u}^{(N)}$ for the BLEAQ-II algorithm. For this given set of upper level members, we know that at least one feasible lower level member exists and we still need to solve the lower level problem to find the optimal lower level solutions $x_{u}^{(1)},\dots,x_{u}^{(N)}$. 111In case of upper level constraints containing both upper and lower level variables, one can find it difficult to arrive at a $(x_{u}^{(i)},x_{l}^{(i)})$ pair that is feasible with respect to all the constraints and the lower level vector is optimal for the given upper level vector. Many formulations of bilevel optimization, therefore, do not consider lower level variables in upper level constraints.

VI-B Constraint handling and Fitness Assignment

The proposed approach always assigns higher fitness to a feasible member over a non-feasible member. For two given members, $(x_{u}^{(i)},x_{l}^{(i)})$ and $(x_{u}^{(j)},x_{l}^{(j)})$, if both members are feasible with respect to constraints then it looks at the function value. If both members are infeasible, then it looks at the overall constraint violation. This fitness assignment scheme is similar to the one proposed in [16]. At the lower level the idea can be implemented directly using lower level constraints and lower level function value. At the upper level, for any given upper and lower level pair, we only consider upper level function and constraints, without considering if the corresponding lower level vector is optimal. The information about a lower level vector corresponding to an upper level vector being optimal is stored using tagging (0 or 1).

VI-C Genetic operations

Offspring are produced in the BLEAQ-II approach using standard crossover and mutation operators. Genetic operations at the upper level involve only upper level variables, and the operations at the lower level involve only lower level variables. We utilize parent centric crossover (PCX) and polynomial mutation for generating the offspring. The crossover operator used in the algorithm is similar to the parent-centric crossover (PCX) operator proposed in [43]. The operator uses three parents and produces offspring around the index parent as described below.

[TABLE]

where,

• •

$z^{(p)}$ is the index parent (the best parent among three parents)

• •

$d=z^{(p)}-g$, where $g$ is the mean of $\mu$ parents

• •

$p^{(1)}$ and $p^{(2)}$ are the other two parents

• •

$\omega_{\xi}=0.1$ and $\omega_{\eta}=0.1$ are the two parameters.

VI-D Termination Criteria

A variance based termination criterion has been used at both levels; some other termination criterion like termination based on no improvement may also be used. Variance based termination allows the algorithm to terminate automatically when the variance of the population becomes small. At the upper level, the variance of the population at any generation, $T$, is computed as follows:

[TABLE]

When the value of $\alpha_{u}^{T}$ at any generation $T$ becomes less than $\alpha_{u}^{stop}$ then the algorithm terminates. In the above equation, $n$ is the number of upper level variables, $\sigma^{2}(x_{i})|_{T}$ is the variance across dimension $i$ at generation $T$ and $\sigma^{2}(x_{i})|_{0}$ is the variance across dimension $i$ in the initial population. A similar termination scheme can be used when the evolutionary algorithm is executed at the lower level.

VI-E Lower Level Optimization

At the lower level, we utilize SQP if the problem is convex, otherwise we use the lower level evolutionary algorithm described in Table VI that uses similar genetic operations as used at the upper level.

VI-F Offspring Update

For an offspring $x^{(j)}=(x^{(j)}_{u},x^{(j)}_{l})$, the lower level vector $x^{(j)}_{l}$ is updated either using $\Psi$-approximation or $\varphi$-approximation. An update using $\Psi$-approximation is straightforward. However, if an update has to be done using $\varphi$-approximation, it requires to solve the following auxiliary optimization problem. In the auxiliary problem $x_{u}$ is fixed as $x^{(j)}_{u}$ and the problem is solved only with respect to $x_{l}$. The optimal $x_{l}$ replaces the lower level vector $x^{(j)}_{l}$ of the offspring.

[TABLE]

In the above formulation we use hat for all the functions and constraints as we solve the auxiliary problem on approximated functions and constraints. We use linear approximations for all the constraints, while quadratic approximation is used for the other functions. The auxiliary problem may have to be solved frequently if the lower level problem contains multiple optimal solutions. Solving the auxiliary problem with approximated functions helps in saving actual function evaluations. Note that in the ideal case the auxiliary problem will lead to an optimistic lower level solution corresponding to the fixed $x^{(j)}_{u}$.

VI-G Local Search

The algorithm utilizes local search after every $k$ generations of the algorithm. The local search is performed by meta-modeling the upper and lower level functions and constraints along with the $\Psi$ and the $\varphi$-mappings in the vicinity of the best member in the population. Once the $\Psi$ and the $\varphi$-mappings are available, the quality of the two mappings are assessed by the mean square error of the approximations (i.e. $e_{mse}^{\Psi}$ and $e_{mse}^{\varphi}$). The better mapping and the corresponding single level reduction (described in Section III-A and III-B) with approximated functions is solved using SQP to arrive at $x_{u}^{(LS)}$. A lower level optimization corresponding to $x_{u}^{(LS)}$ is solved and if the member is found to be better than $x_{best}^{(j)}$ then $x_{best}^{(j)}$ is updated. In case the member is not better than the best member found so far, then the next local search is performed using the exact upper/lower level objective functions and constraints.

VI-H Parameters and Platform

The algorithm has been implemented in MATLAB. At the upper and lower level, the parameters used in the algorithm are:

$\mu=3$ 2. 2.

$\lambda=2$ 3. 3.

$r=2$ 4. 4.

Probability of crossover = 0.9 5. 5.

Probability of mutation = 0.1 6. 6.

N = 50 (Population size at upper level) 7. 7.

n = 50 (Population size at lower level)

VII Results

In this study, we consider three algorithms, the nested approach described in Figure 4, BLEAQ [44, 36, 37], and our proposed BLEAQ-II. To assess the performance of each algorithm, 31 runs have been performed for each test instance. During every simulation process, the algorithm is terminated if the function value accuracy reaches the objective function accuracy of $10^{-2}$ at both levels. For each run, the upper and lower level function evaluations required until termination is recorded separately.

VII-A Standard test problems

We first present the empirical results on 8 standard test problems selected from the literature (referred to as TP1-TP8). The description for these test problems has been provided in the Appendix A. Table VII contains the median upper level (UL) function evaluations, lower level (LL) function evaluations and BLEAQ-II’s overall function evaluation savings as compared to other approaches from 31 runs of the algorithms. The overall function evaluations for any algorithm is simply the sum of upper and lower level function evaluations. For instance, for the median run with TP1, BLEAQ-II requires $63\%$ less overall function evaluations as compared to BLEAQ, and $98\%$ less overall function evaluations as compared to the nested approach.

All these test problems are bilevel problems with small number of variables, and all the three algorithms were able to solve the 8 test instances successfully. A significant computational saving can be observed for both BLEAQ-II and BLEAQ, as compared to the nested approach as shown in the Savings column of Table VII. The performance gain going from BLEAQ to BLEAQ-II is quite significant for these simple test problems even though none of them lead to multiple lower level optimal solutions. Detailed comparison between BLEAQ and BLEAQ-II in terms of upper and lower level function evaluations is provided through Figures 8 and 8.

VII-B Scalable test problems

Next, we compare the results for the three algorithms on the scalable SMD test suite which contains 12 test problems in the original paper [38]. We extend this test suite in this paper to a set of 14 test problems by adding two additional scalable test problems. The description for the additional SMD test problems can be found in Appendix B. First we analyze the performance of the algorithms on a smaller version of the test problems which consists of 5 variables, and then we provide the comparison results on 10-variable instances of the SMD test problems. For the 5 variable version of the SMD test problems, we used the settings as $p=1$, $q=2$ and $r=1$ for all SMD problems except SMD6 and SMD14. For the 5 variable version of SMD6 and SMD14, we used $p=1$, $q=0$, $r=1$ and $s=2$. For the 10 variable version of the SMD test problems, we used the settings as $p=3$, $q=3$ and $r=2$ for all SMD problems except SMD6 and SMD14. For the 10 variable version of SMD6 and SMD14, we used $p=3$, $q=1$, $r=2$ and $s=2$.

Table VIII provides the median function evaluations and overall savings for the three algorithms on the set of 14 SMD problems. These test problems contain 2 variables at the upper level and 3 variables at the lower level and offer a variety of tunable complexities to the algorithms. For instances, the test set contains problems which are multimodal at the upper and the lower levels, contain multiple optimal solutions at the lower level, contain constraints at the upper and/or lower levels etc. It can be found that BLEAQ-II is able to solve the entire set of 14 SMD test problems, while BLEAQ fails on 2 test problems. The overall savings with BLEAQ-II is higher as compared to BLEAQ for all the test problems. The test problems that contain multiple lower level solutions include SMD6 and SMD14, for which BLEAQ is unable to handle the problem. Further details about the required overall function evaluations from 31 runs can be found in Figures 9.

Results for the high dimensional SMD test problems have been provided in Table IX. BLEAQ-II leads to much higher savings as compared to BLEAQ, and with higher dimensions BLEAQ is found to once again fail on SMD6 and also on SMD7 and SMD8. Both methods outperform the nested method on most of the test problems. We do not provide results for SMD9 to SMD14 as none of the algorithms were able to handle these problems. It is noteworthy that SMD9 to SMD14 offer difficulties like multi-modalities and highly constrained regions, which none of the algorithms were able to handle with the parameter setting used in this paper. Details for the 31 runs on each of these test problems can be found in Figure 10.

Through Figures 12 and 12, we provide the quality of prediction of the lower level optimal solution made by the $\Psi$-mapping and $\varphi$-mapping approach over the course of the algorithm. It is interesting to note that the quality of $\varphi$-approximation is better in the case of SMD1 test problem in Figure 12, therefore, the prediction decisions are mostly made using the $\varphi$-approximation approach. However, for SMD13 in Figure 12, which involves a difficult $\varphi$-mapping, the prediction decisions are made using the $\Psi$-approximation approach. Both these mappings are found to be improving with an increase in generations of the algorithm. The two figures show the adaptive nature of the BLEAQ-II algorithm in choosing the right approximation strategy based on the difficulties involved in a bilevel optimization problem.

VIII Conclusions

In this paper, we have presented a computationally efficient evolutionary algorithm for solving bilevel optimization problems. The algorithm is based on iterative approximations of two important theoretically motivated mappings; namely, the lower level rational reaction mapping and the lower level optimal value function mapping. The paper discusses about the pros and cons of utilizing these mappings in an evolutionary bilevel optimization algorithm by embedding them in a nested approach. Thereafter, an algorithm is developed that adaptively decides to use one of the mappings during the execution based on the characteristics of the bilevel optimization problem being solved. The proposed algorithm has been tested on a wide variety of bilevel test problems and it has been able to perform significantly better than other approaches in terms of computational requirements.

Appendix A Standard Test Problems

In this section, we provide some of the standard bilevel test problems chosen from the literature. Most of these test problems are small with only small number of variables at both levels.

Appendix B Additional SMD Test Problems

SMD test problems [38] are a set of 12 scalable test problems that offer a variety of controllable difficulties to an algorithm. We add two more test problems to the previous test-suite in this paper. Both these problems contain a difficult $\varphi$-mapping, among other difficulties. The upper and lower level functions follow the following structure to induce difficulties due to convergence, interaction, and function dependence between the two levels. The vectors $x_{u}$ and $x_{l}$ are further divided into two sub-vectors. The $\varphi$-mapping is defined by the function $f_{1}$.

[TABLE]