The orthogonal reflection method for the numerical solution of linear systems

1.   Introduction

A nested optimization problem with two levels in a hierarchy, i.e the upper-level and lower-level decision-making, is known as the bilevel programming problem. Each level has distinct constraints and objective functions. Both of them have their objective functions and constraints. The decision-maker at the upper level always takes the lead, followed by those at the lower level. The objective function and constraint of the upper-level programming depend on their decision variables and the optimum solution of the lower-level programming. The decision maker at the lower level must maximize its objective function by using the variables provided by the decision maker at the upper level, who in turn chooses the variables after having full knowledge of the lower level's potential responses. Bilevel mathematical programming stresses the system's non-cooperative nature, in contrast to many objective mathematical programming methodologies. This hierarchical model has several applications, including resource allocation, decentralized control, and network design issues. The successful application of this hierarchical model depends on how well it is solved when handling realistic complications, etc. ^[3]. Bilevel mathematical programming has attracted a lot of attention, and many effective algorithms have been presented. There are currently several methods for solving NBLP problems and they can be divided into four categories: the Karush-Kuhn-Tucker condition approach ^{[1,15,17,19,21]}, the penalty function approach ^[24,29], the descent approach ^[22,36], and the evolutionary approach ^[41].

In this paper, a class of NBLP problems is reduced to a traditional nonlinear programming problem with complementary constraints. The lower-level problem is then substituted by its Karush-Kuhn-Tucker optimality conditions. The CHKS smoothing function is then used to smooth it.

The smoothed nonlinear programming problem (SNLP) is solved by using the nonmonton active interior-point trust-region technique to obtain an approximately optimal solution to the nonlinear bilevel programming problem. In the nonmonton active interior-point trust-region technique, the smoothed nonlinear programming problem is transformed into an equality-constrained optimization problem with bounded variables by using an active-set approach and the penalty method; for more details see ^{[12,13,14,15,18,19]}. To solve the equality- constrained optimization problem with bounded variables, Newton's interior-point approach ^[6] is used. Because Newton's interior-point method is a local method, it might not converge if the starting point is far from a stationary point. A trust-region strategy is used to treat this problem and ensure convergence to the stationary point from any starting point. A trust-region technique can induce strong global convergence and it is a very important method for solving unconstrained and constrained optimization problems; see ^{[10,11,12,13,14,16,17,18,19]}. One advantage of the trust-region technique is that it does not require the models objective function to be convex. However, in the traditional trust-region strategy, we must use some criteria to determine whether the trial step is acceptable after solving a trust-region subproblem. Otherwise, the trust-region radius needs to be reduced. A method for calculating the trust-region radius $\Delta_k$ at each iteration is an important part of trust-region techniques. The standard trust-region strategy is predicated on the objective function and the model agreement. The trust region's radius is updated by paying attention to the ratio $t_k = \frac{Ared_k}{Pred_k}$ where $Ared_k$ refers to the actual reduction and $Pred_k$ refers to the predicted reduction. It can be deduced that whenever $t_k$ is close to 1, there will be a good agreement between the model and the objective function over a current trust region. It is well-known that the standard trust-region radius $\Delta_k$ is independent of the gradient and Hessian of the objective function, so we are not able to know if the radius $\Delta_k$ is convenient for the whole implementation. This condition might lead to an increase in the number of subproblems that must be resolved in the method's inner phases which reduces the method's efficiency.

To overcome this problem many authors proposed various nonmonotone trust-region methods, for example, see ^{[9,31,37,38,43,44]}. Motivated by the nonmonotone trust-region strategy in ^[31], we use it in our method. It is generally promising and efficient, and it can overcome the aforementioned shortcomings.

Furthermore, the usefulness of the CHKS smoothing function with the nonmonton active interior-point trust-region algorithm to solve the NBLP problems was examined by using several benchmark problems and a real-world case about a watershed trading decision-making problem. Numerical experiments show that the suggested method surpasses rival algorithms in terms of efficacy.

This paper is organized as follows: Section 2 introduces the mathematical formulation of the NBLP problem, basic definitions of the CHKS smoothing functions, and the smoothing method for the nonlinear complementarity problem to obtain the SNLP problem. Section 3 introduces the nonmonton active interior-point trust region algorithm for solving the SNLP problem. Results from simulations on several benchmark problems and a real-world case about a watershed trading decision-making problem are reported in Section 4. The conclusion is given in Section 5.

