An interior-point trust-region algorithm to solve a nonlinear bilevel programming problem

1.   Introduction

Bilevel programming problem has increasingly been addressed in the literature, both from the theoretical and computational points of view ^[14]. This model has been widely applied to decentralized planning problems involving a decision progress with a hierarchical structure. It is characterized by the existence of two optimization problems in which the constraint region of the first-level problem is implicitly determined by another optimization problem. The NBLP problem is hard to solve. In fact, the problem has been proved to be NP-hard ^[8]. However, the NBLP problem is used so extensively in resource allocation, finance budget, price control, transaction network etc. ^{[1,7,28,29,39]} that many researches have been devoted to this field, which leads to a rapid development in theories and algorithms. For the detailed expositions, the reader may consult ^[21,33].

In this paper we will consider the following NBLP problem

where $t\in \Re^{n_1}$ and $y\in \Re^{n_2}$. The functions $f_u: \Re^{n_1+n_2} \rightarrow \Re$, $f_l: \Re^{n_1+n_2} \rightarrow \Re$, $g_u: \Re^{n_1+n_2} \rightarrow \Re^{m_1}$, and $g_l: \Re^{n_1+n_2} \rightarrow \Re^{m_2}$ are assumed to be at least twice continuously differentiable function.

There are several approaches have proposed to solve problem 1.1, see ^{[2,3,25,35,40]}. One of these approaches and used in this paper, is converted the original two level problems to a single level one by replacing the lower level optimization problem with its Karush-Kuhn-Tucker (KKT) conditions, see ^[24,41]. By KKT optimality conditions for the lower-level problem, then we can reduce the NBLP problem 1.1 to one-level programming problem. This problem is non-convex and non-differentiable, moreover the regularity assumptions which are needed to successfully handle smooth optimization problems are never satisfied and it is not good to use our approach. So, we add slack variables for inequalities constraints in problem 1.1.

By adding slack variables $s_u\in \Re^{m_1}$ and $s_l\in \Re^{m_2}$ to the upper inequality constraint $g_u(t, y)$ and the lower inequality constraint $g_l(t, y)$ respectively, then NBLP problem 1.1 can be written as follows

The above NBLP problem can be simplified as follows

where $\tilde{g}_u(t, y, s_u) = g_u(t, y)+s_u\in \Re^{m_1}$, $\tilde{g}_l(t, y, s_l) = g_l(t, y)+s_l\in \Re^{m_2}$, and $s = (s_u, s_l)^T\in \Re^{m_1+m_2}$.

Applying KKT conditions only on the lower-level problem without the constraint $s\geq 0$, then we can reduce the NBLP problem 1.2 to the following smooth SONP problem:

where $\mu_l \in \Re^{m_2}$ is a Lagrange multiplier vector associated with equality constraint $\tilde{g}_l(t, y, s_l)$, see ^[5].

Using problem 1.3, to overcome the difficulty that problem 1.1 does not satisfy any regularity assumptions, which are needed for successfully handling smooth optimization problems, and pave the way for using the proposed approach to solve problem 1.1. To simply our discussion, we introduce the following notations. $x = (t, y, s)^T\in \Re^n$, $n = {n_1+n_2+m_1+m_2}$, $h(x)\in\Re^m$ represents the vector of equality constraints such that $m = m_1+m_2+n_2$. Then problem 1.3 can be written as follows

where $v\in\{\Re\bigcup\{-\infty\}\}^{n}$, $w\in\{\Re\bigcup\{+\infty\}\}^{n}$, and $v < w$.

Various approaches have been proposed to solve the SONP problem 1.4, see ^{[5,9,10,11,15,16,17,18,19]}. In this paper, we use Newton's interior point method with Das scaling matrix ^[12] to solve problem 1.4. Newton's method converges quadratically to a stationary point under reasonable assumptions if the starting point sufficiently closed to the stationary point. It may not converge if the starting point is far away from the stationary point. To guarantee convergence from any starting point, a trust-region strategy is used. The trust-region strategy can induce strongly global convergence, which is very important method for solving SONP problem and is more robust when it deals with rounding errors. It does not require the objective function of the model be convex or the Hessian of the objective function must be positive definite. Also, some criteria are used to test the trial step is acceptable or no. If it is not acceptable, then the subproblem must be resolved with a reduced the trust-region radius. For the detailed expositions, the reader may consult ^{[4,17,20,21,22,23,30,32,36,42,43,45,46]}.

A reduced hessian technique is used in this paper to overcome some difficulties in trust-region subproblem. This technique was suggested by ^[6,37] and used by ^[19,20].

In this paper, we use the symbol, $f_{u_k} \equiv f_u(x_k)$, $h_k \equiv h(x_k)$, $P_k \equiv P(x_k)$ $\ell_k \equiv \ell(x_k, \lambda_k)$, $\nabla_x\ell_k \equiv \nabla_x\ell(x_k, \lambda_k)$, and so on to denote the function value at a particular point. Finally, all norms are $l_2$-norms.

The rest of the paper is organized as follows. In section 2, we introduce detailed description for the proposed method to solve problem 1.4. Section 3 is devoted to analysis of the global convergence of the proposed algorithm. Section 4 contains implementation of the proposed algorithm and the results of test problems. Section 5 contains concluding remarks.

2.   An interior-point method with trust-region algorithm

