A two-grid $P_0^2$-$P_1$ mixed finite element scheme for semilinear elliptic optimal control problems

1.   Introduction

Consider the constrained system

where the function $F:\Psi \longleftrightarrow \mathbb{R}^n$ is a continuous, monotone, and $\Psi\subset\mathbb{R}^n$ is nonempty, closed and convex. The monotonicity of $F$ implies

This problem has applications in financial forecasting problems ^[1] and compressive sensing ^[2]. The Newton method ^[3], quasi-Newton methods ^[4], and conjugate gradient (CG) methods ^[5] are among the most often used approaches for solving it. When dealing with large-scale nonsmooth algebraic systems, CG methods are seen to be the most successful. The iterative technique for the CG algorithm is as follows:

the vector $d_k$ is given by

where $F_k = F(x_k)$, and $\beta_{k}$ is the CG formula that differentiates CG methods. One of the most efficient CG formulas with restarting feature is the Polak-Ribi$\grave{e}$re-Polyak (PRP) ^[6] formula defined by

where $y_k = F_{k+1}-F_{k}$. Despite the good features of the formula (1.5), it doesn't satisfy the sufficient descent condition (SDC). This opens up the possibility for future investigations on the modification of the PRP formula to satisfy the SDC. There are different PRP algorithms for solving the general monotone system of nonlinear equations in the literature. The first PRP algorithm for solving monotone algebraic systems was proposed by Cheng ^[7] where he extended the default PRP using the hyperplane technique discussed in ^[8]. Subsequently, Yu ^[9] proposed the derivative-free PRP for a large-scale monotone system that uses the backtracking line search procedure. Ahookhosh ^[10] proposed a descent three-term PRP algorithm for solving large-scale nonlinear monotone equations that combine extended PRP formula with the projection technique. Yuan and Zhang ^[11] proposed another descent-modified PRP CG algorithm for monotone nonlinear equations. Sabi'u and Gadu ^[12]proposed a projection-based hybrid FR and PRP algorithm for the monotone system. Others are; the improved three-term PRP CG method ^[13], spectral modified Polak–Ribiére–Polyak projection conjugate gradient method for solving monotone systems of nonlinear equations ^[14], the new hybrid algorithm for solving large-scale monotone nonlinear equations ^[15], and the accelerated CG algorithm for solving nonlinear monotone equations and image restoration problems ^[16].

Moreover, for the convex-constrained system (1.1), Jinkui ^[17] proposed a spectral PRP projection algorithm for solving nonlinear monotone equations with convex constraints. Followed by; the projection-based PRP-like algorithm ^[18], the modified spectral PRP conjugate gradient projection method for solving large-scale monotone equations and its application in compressed sensing ^[19], the new hybrid prpfr CG method for solving nonlinear monotone equations and image restoration problems ^[20], the modified PRP CG projection method for solving large scale nonlinear convex constrained monotone equations ^[21], the modified PRP CG method with weaker conditions ^[22], the PRP-like algorithm for monotone operator equations ^[23], and the modified PRP-type conjugate gradient projection algorithm for solving large-scale monotone nonlinear equations with convex constraint ^[24]. In this work, we will propose a scaled PRP CG gradient algorithm for constrained nonlinear systems in which the scaling parameters are determined by approaching the quasi-Newton direction and exploiting the advantages of the well-known Barzilai-Borwein technique. The appealing features of our algorithm are: (i) derivative and matrix-free algorithm, (ii) exploiting best out of the PRP formula numerically and theoretically in contrast to conventional scaled CG algorithms, (iii) the application of the proposed algorithm in solving the two-joint planar robotic manipulator problems, and finally, (iv) the proposed directions will satisfy the SDC independent of the line search procedure.

The next section will discuss the proposed algorithm using the scaling strategy, Section 3 contains the global convergence of the method. The numerical experiment will be presented in Section 4, and the last section will contain the concluding remarks.

2.   A scaled PRP CG method

The scaling strategy is proved to be an efficient strategy for enhancing numerical and theoretical performances of the CG methods, see ^[25,26]. In line with this, we scaled the PRP CG (SPRP) formula to get

such that $\left|\gamma\right|\leq1$. We are going to propose two effective ways of computing the parameter $\gamma$ at each iteration.

2.1.   Based on quasi-Newton method

We are going to propose an optimal choice for the scaling parameter $\gamma$ at every iteration by tending the proposed direction to approach the general quasi-Newton direction. Starting with the quasi-Newton equation

