Modified Gradient Flow Method for Solving One-Dimensional Optimal Control Problem Governed by Linear Equality Constraint

This study presents a computational technique developed for solving linearly constraint optimal control problems using the Gradient Flow Method. This proposed method, called the Modified Gradient Flow Method (MGFM), is based on the continuous gradient flow reformulation of constrained optimization problem with three-level implicit time discretization scheme. The three-level splitting parameters for the discretization of the gradient flow equations are such that the sum of the parameters equal to one (θ1 + θ2 + θ3 = 1). The Linear and quadratic convergence of the scheme were analyzed and were shown to have first order scheme when each parameter exist in the domain [0, 1] and second order when the third parameter equal to one. Numerical experiments were carried out and the results showed that the approach is very effective for handling this class of constrained optimal control problems. It also compared favorably with the analytical solutions and performed better than the existing schemes in terms of convergence and accuracy. DOI:10.46481/jnsps.2022.589


Introduction
Many authors, in literature, have attempted the class of linearly constrained optimal control problems using the classical analytical techniques of the Hamilton-Jacobi-Bellman (HJB) or "Pontryagin Maximum Principle" by Pontryagin et al. in [14]; though this approach may proof difficult in some cases when inclusive of non-differentiable delay terms. The numerical approach has also been fully handled by several authors like Betts * Corresponding author tel. no: +2348032541380 Email address: dawodukazeem@futa.edu.ng (Kazeem Adebowale Dawodu) in [4], Olotu and Dawodu in [11,12], Adekunle [1] and Akeremale [2], using the direct numerical methods. This method deploys the "first discretized -then -optimize" technique which is known for its amenability to well-structured discretization schemes and optimization solvers after the unconstrained formulation of the control problem using any of the penalty function methods or method of multipliers. The direct method of nonlinear programming formulation tends to eradicate the problem of ill-conditioning associated with control problems with delay dynamical systems. However, this research paper seeks for a better approach in the formulation of the control operators in the optimization of the discretized optimal control problem. The gradient flow method was first applied to nonlinear programming problems by Evtushenko [5] and later improved by Evtushenko and Zhadan [6]; while the application to unconstrained problems was by Behrman [3] and the unified approach to nonlinear optimization problem was by Wang et al. [15]. The unified gradient flow algorithm exhibits linear convergence for θ ∈ [0, 1) and quadratic convergence at θ = 1. As an improvement over [15], a modified gradient flow method MGFM was proposed using the gradient formula for Lagrange multiplier formulation. The development of this method requires a continuous gradient flow approach based on a three level implicit time discretization scheme with splitting parameter θ 1 , θ 2 and θ 3 . The convergence of the solution of the modified discretized gradient flow equation to a local minimum of the original problem, either linearly or quadratically depending on the choice of θ, will as well be demonstrated. However, the technique can also be adapted to handle inequality constraint problems by introducing slack variables. Three test problems were solved and their numerical results presented in comparison with existing methods. This is to demonstrate the effectiveness of the method in solving constrained nonlinear programming problems with either continuous variables or a mixture of continuous and discrete variables.
The structure of this paper is as follows: Section 3 is dedicated to the idea of gradient flow approach for unconstrained optimization while Section 4 presents the modified gradient flow algorithm for constrained optimization in Section 2 using the 3level splitting parameters and as well as the convergence analysis. The Algorithm was experimented on two numerical examples and comparison made with the Conjugate Gradient Method (CGM) in Section 4.

Statement of problem
The general form of the One-Dimensional Optimal Control Problem is given as where Thus, x(t) is the state variable which describes the state of the system at any point in time and u(t) is the control variable in the optimization problem. However, for linear-quadratic problem, the objective functional is written as F(t, x(t), u(t)) = (px 2 (t) + qu 2 (t)) while the linear differential constraint is written as

Research Methodology
This paper considered both the theoretical and numerical analyses of the problem in the development of the algorithm for solving eqns. (1) to (3). The discrete form of the quadratic objective functional and the linear differential constraint are constructed using the Trapezoidal and Euler method respectively. The recurrence relation for both the objective functional and constraint were developed to help construct their respective matrix operators. Furthermore, the symmetry and positive definite properties of the matrix operators were analyzed to guarantee their consistency and invertibilty in the developed MGFM.