In this section, firstly, we will consider the detailed description for the Newton's interior-point method with Das scaling matrix to solve SONP problem 1.4. Secondly, to guarantee convergence from any starting point, we will introduce the detailed description for trust-region strategy. Finally, we clarify main steps for general algorithm to solve NBLP 1.1.

2.1.   Newton's method with scaling matrix

Motivated by the impressive computational performance of Newton's interior-point method for solving SONP problem 1.4, let

be a Lagrangian function associated with problem 1.4 without the constraints $v \leq x \leq w$, and let

be a Lagrangian function associated with problem 1.4 with the constraints $v \leq x \leq w$. The vectors $\lambda\in \Re^m$, $\mu^v\in \Re^n$, and $\mu^w\in \Re^n$ represent Lagrange multiplier vectors associated with the constraints $h(x) = 0$, $0\leq(x-v)$, and $0\leq(w- x)$ respectively. Let $\hat{\boldsymbol{G}} = \{x: v \leq x \leq w \}$ and $int(\hat{\boldsymbol{G}}) = \{x : v < x < w \}$.

The first-order necessary conditions for the point $x_*$ to be a local minimizer of problem 1.4 are the existence of multipliers $\lambda_* \in \Re^m$, $\mu^v_* \in \Re_{+}^n$, and $\mu^w_* \in \Re_{+}^n$, such that $(x_*, \lambda_*, \mu^v_*, \mu^w_*)$ satisfies

and for all $i$ corresponding to $x^{(i)}$ with finite bound, we have

where $\nabla_x \ell(x_*, \lambda_*) = \nabla f_u(x_*)+\nabla h(x_*) \lambda_*$.