2.   Mathematical model for SNLP problem

In this paper, we will consider the following NBLP problem

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

The NBLP problem (2.1) was reduced to the following one-objective optimization problem using Karush-Kuhn-Tucker optimality assumptions for the lower level problem:

where $\lambda\in \Re^{m_2}$ is a multiplier vector associated with the inequality constraint $c_l(y, z)$. Problem (2.2) with the nonlinear complementarity condition is non-convex and non-differentiable; moreover, the regularity assumptions required to handle smooth optimization problems are never satisfied and it is not good to use our approach to solve problem (2.2). Due to this, we use the CHKS smoothing function to overcome this problem, see ^[5,26,35].

Definition 2.1. The Chen-Mangasarian smoothing function is represented by the notation $\hat{\phi}(g, h): \Re^2 \rightarrow \Re$ and defined by the expression $\hat{\phi}(g, h) = g+h-\sqrt{(g-h)^2}$. By introducing a smoothing parameter $\tilde{\epsilon}\in R$ into the the smoothing function $\hat{\phi}(g, h)$, we obtain the CHKS smoothing function

The Chen-Mangasarian smoothing function has the property $\hat{\phi}(g, h) = 0$ if and only if $g \geq 0$, $h \geq 0$, $gh = 0$ but it is non-differentiable at $g = h = 0$. But, the CHKS smoothing function has the property $\hat{\phi}(g, h, \tilde{\epsilon}) = 0$ if and only if $g\geq 0$, $h \geq 0$, and $gh = \frac{\tilde{\epsilon}}{2}$ for $\tilde{\epsilon} \geq 0$, and the function is smoothing for g, h, and $\tilde{\epsilon} \geq 0$.

Consequently, by using the CHKS smoothing function (2.3), the problem (2.2) can be approximated as follows:

The above smoothed nonlinear programming problem can be summarized as follows

where $x = (y, z, \lambda)^T\in \Re^n$ where $n = {n_1+n_2+m_2}$, $c_e(x) = [\nabla_{z} f_l(y, z)+ \nabla_{z} c_l(y, z) \lambda, \lambda_j-c_{l_j}-\sqrt{(\lambda_j+c_{l_j})^2+4 \tilde{\epsilon}^2}]$, $j = 1, ..., m_2$. The functions $f_u(x): \Re^n \rightarrow \Re$, $c_u (x) : \Re^n \rightarrow \Re^{m_1}$, and $c_e (x) : \Re^n \rightarrow \Re^{n_2+m_2}$ are twice continuously differentiable and $m_1 < n$. We denote the feasible set ${\boldsymbol{{E}}} = \{x: \beta_a \leq x \leq \beta_b \}$ and the strict interior feasible set $int ({\boldsymbol{{E}}}) = \{x : \beta_a < x < \beta_b \}$ where $\beta_a \in\{\Re\bigcup\{- \infty\} \}^n$, $\beta_b \in\{\Re\bigcup\{\infty\} \}^n$, and $\beta_a < \beta_b$.

Several methods that have been suggested to solve the smoothed nonlinear programming problem (2.5), see ^{[11,12,16,17,19]}. The nonmonton active interior-point trust-region algorithm is proposed in this paper to solve problem (2.5) and a detailed description of this algorithm is clarified in the following section.

3.   Nonmontone active-set interior-point trust-region algorithm

In this section, firstly, we will offer a detailed description of the active-set strategy with the penalty method to convert problem (2.5) to an equality-constrained optimization problem with bounded variables. Second, the basic steps for using Newton's interior-point method to solve the equality-constrained optimization problem are presented clearly. Thirdly, the main steps for the nonmontone trust-region algorithm are presented. Finally, the main steps for solving problem 2.1 are introduced.

3.1.   An active-set strategy

Motivated by the active-set strategy proposed in ^[8] and used in ^{[10,11,12,13,14]}, we define a 0-1 diagonal matrix $D(x) \in \Re^{m_1\times m_1}$, whose diagonal entries are defined as follows:

Problem (2.5) is transformed into the following equality-constrained optimization problem using the previous matrix.

The previous problem is transformed into the following problem by using a penalty method

where $\sigma_u$ is a positive parameter.

3.2.   Newton's interior-point method

Motivated by the interior point method in ^[6], we let $L(x, \mu_e, \lambda_a, \lambda_b)$ be a Lagrangian function associated with problem (3.2) and it is defined as follows