Preamble
Let f : X → R be a functional such that the solution x(t) is a minimizer of f . The gradient flow method, with an initial point x(0) ∈ X, seeks to find a minimizer of f by the curve x(t) defined by the differential equation where ∇ f is called the gradient of f and the solution is the called the integral curve such that at each instance proceeds in the direction of the steepest descent of f . In an unconstrained optimization problem, the equilibrium points of the gradient flow are the zeros of ∇ f , which are exactly the minimizers of f . The gradient flow solves the problem min f (x), x(0) = x 0 such that in every trajectory x(t) of the gradient flow, then f (x(t)) → p * for p * is the minimizer of f . For clarity, the derivative x (t) is descretized using the forward Euler step which yields or the backward discretization scheme which yields The backward Euler method can be applied for the numerical integration of gradient flow differential equation, also known as the Proximal minimization method expressed thus: However, the proximal gradient flow algorithm can be interpreted as gradient flows, where g is differentiable. Using the forward-backward integration of the gradient flow, eqn. (7) can be expressed as where the forward-backward splitting method above numerically integrates the gradient flow differential equations. The time-stepping discretization of eqns. (5) and (6) is expressed as follows: where θ ∈ [0, 1] is a parameter. When θ = 0, the above discretization is the explicit Euler's scheme and on the other hand, the implicit backward Euler when θ = 1.

Discretization of the Constraint function
The differential constraintẋ(t) = f (t, x(t), u(t)) = ax(t) + bu(t) was discretized by slotting it into the Euler's scheme x i+1 = x i +ẋh with t i = t 0 + ih, to arrive at where c = 1 + ha, d = hb and i = 0, 1, 2, ..., N − 1. And upon expansion and re-arrangement of eqn. (13) into matrix form yields GZ = k express as where Z is of dimension (2N + 1) × 1 and G is of dimension N × (2N + 1) as shown below; with the entries defined as follows: g i j = 1 for 1 ≤ i ≤ N and j = i; g i j = −d for 1 ≤ i ≤ N and j = N + i; g i j = −c for 2 ≤ i ≤ N and j = i − 1 while g i j = 0 otherwise. Hence, eqns. (11) and (14) give the quadratic representation for the constrained minimization problem as

Lagrangian formulation of the Constraint problem
The unconstrained form of the optimal control problem can be obtained using the Lagrangian function on eqn. (15) as stated below: where λ = (λ 1 , λ 2 , ..., λ N ) T ∈ R N is the Lagrangian multiplier. The pair of points (Z * , λ * ) is said to be a stationary pair points of eqn. (17) if the following first order necessary optimality conditions hold: where denote the gradient vectors of L and J with the dimensions (2N + 1) × 1, N × 1 and (2N + 1) × 1 respectively. Similarly, we have the Jacobian of g with dimension N × (2N + 1) represented below: In order to fulfill the above optimality conditions, the continuous gradient flow reformulation with the given initial condition Z 0 is derived as follows; Substituting eqn. (18) into (21) yields, such that the choice of λ 0 (Z) is to satisfy in the least square sense, where τ > 0 is a parameter. Substituting eqn. (18) into (24) yields, Obtaining the least squares solution of λ 0 (Z) in en. (25) yields, where the superscript ' †' refers to the pseudo inverse of g T Z (Z) in Layton [9].
Theorem 3.1. The constraint equations for the discretized optimal control problem g(Z) satisfies the linearly independent constraint qualification (LICQ) for any feasible point Z * . The LICQ states that the gradient of the active inequality constraints and the gradient of the equality constraints are linearly independent at Z * . eyiu LICQ holds if g 1Z (Z), g 2Z (Z), ..., g NZ (Z) are linearly independent for any feasible point Z * . Thus implies that α i = 0. Hence, the dimension of the row space of G must be N.
Applying the row-reduction process to G in eqn. (3.3), we havē G * = [I|V] written as where I ∈ R N×N and V ∈ R N×(N+1) . The entries g * i j ofḠ * can be defined as follows: Thus, rank ofḠ * equals rank of G and α i = 0 are the parameters that make g 1Z (Z), g 2Z (Z), ..., g NZ (Z) linearly independent for any feasible point Z * and G(Z) satisfies LICQ. From theorem (3.1), we have for any feasible point Z * thatḠ * ∈ R N×(2N+1) is re-structured into a N × N non-singular (invertible) Gram matrix G * given as Multiplying both sides of eqn. (25) by g Z (Z) arrives at and simplified to get Pre-multiplying eqn. (33) by G −1 * yields and solving for λ 0 (Z) to obtain where eqns. (27) and (35) . However, eqn. (35) is true when Z is sufficiently close to an equilibrium point Z * . The choice of λ 0 (Z) has the advantage of reducing the complexity in the evaluation of the gradient of λ 0 with respect to Z. Thus substituting λ 0 (Z) in eqn. (35) above into eqn. (23) yields; In testing the convergence of the gradient flow formulation expressed in eqns. (23) and (24), we then considered the LICQ theorem below.