The proposed algorithm here, like its predecessors in ^[12,18,19], starts at a point strictly feasible with respect to the bounds on the variables and produces iterates that are strictly feasible with respect to the bounds (i.e. 'in the interior'). Define a diagonal scaling matrix $P(x) = diag(p(x)$ whose diagonal elements $p(x)$ are given by

Using the matrix P(x), then $(x_*, \lambda_*, \mu^v_*, \mu^w_*)$ satisfy the KKT conditions [2.3-2.7] if and only if

For more details about the proof, see ^[12].

Applying Newton's method on the nonlinear system [2.9-2.10], then we have

where $\theta(x)$ is a vector whose components are given by

For more details see ^[18].

In our method, the matrix $P(x)$ must be nonsingular, so we restrict the point $x\in int (\hat{\boldsymbol{G}})$. Multiplying both sides of equation 2.11 by $P^{-1}(x)$, then we have

Substituting $\Delta x = P(x)\; d$ in the above two system, then we have

where $H(x, \lambda) = \nabla^2_x \ell(x, \lambda)$ represents the Hessian of the Lagrange function 2.1 or an approximation to it. It is easy to see that the step generated by the above system coincides with the solution of the following quadratic programming subproblem

where $B = P(x) H(x, \lambda) P(x) + diag(\nabla_x \ell(x, \lambda)) diag(\theta(x))$. This means that, the point $(x_*, \lambda_*)$ that satisfies the KKT conditions for subproblem 2.16 will satisfy the KKT conditions for problem 1.4.

Although Newton's method converges quadratically to a stationary point under reasonable assumptions, it may not converge to a stationary point if the starting point is far away from the solution. To overcome this disadvantage and to guarantee convergence from any starting point, we use the trust-region technique.

2.2.   Trust-region technique

Trust-region methods can induce strongly global convergence, which are very important methods for solving a smooth nonlinear programming problem and are more robust when they deal with rounding errors. It does not require the objective function of the model be convex. Also, it does not demand the Hessian of the objective function must positive definite.

The trust-region subproblem associated with problem 2.16 is

where $\delta_k > 0$ is the radius of the trust-region.

Subproblem 2.17 may be infeasible, because there may be no intersecting points between the constraint $\|d\|\leq \delta_k$ and $h_k+(P_k \nabla h_k)^T d = 0$ constraints. Even if they intersect, there is no warranty that this will continue true if $\delta_k$ is decreased. For more details see ^[13]. To overcome this difficulty, we use a reduced hessian technique. This technique was suggested by ^[6,37] and used by ^[19,20]. In this technique, the trial step $d$ is decomposed into two orthogonal components: the normal component $d^n$ to improve feasibility and the tangential component $d_k^t$ to improve optimality. Each of components is computed by solving unconstrained trust-region subproblem.

● How to estimate the normal component $d_k^n$

The normal component $d_k^n$ is computed by solving the following trust-region subproblem

for some $0 < \zeta < 1$. To solve the subproblem 2.18, we use a conjugate gradient method which is introduced by ^[38] and used by ^[21], see algorithm 2.1 in ^[21]. It is very cheap if the problem is large-scale and the Hessian is indefinite. By using the conjugate gradient method, the normal predicted decrease obtained by $d_k^n$ is greater than or equal to a fraction of the normal predicted decrease obtained by the Cauchy step $d_k^{ncp}$. This means that

such that $d_k^{ncp}$ is defined as follows

where the parameter $\varphi_k^{ncp}$ is given by

Once $d_k^n$ is obtained, we will evaluate $d_k^t = Z_k \bar{d}_k^t$ where $Z_k$ is the matrix whose columns form a basis for the null space of $(P_k \nabla h_k)^T$.

● How to estimate the tangential component $d_k^t$

To obtain the tangential component $d_k^t$, we use the conjugate gradient method ^[21] to solve the following trust-region subproblem

where $\nabla q_k(P_k d_k^n) = P_k \nabla_x \ell_k+ B_k d_k^n$ and $\Delta_k = \sqrt{\delta_k^2 - \| d_k^n \|^2}$.

By using the conjugate gradient method, the tangential predicted decrease which is obtained by tangential step $\bar{d}_k^t$ is greater than or equal to a fraction of the tangential predicted decrease obtained by a tangential Cauchy step $\bar{d}_k^{tcp}$. This means that

for some $0 < \vartheta_2 \leq 1$ and $\bar{d}_k^{tcp}$ is defined as follows

where the parameter $\varphi_k^{tcp}$ is given by

where $\bar{B}_k = Z_k^T B_k Z_k$.

● How to estimate a parameter $\gamma_k$

Once obtaining $d_k^t$, we set $d_k = d_k^n+d_k^t$ and $x_{k+1} = x_k+P_k d_k$. To ensure $x_{k+1}\in int(\hat{\boldsymbol{G}})$, we need to evaluate the parameter $\gamma_k$. To do this, evaluate

and

Compute

Once the trial step $\gamma_k P_k d_k$ is evaluated, it needs to be tested to decide whether it will be accepted or not. To do this, we need to a merit function which is ties the objective function and the constraints in such a way that progress in the merit function means progress in solving problem. In our method, we use the following merit function which is introduced by ^[26] and known as an augmented Lagrange function

where $\ell(x, \lambda)$ is defined in 2.1 and $\rho > 0$ represents the penalty parameter.

● How to estimate $\lambda_{k+1}$

The Lagrange multiplier vector $\lambda_{k+1}$ will be estimated as follows

To test whether the point $(x_{k+1}, \lambda_{k+1})$, will be accepted in the next iterate or no we need to define the following actual reduction $Ared_k$ and the predicted reduction $Pred_k$.

The actual reduction $Ared_k$ in the merit function 2.27 in moving from $(x_k, \lambda_k)$ to $(x_k+\gamma_k P_k d_k, \lambda_{k+1})$ is defined as follows

Also we can write the actual reduction $Ared_k$ as follows,

where $\Delta \lambda_k = (\lambda_{k+1} - \lambda_k)$.

The predicted reduction $Pred_k$ is defined as follows

Since $q_k(\gamma_k P_k d_k) = \ell_k + (P_k \nabla_x \ell_k)^T \gamma_k d_k + \frac{1}{2}\gamma_k^2 d_k^T B_k d_k$, then $Pred_k$ can be written as follows,

● How to update the penalty parameter $\rho_k$

To ensure $Pred_k$ is strictly positive, we use the following scheme to update the positive penalty parameter $\rho_k$

Algorithm 2.1.: (Updating the penalty parameter $\rho_k$)

Set $\rho_{k+1} = \rho_k$.

If

then set

End if.

● How to test the step $\gamma_k P_k d_k$ and update $\delta_k$

To decide the trial step $\gamma_k P_k d_k$ will be accepted in the next iteration or no, we use the following algorithm.

Algorithm 2.2.: (Testing the step $\gamma_k P_k d_k$ and updating $\delta_k$)

Step 0. Choose $0 < \tau_1 < \tau_2 < 1$, $0 < \beta_1 < 1 < \beta_2$, and $\delta_{min}\leq \delta_0 \leq \delta_{max}$.

Step 1. While $\frac{Ared_k}{Pred_k} < \tau_1$ or $Pred_k \leq 0$.

Do not accept the step and set $\delta_k = \beta_1\|d_k\|$.

Compute a new trial step.

End while.

Step 2. If $\tau_1\leq\frac{Ared_k}{Pred_k} < \tau_2$.

Accept the step: $x_{k+1} = x_k + \gamma_k P_k d_k$.

Set $\delta_{k+1} = \max(\delta_k, \delta_{min})$.

End if.

Step 3. If $\frac{Ared_k}{Pred_k}\geq \tau_2$.

Accept the step: $x_{k+1} = x_k + \gamma_k P_k d_k$.

Set $\delta_{k+1} = \min\{\delta_{max}, \max\{\delta_{min}, \beta_2\delta_k\}\}$.

End if.

Finally, the algorithm is stopped when either $\| Z_k^T P_k \nabla_x \ell_k \| + \| h_k\| \leq \varepsilon_1$, for some $\varepsilon_1 > 0$ or $\|d_k\|\leq \varepsilon_2$ for some $\varepsilon_2 > 0$.

Main steps of the trust-region algorithm for solving subproblem 2.17 are summarized in the following algorithm.

Algorithm 2.3. (Trust-region algorithm)

Step 0. Starting with $x_0\in int (\hat{\boldsymbol{G}})$. Evaluate $\lambda_0$, $P_0$, and $\beta_0$. Set $\rho_0 = 1$ and $c_0 = 0.1$.

Choose $\varepsilon_1 > 0$, $\varepsilon_2 > 0$, and set $k = 0$.

Step 1. If $\| Z_k^T P_k \nabla_x \ell_k \| + \| h_k\| \leq \varepsilon_1$, then stop.

Step 2. (To compute $d_k$)

a) Compute $d_k^n$ by solving trust-region subproblem 2.18.