where

and

such that $\mu_e$, $\lambda_a$, and $\lambda_b$ represent the Lagrange multiplier vectors associated with the equality constraint $c_e(x)$ and inequality constraints $(x-\beta_a)$ and $(\beta_b- x)$ respectively. A point $x_*\in {\boldsymbol{E}}$ will be a local minimizer of problem (3.2) if there exists multiplier vectors $\mu_{e_*} \in \Re^{m_1}$, $\lambda_{a_*} \in \Re_{+}^n$, and $\lambda_{b_*} \in \Re_{+}^n$ such that $(x_*, \mu_{e_*}, \lambda_{a_*}, \lambda_{b_*})$ satisfies the following Karush-Kuhn-Tucker conditions,

where

and

Let $V(x)$ be a diagonal matrix whose diagonal elements are as follows:

For more details see ^[11,12,20]. Using the scaling matrix $V(x)$, conditions (3.6)–(3.9) are equivalent to the following nonlinear system,

For the following reasons, the nonlinear systems (3.13) and (3.14) is continuous but not everywhere differentiable.

● It may be non-differentiable when $v^{(j)} = 0$. To overcome this problem, restricting $x \in int({\boldsymbol{E}})$.

● It may be non-differentiable when $v^{(j)}$ has an infinite upper bound and a finite lower bound, and $(\nabla_x \ell(x_*, \mu_{e_*}; \sigma_{u_*}))^{(j)} = 0$. To overcome this problem, we define a vector $\theta(x)$ whose components $\theta^{(j)}(x) = \frac{\partial({(v^{(j)})}^2)}{\partial x^{(j)}}$, $j = 1, ..., n$ are defined as follows

If we use Newton's method to solve the nonlinear systems (3.13) and (3.14), we get

where

and $H$ is the Hessian of the Lagrangian function (3.4) or an approximation to it.

Restricting $x\in int ({\boldsymbol{E}})$ is necessary to ensure that the matrix $V(x)$ is nonsingular. Therefore, putting $\Delta x = V(x)s$ in both Eqs (3.16) and (3.17) and multiplying both sides of equation (3.16) by $V^{-1}(x)$, we get

It should be noted that the step $s_k$ produced by the systems (3.19) and (3.20) is equivalent to the step produced by resolving the following quadratic programming subproblem,

where

and

While Newton's approach has the advantage of being quadratically convergent under reasonable assumptions, it also has the disadvantage of requiring that the starting point be close to the solution. The nonmonotone trust-region globalization approach is used to ensure convergence from any starting point. The nonmonotone trust-region globalization strategy is a crucial method for solving a smooth nonlinear unconstrained or constrained optimization problem that can produce substantial global convergence. In the following section, we present the main steps of the nonmontone trust-region algorithm to solve quadratic subproblem (3.21).

The main steps of the process to apply the nonmontone trust-region algorithm to solve the quadratic subproblem (3.21) are described in the section that follows.

3.3.   The nonmontone trust region algorithm

The trust-region subproblem associated with the problem (3.21) is given by

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

Due to the possibility of no intersecting points existing between the hyperplane of the linearized constraints $c_e(x) + V(x)\nabla c_e(x)^T s = 0$ and the constraint $\|s\|\leq \delta_k$, subproblem (3.26) may be infeasible. There is no guarantee that this will remain true even if they do intersect if $\delta_k$ is reduced; see ^[7]. To solve this problem, we used a reduced Hessian strategy. This strategy was suggested in ^[2,32] and subsequently implemented in ^{[12,13,16,17,18,20]}. This strategy divides the step $s_k$ into two orthogonal components: the normal component $s_k^n$ for improving feasibility and the tangential component $s_k^t$ for improving optimality. That is, $s_k = s_k^n + s_k^t$ and $s_k^t = \tilde{Z}_k \bar{s}_k^t$ where $\tilde{Z}_k$ is a matrix whose columns form a basis for the null space of $(V_k\nabla c_{e_k})^T$.

● To compute $s_k^n$

Obtaining the normal component $s_k^n$ by solving the following trust-region subproblem

for some $\zeta\in (0, 1)$.

A conjugate gradient method ^[34] which is very cheap if the problem is large-scale and the Hessian is indefinite, is used to compute the normal component $s_k^n$. The main steps involved in applying the conjugate gradient method to solve subproblem (3.25) are presented in the following algorithm.