Theorem 3.4 (Sylvester's Criterion).
A real, symmetric matrix is positive definite if and only if all the principal minors are positive definite. Gilbert [7].
Considering the implementation of Theorem 3.4 on the matrix Q, we seek to establish that the matrix is positive definite. Firstly, Q is real since for every q i, j ∈ Q, q i, j ∈ R (i.e. all the diagonal entries in eqn. (67) are real).
However, we also seek to establish that the principal minors Q i , ∀i = 1, 2, 3, ...N, of Q are all positive definite. When i = 1, 2, 3, and 4, Q 1 , Q 2 , Q 3 and Q 4 are 1, 1 + c 2 , 1 + c 2 + c 4 and 1 respectively; and are all positive ∀c, d > 0. We then assumed, by mathematical induction, that it is true for Q k when i = k and Q k+1 when i = k + 1. And by Cholesky, it was shown that Q k+1 is positive definite and that all the principal minors of Q are positive definite [See Appendix]. We, therefore, concluded that the matrix Q = GG T is positive definite.

Construction of the Modified Gradient Flow Method (MGFM)
We now considered the discretization of the formulation of eqns. (23) and (24) using three level implicit time discretization scheme with splitting parameters θ 1 and θ 2 where θ 3 = (1 − θ 1 − θ 2 ). Let 0 = t 0 < t 1 < t 2 < t 3 < ... < t k = t be a sequence time points for the time t ≥ t 0 , h k = t k+1 − t k and δ Z k+1 = Z k+2 − Z k being the sequence of solution distance between two successive 151 equilibrium points. We discretize eqn. (23) by the following implicit time stepping algorithm; Simplified to Employing Taylor's expansion about (Z k+2 , λ 0 k+2 ) arrives at where, Neglecting the higher order terms in eqn. (73) and solving the resultant equation for Z k+2 yields For equal step size, h k+1 = h k with θ 1 = 0 we obtain Z k+2 expressed below In collecting like terms and factorizing eqn. (75), it becomes Making Z k+2 the subject of formula in eqn. (76) yields for any θ 2 , θ 3 ∈ [0, 1] and θ 2 + θ 3 = 1 by the time-stepping algorithm for convex functions. Also, from eqn. (76), we obtained Given that δ k+1 = Z k+2 − Z k and θ 2 = 0, eqn. (77) becomes Thus, for any initial guess Z 0 , eqn. (77) defines a series converging to the solution of the optimal control problem stated in eqns.
(1) to (3). Hence, the Modified Gradient Flow (MGF) Algorithm for the Discretized Optimal Control Problem (OCP) is stated as follows:

Error and Convergence Analysis
The convergence analyses of the MGFM requires us to shown the rate and type of convergence of the algorithm for varying values of parameters θ1, θ2 and θ3. We shall show that the algorithm is super-linearly convergent for some θ 2 and θ 3 and quadratically convergent if θ 1 = θ 2 = 0 and θ 3 = 1.  1) and (2). If the initial guess Z 0 is sufficiently close to Z, then we have 1. Z k converges linearly to Z˚a if θ 1 , θ 2 , θ 3 ∈ [0, 1] and h k > 0 is sufficiently small; 2. Z k converges to Z˚a quadratically if θ 2 = 0, θ 3 = 1 and h k → ∞.
We shall prove the two cases separately.
• Case 1. θ 1 , θ 2 , θ 3 ∈ [0, 1]. From eqn. (77) we have: Hence, Let (Z * , λ * ) be the minimum points of eqns. (16) and (17) and e k+1 = Z k+1 − Z * . Subtracting Z * from both sides of the above equation and noting that we have Thus, by the Mean Value Theorem, where ζ k denotes a point in the line segment connecting Z k and Z * . Taking the norm of both sides, we have where From eqn. (81), we see that if γ(Z k , λ 0 k , ζ k , θ 1 , θ 2 , θ 3 , h k ) < 1, then e k converges to zero linearly. To show this, we need to prove that H Z k is positive definite. But by theorems and , H Z k is positive definite. With these, when Z k is sufficiently close to Z * , from eqn. (82), it is seen that if θ 1 , θ 2 , θ 3 ∈ [0, 1] and h k sufficiently small, we have where λ k min and λ k max denote the minimum and maximum eigenvalues of H Z k , respectively. Therefore, eqn. (81) implies that lim k→∞ e k = 0 is linear. This proves that Z k converges linearly to Z˚a if θ 1 , θ 2 , θ 3 ∈ [0, 1] and h k > 0 is sufficiently small.
From eqn. (77), when θ 3 = 1 and h k = h, we have When h → ∞, the above equation becomes Eqn. (84) coincides with the equation resulting from the application of the Newton's method to L Z = 0. Hence, Z k converges quadratically to Z * if Z 0 is sufficiently close to Z * .