b) Compute $\bar{d}_k^t$ by solving trust-region subproblem 2.22.

c) Set $d_k = d_k^n + Z_k \bar{d}_k^t$.

Step 3. If $\| d_k\| \leq \varepsilon_2$, then stop.

Step 4. Compute $\gamma_k$ using equation 2.26.

Step 5. Update $\lambda_{k+1}$ using subproblem 2.28.

Step 6. Update the penalty parameter using scheme 2.1.

Step 7. Test the step $\gamma_k P_k d_k$ and update $\delta_k$ by using algorithm 2.2.

Step 8. Compute $P_{k+1}$ and $\alpha_{k+1}$ using definitions 2.8 and 2.13 respectively.

Step 9. Set $k = k+1$ and go to Step 1.

Main steps for solving NBLP problem1.1 are summarized in the following algorithm.

Algorithm 2.4. (Interior-point trust-region (IPTR) algorithm)

Step 1. Adding slack variables to inequalities in NBLP problem1.1 and convert it to problem 1.2.

Step 2. By KKT optimality conditions for the lower-level problem, NBLP problem 1.2 is equivalent to the one level problem 1.3 which can be written in the form 1.4.

Step 3. Using Newton's method and Das strategy to transform problem 1.4 to subproblem 2.16.

Step 4. Using trust-region algorithm 2.3 to solve subproblem 2.16.

The following section is devoted to global convergence analysis for IPTR algorithm 2.4.

3.   Global convergence theory

We state the general assumption under which the global convergence theory for IPTR algorithm 2.4 is proved.

3.1.   A general assumptions

Let $\Omega$ be a convex subset of $\Re^n$ that contains all points $x_k \in int (\hat{\boldsymbol{G}})$ and $(x_k+\gamma_k P_k d_k) \in int (\hat{\boldsymbol{G}})$. On the set $\Omega$ we state the following general assumptions under which the global convergence theory of IPTR algorithm is proved

[$GS_1$.] The functions $f_u(x)$, $h(x)\in \mathcal{C}^2$ for all $x \in \Omega$.

[$GS_2$.] The matrix $P_k \nabla h_k$ has full column rank.

[$GS_3$.] All of $f_u(x)$, $\nabla f_u(x)$, $\nabla^2 f_u(x)$, $h(x)$, $\nabla h(x)$, $\nabla^2 h_i(x)$ for $i = 1, ..., m$ and $(P_k\nabla H_k)((P_k\nabla h_k)^T(P_k\nabla h_k))^{-1}$ are uniformly bounded in $\Omega$.

[$GS_4$.] The sequence of Lagrange multiplier vectors $\{\lambda_k\}$ is bounded.

[$GS_5$.] The sequence of approximate Hessian matrices $\{H_k\}$ is bounded.

An immediate consequence of the above general assumptions is that the existence of positive constant $b_1$, such that

3.2.   Technical lemmas

In this section, we introduce some important results which are needed in the subsequent proof.

The following lemma shows how accurate the definition of $Pred_k$ is as an approximation to $Ared_k$.

Lemma 3.1. Under assumptions $GS_1$-$GS_5$, there exists a positive constant $K_1$, such that

Proof. From Equations 2.29, 2.30, and using the inequality of Cauchy-Schwarz, we have

for some $\xi_1$ and $\xi_2 \in (0, 1)$. Using the general assumptions $GS_1-GS_5$ and $0 < \gamma_k \leq 1$, we have

where $\kappa_1$, $\kappa_2$, and $\kappa_3$ are positive constants. Since $\rho_k \geq 0$, $\| d_k\|\leq \delta_{max}$, and $\| h_k \|$ is uniformly bounded, then inequality 3.2 hold.

The following lemma obviously that the normal predicted reduction at any iteration k, is at least equal to the decrease in the 2-norm of the linearized constrained by the Cauchy step

Lemma 3.2. Under assumptions $GS_1$-$GS_5$, there exists a constant $K_2 > 0$, such that

Proof. Since $d_k^n$ is obtained by approximating the solution of subproblem 2.18 using the conjugate gradient method ^[21], then the fraction of Cauchy decrease condition 2.19 is hold. We will consider two cases:

Firstly, if $d^{ncp} = - \frac{\delta_k }{\| P_k\nabla h_k h_k \| }(P_k \nabla h_k h_k)$ and $\|\delta_k \|(P_k\nabla h_k)^T P_k \nabla h_k h_k \|^2 \leq \|P_k \nabla h_k h_k \|^3$ then

Secondly, if $d^{ncp} = - \frac{\|P_k\nabla h_k h_k \|^2}{\|(P_k\nabla h_k)^T P_k \nabla h_k h_k \|^2}(P_k \nabla h_k h_k)$ and $\delta_k \|(P_k\nabla h_k)^T P_k \nabla h_k h_k \|^2\geq\|P_k \nabla h_k h_k \|^3$, then

Using assumption $GS_2$, we have

Then, from the above inequality, inequalities 3.5, 3.6, and using assumption $GS_3$, we have

From the above inequality and 2.19, we have

Since $0 < \gamma_k \leq 1$, then we have

From inequality 3.7 and the above inequality, we have

From inequalities 2.32 and 3.8 we have