Algorithm 3.1. : (A conjugate gradient method to calculate $s_k^n$)

Step 1. Set $s_0^n = 0$, $r_{n_0} = -V_k\nabla c_{e_k} c_{e_k}$ and $p_{n_0} = r_{n_0}$; pick $\epsilon > 0$.

Step 2. For i = 0, 1, .... do

Compute $\; \gamma_{n_i} = \frac{ {r_{n_i}}^T r_{n_i}}{{p_{n_i}}^T V_k \nabla c_{e_k} \nabla c_{e_k}^TV_k p_{n_i} }$.

Compute $\tau_{n_i}$ such that $\|s_i^n+\tau_n p_{n_i}\| = \delta_k$.

If $\gamma_{n_i}\leq 0$, or if $\gamma_{n_i} > \tau_{n_i}$, then set $s_k^n = s_i^n+\tau_{n_i} p_{n_i}$ and stop.

Otherwise set $s^n_{i+1} = s_i^n+ \gamma_{n_i} p_{n_i}$ and

$r^n_{i+1} = r_{n_i}- \gamma_{n_i} V_k \nabla c_{e_k} \nabla c_{e_k}^T V_k p_{n_i}$.

If $\frac{\| r_{n_{i+1}}\|}{\| r_{n_0}\|}\leq\; \epsilon_0\;$,

set $s_k^n = s^n_{i+1}$ and stop.

Compute $\; \hat \beta_i = \frac{ {r_{n_{i+1}}}^T r_{n_{i+1}}}{ {r_{n_i}}^T r_{n_i}}\;$, and the new direction $p_{n_{i+1}} = r_{n_{i+1}}+\hat \beta_i p_{n_i}$.

Given $s_k^n$, let $q(V_k s_k)$ be the quadratic form of the function (3.5) and let it be defined as follows

● To compute $s_k^t$

To compute the tangential component $s_k^t = \tilde{Z}_k \bar{s}_k^t$, solve the following trust-region subproblem

where $\nabla q_k(V_k s_k^n) = V_k \nabla_x \ell(x_k, \mu_{e_k}; \sigma_{u_k})+ B_k s_k^n$ and $\Delta_k = \sqrt{\delta_k^2 - \| s_k^n \|^2}$. The main steps involved in applying the conjugate gradient method to solve subproblem (3.27) are presented in the following algorithm.

Algorithm 3.2. (A conjugate gradient method to compute $s_k^t$)

Step 1. Set $\bar{s}_0^t = s_k^n$, $r_{t_0} = -\tilde{Z}_k^T \nabla q_k(V_k s_k^n) + B_k s_k^n$, $p_{t_0} = r_{t_0}$.

Step 2. For i = 0, 1, .... do

Compute $\gamma_{t_i} = \frac{ r_{t_i}^T r_{t_i}}{p_{t_i}^T \tilde{Z}_k^T B_k \tilde{Z}_k p_{t_i}}$.

Compute $\tau_{t_i}$ such that $\|\bar{s}_i^t+\tau_t\; p_{t_i}\| = \Delta_k\;$.

If $\gamma_{t_i}\leq 0$, or if $\gamma_{t_i} > \tau_{t_i}$, then set $\bar{s}_k^t = \bar{s}_i^t+\tau_{t_i} p_{t_i}$ and stop.

Otherwise set $\bar{s}_{i+1}^t = \bar{s}_i^t+ \gamma_{t_i} p_{t_i}$ and

$r_{t_{i+1}} = r_{t_i}-\gamma_{t_i} \tilde{Z}_k^T B_k \tilde{Z}_k p_{t_i}$.

If $\; \frac{\| r_{t_{i+1}}\|}{\| r_{t_0}\|}\leq\; \epsilon\;$,

set $\bar{s}_k^t = \bar{s}_{i+1}^t$ and stop.

Compute $\bar{\beta}_i = \frac{r_{t_{i+1}}^T r_{t_{i+1}}}{r_{t_i}^T r_{t_i}}\;$, and $p_{t_{i+1}} = r_{t_{i+1}}+\bar{\beta}_i p_{t_i}$.

After computing $s_k$, we set $x_{k+1} = x_k+V_k s_k$. To ensure that the matrix $V_k$ is nonsingular, we need to guarantee that $x_{k+1}\in int {\boldsymbol{E}}$. So, the damping parameter $\varphi_k$ is required at every iteration $k$.