where $B_{k+1}$ is a positive definite and symmetric Jacobian approximate. Now, the general quasi-Newton direction is as follows:

where the matrix $H_{k+1}$ is the inverse of $B_{k+1}$. Now, using (1.4) and (2.1), the SPRPR search direction can be written as

In order to incorporate Jacobian approximation information in the quasi-Newton direction into our scheme, since the quasi-Newton schemes employ the actual Jacobian approximations, from (2.3), (2.4), and also by the virtue of equality relation, we get

Multiplying (2.5) by $-B_{k+1}s_k^T$ to obtain

Solving for $\gamma$ after eliminating the matrix $B_{k+1}$ in (2.6) by using (2.2) to get

We further achieved $\left|\gamma_k^1\right|\leq1$ by using the Polak-Ribi$\grave{e}$re-Polyak non-negative restriction to obtain

2.2.   Based on Barzilai-Borwein approach

We will again propose another optimal choice for the scaling parameter $\gamma$ by using some prominent features of the Barzilai-Borwein approach. The Barzilai-Borwein ^[27] proposed some prevalent choices for the scaling parameters given by

By considering Perry's point of view, the scaled SPRP search direction can be written as

where $\hat{H_{k+1}}$ is defined by

Now, we aim at taking the advantage of the Barzilai-Borwein ^[27] approach to propose the parameter $\gamma$ by minimizing

where $\left\|.\right\|_F$ stands for the Frobenius matrix norm. Moreover, if we assigned $Q = \hat{H_{k+1}}-\omega_kI$ and utilized the relation $\left\|Q\right\|^2_F = {\text{Trace}}(Q^TQ)$, we arrived at

In addition, we ensure that all of the prevalent choices for the Barzilai-Borwein scaling parameters are properly integrated into our scheme by suggesting that

with $0 < \omega_{\min} < \omega_{\max} < \infty$. We proposed the following modified version of (2.13) that satisfies $\left|\gamma_k^2\right|\leq1$ as

Furthermore, to guarantee that the sufficient descent condition is fulfilled independent of the line search process, we provided the following modified SPRP directions

and,

where

Below we present the proposed algorithm.

Algorithm 2.1. The scaled PRP projection-based CG algorithm (SPRPCG)

Step 0. Start by initializing: $\epsilon\geq0$, $b\in(0, 1)$, $\theta > 0$, $\tau > 0$, $a > 0$, $0 < \omega_{\min}\leq\omega_{\max}$ and $x_0\in \bf{R^n}$. Set $k = 0$ and $d_0 = -F_0$.

Step 1. If $\left\|F_k\right\|\leq\epsilon$, terminate, otherwise continue with Step 2.

Step 2. Determine the scaled PRP CG direction $d_k$ via (2.16) or (2.17), where $s_k = u_k-x_k$, and $y_k = F(u_k)-F_k+bs_k$.

Step 3. Set $u_k = x_k+\alpha_kd_k$ and compute $\alpha_k = \max \left\{ \tau\theta^i: i = 0, 1, 2, \cdots\right\}$ satisfying

Step 4. If $u_k\in\Psi$ and $\left\|F(u_k)\right\| = 0$ stop, otherwise

where $A$ is the projection operator, and

Step 5. Set $k = k + 1$ and repeat from Step 1.

3.   Convergence analysis

This section describes the global convergence of the proposed algorithm under the Lipschitz continuous and monotonicity assumptions on $F$.

Lemma 3.1. For all $k\geq0$, the line search (2.19) is well-defined.

Proof. Suppose there exists $k_0\geq0$ such that for $i = 0, 1, 2, \cdots$,

Let $i\rightarrow \infty$ in (3.1), we have

Since the proposed directions satisfiy the sufficient condition, we have

Thus, the relations (3.2) and (3.3) yield a contradiction. Hence, the proof is complete.

Lemma 3.2. ^[25] If the sequences $\{x_{k}\}$ and $\{v_{k}\}$ are produced by the SPRPCG algorithm, then for some $M > 0$ we have

Proof. Now, from the line search (2.19), we have

Assume that $x^{*}\in \Omega$ such that $F(x^{*}) = 0$, using the monotonicity of $F$, we have

Using (3.5) and (3.6) to get

By the definition of $v_k$, and the Cauchy Schwarz inequality in (3.7), we obtain

thus, we have

Therefore, $\{\|x_k-x^{*}\|\}$ is a decreasing sequence and hence convergent. Now, utilizing (3.9) and the Lispchitz continuity of $F$, we get