Examples and Numerical Results
The numerical examples and respective results are presented below based on the construction of modified gradient flow method (MGFM) for discretized optimal control problems. The results show the robustness and efficiency of the scheme (MGFM) as compared with existing schemes.

Example 1
Considering the optimal control problem presented by [1] Minimize The analytical solution is given as and the control variable is given as u(t) = − 1.4770e −5.3852t . The solution of eqns. (85) and (86) were obtained, using h = 0.1 and h k = 20, by the modified gradient algorithm for discretized optimal control problems and the results were compared with the existing algorithms as shown in table 1 below. However, the step-size h = 0.1 was chosen for the purpose of illustration. Meanwhile smaller choices of step-size could be used for better refinement.
Algorithm objective functional value Analytical Method 0.2953894 Adekunle [1] 0.3033885 Akeremale [2] 0.2955932 MGFM 0.2955879   Table 2 above shows the errors in the state and control variables for example 1 while the errors generated with the MGFM were minimal compared with the errors generated from Adekunle [1] and Akeremale [2]. Figures 1 and 2 below clearly show the Algorithm Objective Functional Value Analytical Method 0.5647 Adekunle [1] 0.5746 Akeremale [2] 0.5649 MGFM 0.5648 comparisons of the behaviors of the state and control trajectories to the various algorithms respectively. It was observed that the trajectories for both the state and control variables overlap due to closeness in results.

Example 2
Considering the optimal control problem presented by [1] MinimizeJ(x, u) = 1 0 (x 2 (t) + u 2 (t))dt,    respectively. It was observed that the trajectories for both the state and control variables overlap due to closeness in results.

Conclusion
This paper developed a new scheme for the solution of a quadratic optimal control problem with linear equality constraint. The method is a gradient flow formulation that requires dis-cretization techniques and a splitting parameter approach for the gradient of the equation whose solution converges to a local minimum of the original problem; either linearly or quadratically. Two linear equality constrained optimal control problems were considered and their results compared with the analytical and two known existing schemes. Thus, it was observed that the new scheme compared favorably with the analytic solutions and performed better than the Conjugate Gradient Method used in [2] in terms of convergence and accuracy. It is therefore obvious that this developed algorithm, due to its robustness and efficiency, would provide a new platform for solving optimal control problems with or without box constraints, equality or inequality constraints or with delay or non-delay dynamical systems.
By Cholesky, Q k+1 is said to be positive definite if there exists a lower triangular matrix L i, j such that Q k+1 = LL T . That is, where H k is a lower triangular matrix with positive diagonal entries. However, eqn. (5) such that Q k = H k H T k , β T k+1 = L T k+1,1 H T k , β k+1 = H k L k+1,1 and σ 2 = L T k+1,1 L k+1,1 + L 2 k+1 L k+1 . Thus, since the diagonal entries of H k are greater than zero, it is non-singular and therefore the linear system of the equation has a unique solution given by L k+1,1 = H −1 k β k+1 and a positive value for L k+1,k+1 can be obtained if σ 2 − L T k+1 L k+1 > 0. Hence, Substituting the values of β k+1 and Q k , we have where (H k H T k ) −1 = Q −1 k , Q k = [σ 2 − L T k+1,1 L k+1,1 ].
Thus L i > 0 arises from the fact that G is positive definite (LICQ). Hence, the pair (Z(t), λ 0 (Z(t))) tends to (Z * , λ * ) as t → ∞ provided the starting point Z 0 is sufficiently close to point Z * .