● To obtain the damping parameter $\varphi_k$

The following algorithm clarifies the fundamental steps required to obtain the damping parameter $\varphi_k$.

Algorithm 3.3. (Damping parameter $\varphi_k$)

Step 1. Compute the parameter $\omega_k$ as follows

Step 2. Compute the parameter $\upsilon_k$ as follows

Step 3. Compute the damping parameter $\varphi_k$ as follows

Step 4. Set $x_{k+1} = x_k + \varphi_k V_k s_k$.

To determine whether the scaled step $\varphi_k V_k s_k$ will be accepted or no, we need a merit function that connects the objective function and the constraints such that progress in the merit function equates to progress in solving the problem. The augmented Lagrangian function that follows is used as a merit function,

where $\sigma_e$ is the penalty parameter.

To test the scaled step, we use the following scheme to determine the Lagrange multiplier vector $\mu_{e_{k+1}}$,

The following actual reduction and the predicted reduction must be defined to determine if the point $(x_{k+1}, \mu_{e_{k+1}})$ will be accepted in the next iteration or no. The actual reduction in the merit function in moving from $(x_k, \mu_{e_k})$ to $(x_k+\varphi_k V_ks_k, \mu_{e_{k+1}})$ is defined as

The predicted reduction $Pred_k$ in the merit function (3.29) is defined as follows

where

To ensure that $Pred_k \geq 0$, the penalty parameter $\sigma_{e_k}$ must be updated.

● To update the penalty parameter $\sigma_{e_k}$

To ensure that $Pred_k \geq 0$, we need to update the penalty parameter $\sigma_{e_k}$ by using the following algorithm.

Algorithm 3.4. (Process to update $\sigma_{e_k}$)

If

set

where $\tilde{b}_0 > 0$ is a small fixed constant.

Else, set $\sigma_{e_{k+1}} = \sigma_{e_k}$.

End if.

For more details, see ^[12].

● To test $\varphi_k V_ks_k$ and update $\delta_k$

Motivated by the nonmonotone trust-region strategy in ^[31], we define

where

and

such that $0\leq \xi_{\min}\leq \xi_{k-1}\leq \xi_{\max}\leq 1$ and

The procedure below introduces the trial step $\varphi_k V_k s_k$ for testing and updates the radius $\delta_k$.

Algorithm 3.5. (Test $\varphi_k V_k s_k$ and update $\delta_k$)

Step 0. Choose $0 < \alpha_1 < \alpha_2 \leq 1$, $\delta_{max} > \delta_{min}$, and $0 < \zeta_1 < 1 < \zeta_2$.

Step 1. (Compute $t_k$)

Compute $\xi_k$ by using (3.38).

Compute $Q_k$ as defined in (3.37).

Compute $C_k$ as defined in (3.36).

Compute $t_k = \frac{C_k-\Phi(x_k + \varphi_k V_k s_k, \mu_{e_{k+1}}; \sigma_{u_k}; \sigma_{e_k})}{Pred_k}$.

Step 2. (Update the trial step $s_k$)

While $t_k < \alpha_1$, or $Pred_k \leq 0$,

set $\delta_k = \zeta_1\| s_k \|$.

To evaluate a new step $s_k$, go to Algorithms (3.1) and (3.2).

Step 3. (Update $\delta$)

If $\alpha_1 \leq t_k < \alpha_2$, then

set $x_{k+1} = x_k + \varphi_k V_k s_k$.

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

End if.

If $t_k \geq \alpha_2$, then

set $x_{k+1} = x_k + \varphi_k V_k s_k$.

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

End if.

● To update the parameter $\sigma_{u_k}$

To update $\sigma_{u_k}$, we use a scheme proposed in ^[39]. Let a tangential predicted decrease $Tpred$ which is obtained by the tangential component $s_k^t$ be given by

the main steps used to update the parameter $\sigma_{u_k}$ is clarified by the following algorithm.

Algorithm 3.6. (Updating $\sigma_{u_k}$)

Step 1. Set $\sigma_{u_0} = 1$.

Step 2. If

then $\sigma_{u_{k+1}} = \sigma_{u_k}$,

else, $\sigma_{u_{k+1}} = 2\sigma_{u_k}$.

End if

Finally, the nonmontone trust-region algorithm is terminated when either