Lemma 3.3. Under assumptions $GS_1$-$GS_5$, there exists a positive constant $K_3$, such that

Proof. Since the conjugate gradient method is used to solve subproblem 2.22 to obtain an approximate solution for $\bar{d}_k^t$, then the fraction of Cauchy decrease condition 2.23 is hold and we will consider two cases:

Firstly, if $\bar{d}_k^{tcp} = -\frac{\Delta_k }{\| Z_k^T\nabla q_k(P_k d_k^n) \| }(Z_k^T\nabla q_k(P_k d_k^n))$ and $\Delta_k(Z_k^T\nabla q_k(P_k d_k^n))^T \bar{B}_k Z_k^T\nabla q_k(P_k d_k^n)\leq \\\| Z_k^T\nabla q_k(P_k d_k^n) \|^3$, then

Secondly, if $\bar{d}_k^{tcp} = -\frac{\| Z_k^T\nabla q_k(P_k d_k^n) \|^2}{Z_k^T\nabla q_k(P_k d_k^n))^T \bar{B}_k Z_k^T\nabla q_k(P_k d_k^n)}Z_k^T\nabla q_k(P_k d_k^n)$ and $\Delta_k(Z_k^T\nabla q_k(P_k d_k^n))^T \bar{B}_k Z_k^T\nabla q_k(P_k d_k^n)\geq \| Z_k^T\nabla q_k(P_k d_k^n) \|^3$, then

From inequalities 3.10, 3.11, and using necessary assumptions, we have

From condition 2.23 and the above inequality, we have

Since $0 < \gamma_k \leq 1$, then we have

From the above inequality and inequality 3.12, we have

This completes the proof.

The following lemma is needed in many forthcoming lemmas. In what follows, we use implicitly that $\nabla h_k d_k^n = \nabla h_k d_k$.

Lemma 3.4. Under assumptions $GS_1$-$GS_5$, there exists a positive constant $K_4$, such that

Proof. Since $d_k^n$ is normal to the tangent space, then we have

Using the fact that $\|h_k+(P_k\nabla h_k)^T d_k\|\leq\| h_k\|$, we have

From above inequality and necessary assumptions, we have

Since

From the above inequality and inequality 3.14 and using the fact that $\| d_k^n \|\leq \delta_{max}$, we have

This completes the proof.

Lemma 3.5. Under assumptions $GS_1$-$GS_5$,

Proof. From Equation 2.31, we have

Using inequalities 3.9 and 3.13, we obtain the desired result.

The following lemma shows that, if $\|Z_k^T P_k \nabla_x \ell_k \| \geq \varepsilon_1$ and $\| h_k \| \leq \eta \delta_{k^i}$ at any trial iteration $k^i$, then the penalty parameter $\rho_k$ is not increased.

Lemma 3.6. Under assumptions $GS_1-GS_5$, if $\|Z_k^T P_k \nabla_x \ell_k \| \geq \varepsilon_1$ and $\| h_k \| \leq \eta \delta_{k^i}$ at any trial iteration $k^i$, then there exists a positive constant $K_5$, such that

where $\eta$ is given by

Proof. Since $\|Z_k^T P_k \nabla_x \ell_k \| \geq \varepsilon_1$ and $\| h_k \| \leq \eta \delta_{k^i}$, and using inequalities 3.1 and 3.14, we have

But $\eta \leq \frac{\varepsilon_1}{2 b_1 \kappa_4 \delta_{\max}}$, hence

From inequality 3.14, assumption $\| h_k \| \leq \eta \delta_{k^i}$, and $\eta \leq \frac{\sqrt{3}}{2\kappa_4}$, then we have $\|d_{k^i}^n\|\leq \kappa_4 \eta \delta_{k^i}\leq \kappa_4\frac{\sqrt{3}}{2\kappa_4}\delta_{k^i} = \frac{\sqrt{3}}{2}\delta_{k^i}$. Since $\Delta_{k^i} = \sqrt{\delta_{k^i}^2-\|d_{k^i}^n\|^2}$, then $\Delta_{k^i}\geq\frac{1}{2}\delta_{k^i}$. Hence, from inequalities 3.15 and 3.17, we have

But $\eta \leq \frac{ K_3 \varepsilon_1 }{ 8 K_4 }\min \{\frac{ \varepsilon_1 }{ b_1\delta_{\max} }, 1 \}$, hence

The result follows if we take $K_5 = \frac{K_3 \varepsilon_1}{8} \min \{ \frac{\varepsilon_1}{b_1 \delta_{\max}}, 1\}$.

The following lemma shows that, at any iteration $k$, we can find an acceptable step after finite number of trials, or equivalently, the condition $Ared_{k^j}/Pred_{k^j} \geq \tau_1$ will be satisfied for some finite $j$.

Lemma 3.7. Under assumptions $GS_1-GS_5$, if $\| h_k \| > \varepsilon_1$, where $\varepsilon_1 > 0$, then $Ared_{k^j}/Pred_{k^j} \geq \tau_1$ will be satisfied for some finite $j$.

Proof. From inequalities 3.2, 3.4, and assumption $\| h_k \| > \varepsilon_1$, we have

If the trial step $d_{k^j}$ gets rejected, then $\delta_{k^j}$ becomes small and hence we have