By (2.19), monotonicity of $F$, and the Cauchy Schwarz inequality, we obtain

which gives

From (3.12), we see that the sequence $\{u_k\}$ is bounded. Also, from (3.8), we obtain

This implies

The proof is now complete.

Theorem 3.1. The SPRPCG algorithm converges, i.e.,

Proof. Assume that (3.15) is not true, i.e., there exists $q > 0$ such that

Again from the Lipschitz continuinity of $F$ it is easy to see that $\left\|y_k\right\|\leq (L+a)\left\|s_k\right\| = N$, where $L$ is assumed to be the Lipschitz constant. Thus,

Now, let $h = 2\frac{MN}{\left\|q^2\right\|}$, therefore from (3.17) we get

Hence, for $Z = \max\{\left\|d_{1}\right\|, \left\|d_{2}\right\|, \cdots, \left\|d_{k_0}\right\|, \frac{M}{1-h}+\left\|d_{k_0}\right\|\}$, we obtain the boundedness of the proposed direction. Similarly, we can show that (2.17) is also bounded. From the sufficient condition and the Cauchy-Schwarz inequality we get

From (3.19) we obtain

Therefore, from (3.14) and (3.20) we have

Suppose $\alpha_k^{\prime}$ doesn't safies the line search procedure (2.19), i.e.,

It is obvious that $\left\lbrace x_k\right\rbrace$ is bounded and has an accumulation point $\hat{x}$ and an infinite index set $D_1$ such that $\lim_{k\in D_1}x_k = \hat{x}$. It also follows that $\left\lbrace d_k\right\rbrace_{k\in D_1}$ is bounded. Therefore, there exist an infinite index set $D_2 \subset D_1$ with an accumulation point $\hat{d}$ such that $\lim_{k\in D_2}d_k = \hat{d}$. Now, taking limit in (3.22), we obtain

Also taking the limit in (3.19), we get

The last two inequalities give a contradiction and the proof is complete.

4.   Numerical experiment

The proposed scaled PRP algorithm is numerically compared to the MPRPA ^[24], DFSP ^[26], and STCGM ^[25] algorithms, using published values in the compared papers. Furthermore, we initialized the following variables for the SPRPCG algorithm: $b = 0.2$, $a = 0.0001$, $\theta = 0.99$, $\tau = 1$, $\omega_{\min} = 0.0001$, $\omega_{\max} = 10^4$, $\Psi = R^n_+$, and $\epsilon = 10^{-10}$. Matlab14 is used to run all algorithms on an Intel Core i5 8th Generation personal computer. We carried out our experiment using the following test problems:

Problem 4.1. Reference [25] $F(x_i) = 2x_i-\sin\left|x_i\right|$, for $i = 1, 2, 3, \ldots, n$, and ${\Omega = {\bf {R}}_+^n}$.

Problem 4.2. Reference ^[25] $F_i(x) = 4x_i+(x_{i+1}-2x_i)-\frac{x_{i+1}^2}{3},$ for $i = 1, 2, \ldots, n-1$, $F_n(x) = 4x_n+(x_{n-1}-2x_n)-\frac{x_{n-1}^2}{3}$, and ${\Omega = {\bf {R}}_+^n}$.

Problem 4.3. Reference ^[26] $F_i(x) = e^{x_i}-1, \quad i = 1, 2, \dots, n$, and ${\Omega = {\bf {R}}_+^n}$.

Problem 4.4. Reference [25] $F_1(x) = \cos (x_1)-9+3x_1+8\exp(x_2)$ $F_i(x) = \cos (x_i)-9+3x_i+8\exp(x_{i-2}),$ for $i = 2, 3, \dots, n,$ and ${\Omega = {\bf {R}}_+^n}$

Problem 4.5. Reference [25] $F_1(x) = e^{x_1}-1,$ $F_i(x) = e^{x_i}+x_{i-1}-1, \quad i = 2, 3, \cdots, n-1$, and ${\Omega = {\bf {R}}_+^n}$.