or $\|s_k\|\leq \varepsilon_2$, for some $\varepsilon_1, \varepsilon_2 > 0$.

● Nonmonotone trust-region algorithm

In the algorithm that follows, we will outline the main steps of the nonmonotone trust-region algorithm in order to solve subproblem (3.21).

Algorithm 3.7. (The nonmonotone trust-region algorithm)

Step 0. Start with $x_0\in int {\boldsymbol{E}}$. Compute $D_0$, $V_0$, $\mu_{e_0}$, and $\theta_0$. Set $\sigma_{u_0} = 1$, $\sigma_{e_0} = 1$, and $\tilde{b}_0 = 0.1$.

Choose $\varepsilon_1 > 0$ and $\varepsilon_2 > 0$. Choose $\delta_{min}$, $\delta_{max}$, and $\delta_0$ such that $\delta_{min}\leq \delta_0 \leq \delta_{max}$.

Choose $\alpha_1$, $\alpha_2$, $\zeta_1$, and $\zeta_2$ such that $0 < \zeta_1 < 1 < \zeta_2$, and $0 < \alpha_1 < \alpha_2 < 1$. Set $k = 0$.

Step 1. (Termination)

If $\| \tilde{Z}_k^T V_k \nabla_x \ell(x_k, \mu_{e_k}) \| + \|V_k \nabla c_{u_k} D_k c_{u_k} \| + \| c_{e_k} \| \leq \varepsilon_1$, then stop the algorithm.

Step 2. (Computing step $V_k \varphi_k s_k$)

Evaluate the normal component step $s_k^n$ by using Algorithm (3.1).

Evaluate the tangential step $\bar{s}_k^t$ by using Algorithm (3.2).

Set $s_k = s_k^n + \tilde{Z}_k \bar{s}_k^t$.

If $\| s_k \| \leq \varepsilon_2$, then stop the algorithm.

Else, compute the parameter $\varphi _k$ by using (3.28).

Set $x_{k+1} = x_k + V_k \varphi_k s_k$.

End if.

Step 3. Evaluate $D_{k+1}$ which is defined by (3.1).

Step 4. (Updating the multiplier vector $\mu_{e_{k+1}}$)

Computing the Lagrange multiplier vector $\mu_{e_{k+1}}$ by using (3.30).

Step 5. (Updating the penalty parameter $\sigma_e$)

Using Algorithm (3.4) to update the penalty parameter $\sigma_{e_k}$.

Step 6. Compute $t_k$ as defined in (3.35) by using both $Q_k$ as defined in (3.37) and $C_k$ as defined in (3.36).

Using Algorithm (3.5) to test the scaling step $\varphi_k V_k s_k$ and update the radius of the trust-region.

Step 7. Updating the parameter $\sigma_{u_k}$ using Algorithm (3.6).

Step 8. Utilize (3.12) to evaluate the matrix $V_{k+1}$.

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

The global convergence theory for Algorithm (3.7) is similar to the global convergence theory of the algorithm presented in (^[12]) when used to solve a general nonlinear programming problem.

3.4.   CHKS nonmontone active-set interior-point trust-region algorithm

We will outline the main steps of the nonmontone active-set interior-point trust-region algorithm based on the CHKS smoothing function to resolve problem (2.1) in the following algorithm.

Algorithm 3.8. (CHKS nonmontone active-set interior-point trust-region algorithm)

Step 1. Reduce the nonlinear bilevel programming problem (2.1) to the one-objective optimization problem (2.2) using Karush-Kuhn-Tucker optimality assumptions for the lower level problem.

Step 2. Use the CHKS smoothing function to convert non-deferential problem (2.2) to smoothing problem (2.4).

Step 3. Utilize an active-set strategy to convert the nonlinearly constrained optimization problem (2.4) to the equality constrained optimization problem with bounded variables (3.2).

Step 4. Utilize an interior-point method with the diagonal scaling matrix $V(x)$ given in (3.12) to obtain the nonlinear system [(3.13) - (3.14)].

Step 5. Utilize Newton's method with diagonal matrix $\theta$ as defined in (3.15), to solve the nonlinear system [(3.13)-(3.14)] and obtain the equivalent subproblem (3.21).

Step 6. Solve subproblem (3.21) by using nonmonton trust-region given by Algorithm (3.7).

4.   Computational tests and comparisons