That is the acceptance rule will be met after finite number of trials (i.e., for finite $j$) and this completes the proof.

Lemma 3.8. Under assumptions $GS_1-GS_5$ and at any iteration $k$, if

at the $j^{th}$ trial step, then the step must be accepted.

Proof. Assume that inequality 3.18 holds and the step $d_{k^j}$ is rejected. From the way of updating trust-region radius which is clarified in Algorithm 2.2 we have

From the above inequality and using inequalities 3.2, 3.4, and 3.18 we have

This is a contradiction with the assumption $d_{k^j}$ was rejected. Hence the step must be accepted.

Lemma 3.9. Under assumptions $GS_1-GS_5$ and for all trail iterates $j$ of any iteration $k$ we have

where $b_2 = \sup_{x\in \Omega}\| h_k\|$.

Proof. Consider any trial iterate $k^j$, if $j = 1$, then the step is accepted and hence

such that $b_2 = \sup_{x\in \Omega}\| h_k\|$.

Now, if $j > 1$, then there exists at least one rejected trial step. For all rejected trial steps, we have from Lemma 3.8,

for all $i = 1, 2, ...j-1$. Since $d_k^i$ is rejected trial step, then from the way of updating the radius of trust-region, we have

From the above inequality and inequality 3.20, we obtain the desired results.

The following lemma show that the sequence of trust-region radii $\{ \delta_{k^j} \}$ is bounded away from zero if $\{ \| h_k \|\}$ is bounded away from zero.

Lemma 3.10. Under assumptions $GS_1-GS_5$, if $\| h_k \| \geq \varepsilon_1$ where $\varepsilon_1 > 0$, then there exists a constant $K_6 > 0$ such that

for all trial iterates $j$ of any iteration $k$.

Proof. From Lemma 3.9 and the condition $\| h_k \| \geq \varepsilon_1$, the proof follows directly by taking $K_6 = \min\{ \frac{\delta_{\min}}{b_2}, \beta_1 \frac{(1-\tau_1) K_2}{4 K_1}, \beta_1\}\varepsilon_1$.

Lemma 3.11. Under assumptions $GS_1-GS_5$, there exists a subsequence $\{ k_i \}$ of the iteration sequence at which $\rho_k$ is increased such that at any trial steps $j$ of any iteration $k\in\{ k_i \}$, we have

where $K_7$ is a positive constant.

Proof. At any trial steps $j$ of any iteration $k$, if $\rho_{k^j}$ is increased, then from equation 2.33, we have

Applying inequality 3.8 on the left hand side and inequalities 3.9, 3.13, and 3.14 on the right hand side, we have

From assumptions $GS_2$, $GS_3$, and using the fact that $\|d_{k^j}^n \|\leq \delta_{k^j}\leq \delta_{max}$, we have

From inequalities 3.19 and 3.23, there exists a constant $K_7 > 0$ such that

at any trial steps $j$ for any iteration $k\in \{k_i\}$.

In the following lemma we will prove that the sequence $\{ \|h_k \| \}$ is not bounded away from zero when $\{\rho_k\}$ unbounded sequence.

Lemma 3.12. Under assumptions $GS_1-AS_6$, there exists a subsequence $\{ k_i \}$ of the iteration sequence at which $\rho_k$ is increased such that

Proof: From Lemma 3.11 and the assumption $\rho_k$ is increased, we obtain the desired result.

In the following section, we prove the main global convergence results for IPTR algorithm 2.4.

3.3.   Fundamental convergence theorem

In the following theorem we prove that the sequence of the iterates generated by algorithm 2.4 converges to the feasible set.

Theorem 3.1. Under assumptions $GS_1-GS_5$, the sequence of iterates generated by IPTR algorithm satisfies

Proof. The proof of this theorem is by contradiction, so we assume that $\limsup_{k \rightarrow \infty} \| h_k \| \geq \varepsilon_1$ where $\varepsilon_1 > 0$. This implies the existence an infinite subsequence of indices $\{k_j\}$ indexing iterates that satisfy $\| h_k \| \geq \frac{\varepsilon_1}{2}$, for all $k \in \{ k_j \}$. From Lemma 3.7, there exists a finite sequence of acceptable steps. Without lose of generality, we assume all members of the sequence $\{ k_j \}$ are acceptable iterates. Now we will consider two cases:

Firstly, if the sequence of the penalty parameter $\{\rho_k\}$ is unbounded, then there exists a subsequence $\{ k_i \}$ of the iteration sequence at which $\rho_k$ is increased. Using Lemma 3.12, we have $\lim_{k_i \rightarrow \infty}\| h_{k_i} \| = 0$. Therefore, there are no common elements between $\{k_i\}$ and $\{k_j\}$ at iteration $k$ which is sufficiently large.

From inequality 3.4 and Lemma 3.10, we have

for all $k \in \{k_j\}$, such that $\bar{K}_{6} = \frac{ \varepsilon_1}{2} \min\{ \frac{\delta_{\min}}{b_2}, \beta_1 \frac{(1-\tau_1) K_2}{2 K_1}, \beta_1 \}$. Since

then from 3.25 we have

Hence

for all acceptable steps which are generated by IPTR algorithm 2.4. Let $k \in \{ k_j \}$ be an element between the two elements $k_{\hat{i}}$ and $k_{\hat{i}+1}$ which are consecutive elements of the sequence $\{ k_i \}$. From inequality 3.26, we have