From Tables 1–5, we used eight initial points, namely; $x_0^1 = (2, 2, \ldots, 2)$, $x_0^2 = (1, \frac{1}{2}, \frac{1}{3}, \ldots, \frac{1}{n})$, $x_0^3 = (1, 1, \ldots, 1)$, $x_0^4 = (\frac{1}{n}, \frac{2}{n}, \ldots, 1)$, $x_0^5 = (n-\frac{1}{n}, n-\frac{2}{n}, \ldots, n-1)$, $x_0^6 = (2, \frac{2}{2}, \frac{2}{3}, \ldots, \frac{2}{n})$, $x_0^7 = (1-1, 1-\frac{1}{2}, 1-\frac{1}{3}, \ldots, 1-\frac{1}{n})$ and $x_0^8 = (-3, -3, \ldots, -3)$. In which the terms ITER, FVAL, TIME, and NORM stand for the number of iterations, function evaluation, computing time and the norm of $F_k$ at the stopping criteria respectively.

When compared to the MPRP, STCGM, and DFSP algorithms, the proposed approach solved the Problem 1 with the lowest number of ITER and FVAL. However, for the computing time, the MPRP algorithm won initially and failed to STCGM as the dimension increases. Likewise for Problems 2–5. Moreover, it is remarkable to note that the MPRP algorithm is competing with the STCGM algorithm for computing time with respect to all the considered test problems. Finally, the proposed algorithm outperforms the MPRPA, STCGM, and DFSP algorithms for ITER, and FVAL, but competes with the STCGM algorithm with respect to TIME.

Furthermore, we used the well-known Dolan and Mor$\acute{e}$ ^[28] performance profile to produce a more exact numerical comparison using the bulk data in all the tables. Figure 1 three subfigures indicate that the SPRPCG algorithm with two directions reflects the upper left curves with regards to 1a and 1b. However, for the 1c, the SPRPCG initially won in the fractions of 0 to approximately 0.9, and later failed to STCGM algorithm. As a result, when compared to the MPRPA, STCGM, and DFSP algorithms, the SPRPCG algorithm has fewer iterations, and function evaluations on average and competes with the STCGM algorithm for CPU time.

4.1.   The motion control of robotic manipulator

IIn this part, we will tackle the motion control of a two-joint planar robotic manipulator using the SPRP projection-based algorithm. The following discrete-time kinematics equation describes the position level of the two-joint planar robot manipulator

where the kinematics function $f(.)$ is defined by

in which, $m_j$ for $j = 1, 2$ is the length of the $j-th$ rod, $p_1 = \cos(\eta_1)$, $r_1 = \sin(\eta_1)$, $p_2 = \cos(\eta_1 + \eta_2)$, and $r_2 = \sin(\eta_1 + \eta_2)$. Also, the vector $\sigma_k\in\mathbb{R}^2$ stands for the end effector position vector, and the joint angle vector is represented by $\eta_k\in\mathbb{R}^2$, for more detail about the kinematics equation we refer readers to ^[29] and the references therein. In our experiment we aimed at solving the series of optimization problems at a time interval $t_k\in\left[0, t_f\right]$ defined as follows:

In this experiment, $m_j = 1$ is assigned for $j = 1, 2$, and the end effector is regulated to trace a Lissajous curve as provided by

For SPRP algorithm, we initialized $b = 0.2$, $a = 0.08$, $\theta = 0.2$, $\tau = 1$, $\omega_{\min} = 0.0001$, $\omega_{\max} = 10^4$, $\eta_0 = \left[0, \frac{\pi}{3}\right]$, and the end task duration $t_f = 10s$. We further divided the task duration $[0, 10]$ into 200 equal parts.

Figure 2 depicts the experimental findings generated by the SPRP algorithm, including the robot trajectories synthesized by the SPRP algorithm, the end plots of the effector trajectory and the intended route, and the error of the SPRP algorithm on both the vertical $x$ and $y$ axes. It is not difficult to see from Figures 2a and 2b that the SPRP algorithm completes the specified task satisfactorily. Furthermore, Figures 2c and 2d show that the resultant error is around $10^{-5}$.

5.   Conclusions

A robust descent-scaled PRP projection-based technique for solving monotone nonlinear problems with convex constraints was given. The monotone and Lipschitz continuous assumptions are applied to accomplish the proposed algorithm's global convergence. A detailed numerical comparison with the related algorithms indicated that the proposed technique is much more efficient than the existing algorithms. Furthermore, the proposed technique is used to solve the motion control problem of a two-joint planar robotic manipulator with extremely little error. Finally, the provided technique may be used for various problems, including solving unconstrained optimization problems, generic algebraic nonlinear systems, and so on.

Acknowledgments

This research received funding support from the NSRF via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation, (grant number B05F650018).

Conflict of interest

The authors declare no conflict of interest.