In this section, we introduce an extensive variety of possible numeric NBLP problems to illustrate the validity of the proposed Algorithm (3.8) to solve problem (2.1).

4.1.   Numerical examples

We shall introduced the MATLAB numerical results for Algorithm (3.8) with a starting point $x_0\in int ({\boldsymbol{{E}}})$. The parameter setting that follows was utilized: $\alpha_1 = 10^{-2}$, $\alpha_2 = 0.75$, $\zeta_1 = 0.5$, $\zeta_2 = 2$, $\delta_{min} = 10^{-3}$, $\delta_0 = max(\| s_0^{cp} \|, \delta_{min})$, $\delta_{max} = 10^3\delta_0$, $\varepsilon_1 = 10^{-10}$, and $\varepsilon_2 = 10^{-12}$.

To demonstrate the efficacy of the suggested Algorithm (3.8) in obtaining the solution to NBLP problem (2.1), we offer a wide variety of possible numeric NBLP problems in this section. Fifteen benchmark instances from ^{[4,22,27,30,33]} were used to test the proposed Algorithm (3.8).

To check the consistency of the results, 10 independent runs using different initial starting points were carried out for each test example. Table 1 summarizes the statistical data for all examples and demonstrates that the proposed Algorithm (3.8) results are approximate or equal to those of the compared algorithms in the method proposed in ^[18] and the literature.

For purposes of comparison, in Table 2 we have provided the corresponding results of the average number of iterations (iter), the average number of function evaluations (nfunc), and the average amount of CPU time (CPUs) in seconds as obtained via the method in ^[18]. The results in Table 2 demonstrate that the proposed Algorithm (3.8) results are approximative or the best of those obtained via the comparative algorithms in the literature. Results show that our proposed method is capable of handling the NBLP problem (2.1) whether the upper or lower levels are convex or not, and the calculated results converge to the optimal solution that is approximate or equal to the optimal solution stated in the literature. Finally, it is evident from a comparison of the solutions obtained via the proposed Algorithm (3.8) and those found in the literature that the proposed Algorithm (3.8) is capable of finding the best solution to some problems under conditions of a small amount of time, fewer function evaluations, and fewer iterations.

Test Problem 1 ^[28]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 1 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y_1^*, z_1^*, z_2^*) = (0.8438, 0.7657, 0)$, $f_u = -2.0769$, and $f_l = -0.5863$.

Test Problem 2 ^[28]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 2 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y_1^*, y_2^*, z_1^*, z_2^*, z_3^*) = (0.6111, 0.3890, 0, 0, 1.8339)$, $f_u = 0.64013$, and $f_l = 1.6816$.

Test Problem 3 ^[28]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 3 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y_1^*, y_2^*, z_1^*, z_2^*) = (0.97, 3.14, 2.6, 1.8)$, $f_u = -8.92$, and $f_l = -6.05$.

Test Problem 4 ^[28]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 4 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y_1^*, y_2^*, z_1^*, z_2^*) = (0.5, 0.5, 0.5, 0.5)$, $f_u = -1$, and $f_l = 0$.

Test Problem 5 ^[28]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 5 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y^*, z^*) = (9.9998, 9.9998)$, $f_u = 99.9996$, and $f_l = 3.7160e-07$.

Test Problem 6 ^[28]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 6 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y_1^*, z_1^*, z_2^*) = (1.8884, 0.89041, 0)$, $f_u = -1.2067$, and $f_l = 7.6062$.

Test Problem 7 ^[28]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 7 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y^*, z^*) = (1, 0)$, $f_u = 17$, and $f_l = 1$.

Test Problem 8 ^[28]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 8 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y_1^*, y_2^*, z_1^*, z_2^*) = (0.75;0.75;0.75;0.75)$, $f_u = -2.25$, and $f_l = 0$.

Test Problem 9 ^[23]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 9 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y^*, z^*) = (11.25, 5)$, $f_u = 2250$, and $f_l = 197.753$.

Test Problem 10 ^[23]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 10 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y_1^*, z_1^*, z_2^*) = (1, 0, 0)$, $f_u = 0$, and $f_l = 0$.

Test Problem 11 ^[40]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 11 is converted to the following smooth nonlinear programming problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y_1^*, y_2^*, z_1^*, z_2^*) = (25, 30, 4.9996, 10)$, $f_u = 5.0024$, and $f_l = 1.5000e-06$.