Since the value of $\rho_k$ is the same for all iterates $k_{\hat{i}}, \ldots, {k_{\hat{i}+1}-1}$, we have

Since $\rho_k \rightarrow \infty$ as $k \rightarrow \infty$, and $|\ell_k|$ is bounded, we can write

for $k_{\hat{i}}$ sufficiently large. But this leads to a contradiction with Lemma 3.12.

Secondly, if the sequence of the penalty parameters $\{\rho_k\}$ is bounded, then there exists an integer $\bar{k}$ such that for all $k \geq \bar{k}$, we have $\rho_k = \bar{\rho}$. Since all the iterates of $\{ k_j \}$ are acceptable, then for any $\bar{k}\in \{ k_j \}$ we have

From Lemma 3.10, inequality 3.4, we have for any $\tilde{k}\in \{ k_j \}$ and $\tilde{k}\geq \bar{k}$

such that $K_8 = \frac{K_2 \tau_{\tilde{k}} \bar{\rho}\varepsilon_1}{4} \min \{ \| \frac{\varepsilon_1}{2\delta_{max}}, 1 \}$. From inequalities 3.28 and 3.29, we have

This gives a contradiction with the fact that $\{\Phi_k \}$ is bounded below when $\{\rho_k\}$ is bounded. Hence in both cases, we have a contradiction. Thus, the supposition is not correct and the theorem is proved.

Theorem 3.2. Under assumptions $GS_1-GS_5$, the algorithm is terminated because

Proof. Assume that IPTR algorithm 2.4 does not terminate and that some subsequences of $\{ \| Z_k^T P_k \nabla \ell_k\|\}$ convergence to zero, then the nontermination is immediately contradicted by Theorem 3.1.

Now assume that for $\bar{k}$ sufficiently large, there exists an index $\tilde{k} > \bar{k}$ such that $\| Z_k^T P_k \nabla \ell_k\|\geq \varepsilon_1$. Let $\{ k_j \}$ be a subsequence of iterates that satisfy $\|h_{k_j}\| > \eta \delta_{k_j}$, then $\lim_{k_j \rightarrow \infty} \delta_{k_j} = 0$ such that $\lim_{k_j \rightarrow \infty}\|h_{k_j}\| = 0$. This implies the existence of an infinite sequence $\{ k_j \}$ of rejected trial steps. But this leads to contradiction. To show this, we consider two cases:

Firstly, if the sequence of the penalty parameter $\{\rho_k\}$ is unbounded, then from inequalities 3.3 and 3.4, we have

As $\rho_{k_j} \rightarrow \infty$ and $\delta_{k_j} \rightarrow 0$, then $\frac{\left | Ared_{k_j}- Pred_{k_j} \right |}{Pred_{k_j}}\rightarrow 0$. This means that for $k_j$ large enough, all trial steps $\|d_{k_j} \|$ must be accepted. This leads to a contradiction, so $\delta_{k_j}$ must be bounded away from zero in this case.

Secondly, if the sequence of the penalty parameter $\{\rho_k\}$ is bounded, then there exists an integer $\bar{k}$ such that for all $k\geq \bar{k}$, $\rho_k = \bar{\rho}$. Now, we discuss three cases:

1] If the previous step is accepted ($j = 1$), then from the way of updating the trust-region radius in algorithm 2.2, we have $\delta_{k_j} \geq \delta_{min}$. That is $\delta_{k_j}$ is bounded away from zero in this case.

2] If $j > 1$ and $\|h_{k_r}\| > \eta \delta_{k_r}$ for all $r = 1, \cdots, j-1$. Then

such that all the trial steps on $\{ k_j \}$ are rejected. From above inequality and inequalities 3.2 and 3.4 we have

Hence

But from the way of updating the radius of trust-region in algorithm 2.2, all the rejected trial steps satisfy $\delta_{k_r} = \beta_1\|d_{k_r}\|$, hence

This means that $\delta_{k_r}$ is bounded away from zero in this case.

3] If $j > 1$ and $\|h_{k_r}\| > \eta \delta_{k_r}$ does not hold for all $r$. Hence, there exists an integer $i$ such that $\|h_{k_r}\|\leq \eta \delta_{k_r}$ for all $r = 1, \cdots, i$, and $\|h_{k_r}\| > \eta \delta_{k_r}$ for all $r = i+1, \cdots, j-1$. Since $\|h_{k_r}\| > \eta \delta_{k_r}$ for all $r = i+1, \cdots, j-1$, then as the above case we can prove $\delta_{k_r}$ is bounded away from zero.

The case when $\|h_{k_r}\|\leq \eta \delta_{k_r}$ for all $r = 1, \cdots, i$, then for all rejected trial steps, we have

From inequality 3.2, Lemma 3.6, and the above inequality, we have

Hence

From the way of updating the radius of trust-region, we have for all rejected trial step

Hence, $\delta_{k_r}$ is bounded away from zeros. This leads to a contradiction and then for $k_j$ sufficiently large, all the iterates satisfy $\|h_k\|\leq\eta \delta_{k_j}$.

For all successful steps and from the way of updating the radius of trust-region and Lemma 3.6, we have for all $k\in \{k_j\}$ and $k\geq \bar{k}$