Test Problem 12 ^[23]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 12 is converted to the following SNLP problem

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y^*, z^*) = (3, 5)$, $f_u = 9$, and $f_l = 0$.

Test Problem 13 ^[40]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 13 is converted to the following SNLP problem.

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y_1^*, y_2^*, z_1^*, z_2^*) = (0, 2, 1.875, 0.90625)$, $f_u = -18.679$, and $f_l = -1.0156$.

Test Problem 14 ^[40]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 14 is converted to the following SNLP problem.

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y^*, z^*) = (10.016, 0.81967)$, $f_u = 81.328$, and $f_l = -0.33593$.

Test Problem 15 ^[40]:

Applying the Karush-Kuhn-Tucker condition to the lower level problem and using the CHKS smoothing function, Test Problem 15 is converted to the following SNLP problem,

Applying Algorithm (3.7) to the above nonlinear programming problem we have that $(y_1^*, y_2^*, z_1^*, z_2^*, z_3^*) = (0, 0.9, 0, 0.6, 0.4)$, $f_u = -29.2$, and $f_l = 3.2$.

4.2.   Practical example

In this section, the efficacy of the proposed Algorithm (3.8) was tested by using a watershed water trading decision-making model based on bilevel programming ^[42]. The upper decision-maker is the watershed management agency which acts as the planning, controlling and coordinating center of a watershed, and each user is the lower decision-maker. The mathematical formulation for the watershed water trading decision-making model is formulated as follows:

where $q_1$ and $q_2$ are the actual water intake volumes of water consumer A and water consumer B respectively. $l_1$ and $l_2$ are the wastewater discharge volumes of two users respectively. $r_1$ and $r_2$ are the water rights of the two users respectively. $g_1$ and $g_2$ are the emission rights of two users respectively. $w$ is the ecological water requirement of the watershed. $t$ represents the water resource fees. Applying the Karush-Kuhn-Tucker condition to the lower-level problem and using the CHKS smoothing function, the watershed water trading decision-making model is converted to the following SNLP problem,

Applying Algorithm (3.7) to the above nonlinear programming problem we have the following data at $\tilde{\epsilon} = 0.001$. $g_1 = 9$, $g_2 = 10.450$, $l_1 = 6.8755$, $l_2 = 9.0751$, $q_1 = 41.497$, $q_2 = 43$, $r_1 = 39.332$, $r_2 = 44.269$, $t = 0.3$, $w = 6.2338$, $f_1 = 12.17$, $f_2 = 19.043$, and $F = 59.055$.

Table 3 compares the obtained results with method ^[25] and those of various existing methods(Ref). From Table 3, we can see that the solutions obtained via Algorithm (3.8) are better than those given in ^[25] and those of various existing methods(Ref). Therefore, Algorithm (3.8) can also be applied to the practical problem.

5.   Conclusions

In this paper, Algorithm (3.8) has been presented to solve the NBLP problem (2.1). A Karush-Kuhn-Tucker condition is used with the CHKS smoothing function to convert the NBLP problem (2.1) into a standard SNLP problem. The SNLP problem is solved by using the proposed nonmonton active interior-point trust-region technique to obtain an approximately optimal solution to the NBLP problem. In the proposed nonmonton active interior-point trust-region technique, the smoothed nonlinear programming problem is transformed into an equality-constrained optimization problem with bounded variables by using an active-set approach and the penalty method. To solve the equality-constrained optimization problem with bounded variables, Newton's interior-point method is used. But it is a local method so it might not converge if the starting point is far from a stationary point. To overcome this problem, the nonmonotone trust-region strategy is used. It is generally promising and efficient, and it can overcome the aforementioned shortcomings.

Applications to mathematical programs with equilibrium constraints are given to clarify the effectiveness of the proposed approach. Numerical results reflect the good behavior of the proposed technique and the computed results converge to the optimal solutions. It is clear from the comparison between the solutions obtained by using Algorithm (3.8) and the algorithm of ^[18], that Algorithm (3.8) is able to find the optimal solution of some problems with fewer iterations, fewer function evaluations, and less time.

Furthermore, the usefulness of the CHKS smoothing function with the nonmonton active interior-point trust-region algorithm in efforts to solve NBLP problems was examined by using a real-world case about a watershed trading decision-making problem. Numerical experiments show that the suggested method surpasses rival algorithms in terms of efficacy.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Conflict of interest

The authors declare that they have no competing interests.