We proved in the above cases, that $\delta_{k_j}$ is bounded away from zeros. Then $\Phi_k - \Phi_{k+1} > 0$. This leads to a contradiction with the fact that $\{\Phi_k \}$ is bounded below when $\{\rho_k\}$ is bounded. Hence in both cases, we have a contradiction. Thus, the supposition is not correct and the theorem is proved.

4.   Application

In this section, firstly the proposed algorithm IPTR is applied to the engineering application which is called two-echelon supply chain system with one manufacturer and one retailer.

The manufacturer purchases raw materials from the supplier first, then after the manufacturer's production and processing, the end products are sold to the retailer, this problem is formulated as bilevel models for joint pricing and lot-sizing decisions, see ^[34].

where $\tilde{T} = 52$; $\tilde{P}_s = 4$; $\tilde{T}_c = 0.5$; $\tilde{M}_c = 1$; $\tilde{c}_m = \tilde{c}_r = 0.001$; $\tilde{O}_m = 400$; $\tilde{O}_r = 200$. For more details about the above application and its notations, see ^[34].

We solve this model in case of the manufacturer is the leader, who makes the first decision, and the retailer is the follower. Our results, when applying Algorithm (2.4) is $t_1 = 5.8778$, $t_2 = 6.002$, $t_3 = 19710.195$, $y_1 = 7.691$, $y_2 = 2.6007$, $f_u = 431230$, and $f_l = 8548300$, which is closed to whose reported in ^[34].

Secondly, we introduce an extensive variety of possible numeric bilevel nonlinear programming problems to clarify the effectiveness of our IPTR algorithm, since, Problems 1, 2, 6, 7, 13, and 14 have quadratic functions in both levels. Problems 3, 4, 5, 8, 9 all the inner level functions are convex and Problem 10 ^[27], at fixed $x$, the inner problem is convex. These problems are solved numerically with the help of algorithm (2.4) to clarify the effectiveness of that approach. For each test example, 10 independent runs with different initial starting point are performed to observe the consistency of the outcome. Statistical results of all examples are summarized in Table 1 which shows that the results found by the IPTR algorithm (2.4) are approximate or equal to those by the compared algorithms in the literature.

Table 1 also including the mean number of iterations (iter), the mean number of function evaluations (nfunc), the mean value of CPU time (CPUs) in seconds.

For comparison, we have included the corresponding results of the mean value of CPU time (CPUs) obtained by Method in ^[31](Table 2), ^[27](Table 3), and ^[44](Table 4) respectively. It is clear from the results that our approach is capable for treating nonlinear bilevel programming problems even the upper and the lower levels are convex or not and the computed results converge to the optimal solution which is similarly or approximate to the optimal that reported in literature. Finally, it is clear from the comparison between the solutions obtained using IPTR algorithm with literature, that IPTR is able to find the optimal solution of all problems by a small number of iterations, small number of function evaluations, and less time.

Problem 1 ^[31]:

Problem 2 ^[31]:

Problem 3 ^[31]:

Problem 4 ^[31]:

We offered the numerical results of our algorithm using MATLAB (R2013a)(8.2.0.701)64-bit(win64) and a starting point $x_0\in int (\hat{\boldsymbol{G}})$. The following parameter setting is used: $\delta_{min} = 10^{-3}$, $\delta_0 = max(\| s_0^{cp} \|, \delta_{min})$, $\delta_{max} = 10^3\delta_0$, $\tau_1 = 10^{-4}$, $\tau_2 = 0.75$, $\beta_1 = 0.5$, $\beta_2 = 2$, $\hat{\varepsilon} = 0.01$, $\varepsilon_1 = 10^{-8}$, and $\varepsilon_2 = 10^{-10}$.

Problem 5 ^[31]:

Problem 6 ^[31]:

Problem 7 ^[31]:

Problem 8 ^[31]:

Problem 9 ^[27]:

Problem 10 ^[27]:

Problem 11 ^[44]:

Problem 12 ^[27]:

Problem 13 ^[44]:

Problem 14 ^[44]:

Problem 15 ^[44]:

Problem 16 ^[44]:

5.   Concluding remarks

This paper presented a new technique for solving a nonlinear bilevel optimization problem based on using the slack variable with KKT condition to transform NBLP problem into an equivalent smooth SONP problem. A Newton's interior-point method with Das scaling matrix is utilized to solve the equivalent smooth SONP problem effectively. Newton's method is locally method, so a trust region technique is utilized to ensure global convergence from any starting point. On applying this methodology we overcome some known difficulties on treating such problems, as

● A trust-region technique can induce strongly global convergence, which is very important technique for solving a smooth optimization problems and is more robust when they deal with rounding errors

● Our approach used to transform Problem 1.3 which is not smooth to smooth problem

● Using the interior-point method guarantees the converges quadratically to a stationary point.

On the other hand, the global convergence theorems for the IPTR algorithm is presented and numerical results reflect the good behavior of our algorithm and computed results converge to the optimal solutions. Finally, it is clear from the comparison between the solutions obtained using IPTR algorithm with literature, that IPTR is able to find the optimal solution of all problems by a small number of iterations.

Acknowledgments

The authors would like to thank the anonymous referees for their valuable comments and suggestions which have helped to greatly improve this paper.

Conflict of interest

The authors declare that there is no conflict of interest in this paper.