Direct transcription methods based on fractional integral approximation formulas for solving nonlinear fractional optimal control problems

This paper presents three direct methods based on Gr\"{u}nwald-Letnikov, trapezoidal and Simpson fractional integral formulas to solve fractional optimal control problems (FOCPs). At first, the fractional integral form of FOCP is considered, then the fractional integral is approximated by Gr\"{u}nwald-Letnikov, trapezoidal and Simpson formulas in a matrix approach. Thereafter, the performance index is approximated either by trapezoidal or Simpson quadrature. As a result, FOCP are reduced to nonlinear programming problems, which can be solved by many well-developed algorithms. To improve the efficiency of the presented method, the gradient of the objective function and the Jacobian of constraints are prepared in closed forms. It is pointed out that the implementation of the methods is simple and, due to the fact that there is no need to derive necessary conditions, the methods can be simply and quickly used to solve a wide class of FOCPs. The efficiency and reliability of the presented methods are assessed by ample numerical tests involving a free final time with path constraint FOCP, a bang-bang FOCP and an optimal control of a fractional-order HIV-immune system.

M A N U S C R I P T

Introduction
Fractional calculus may be considered an old and yet an interesting topic [1]. It deals with the investigation of integrals and derivatives of an arbitrary order. At the initial stage, the fractional derivative was a mathematical tool without a tangible application. But presently, fractional calculus has wide applications in various disciplines like viscoelastic materials [2,3], bioengineering applications [4,5], signal processing [6,7], mathematical finance [8,9], control theory [10,11] and many other areas [12,13]. The nonlocal nature of the fractional calculus has given it a unique characteristic in modeling some complex systems with memory and hereditary properties that arise from physical, engineering and even economic processes [14,15]. More elaborate details on the theory and A C C E P T E D M A N U S C R I P T applications of fractional calculus can be found in [2,16,17,18].
Application of fractional calculus to optimal control problems has given birth to another new specialization, known as fractional optimal control [18,19,20,21]. According to [22], a Fractional Optimal Control Problem (FOCP) is an optimal control problem in which the performance index and/or the dynamic equations contain at least one fractional derivative operator. In simple terms, FOCPs are a generalization of the integer-order optimal control problems, which are obtained by replacing integer-order derivatives with fractional ones. Recently, some interesting and real-life models of FOCPs have been presented by the researchers in [23,24,25,26,27].
Like integer-order optimal control problems, the numerical methods for FOCPs can be categorized into "direct" and "indirect" methods. With indirect methods, the solution is obtained by solving a fractional Hamiltonian boundary-value problem, which is derived from the optimality conditions. Accordingly, the first step in indirect methods is the derivation of first-order optimality conditions. On the other hand, direct methods do not rely on optimality conditions. In these approaches, FOCPs are solved by transcribing them into Nonlinear Programming problems (NLP). Thereafter, a NLP-solver is used to solve the resulting problem [28]. The indirect methods, in comparison with the direct, have some disadvantages, including, (i) difficulties in deriving the Hamiltonian boundary-value problem, especially for problems with path constraints, and (ii) sensitivity to initial guess for state functions as well as costate or adjoint functions. Due to the aforementioned reasons, direct methods are easier to employ than the indirect methods to solve both integer and fractional order optimal control problems.
Direct transcription methods based on local methods, such as trapezoidal and Simpson methods, are widely used to solve integer-order optimal control problems and some softwares, such as SOCS [67] and ICLOCS [68], are developed based on these methods. However, it is somewhat surprising that these methods have not been yet fully explored to solve FOCPs. In this paper, we extend these methods to fractional order optimal control problems. For this purpose, the matrix form of Grünwald-Letnikov, trapezoidal and Simpson approximation A C C E P T E D M A N U S C R I P T formulas for fractional integrals, which are called fractional integral matrices, are derived. These fractional integration matrices are used to reduce the FOCP to a NLP. Thereafter, the resulted NLP is solved by using a well-developed solver to obtain an approximate solution of the FOCP. Moreover, in order to increase the speed of the proposed direct methods, the exact gradient of the objective function and the Jacobian of constraints are derived and supplied to the NLP-solver. Finally, the reliability, efficiency and accuracy of the proposed direct method are demonstrated with four test problems, including a nonlinear and complex FOCP, a FOCP with path and terminal constraints, a bang-bang FOCP and an applied optimal control problem of a fractional order HIV-immune system with memory.
The paper is organized as follows. In Section 2, the notions of fractional derivative and integral are reviewed and three approximation methods for the fractional integral, in matrix form, are introduced. The considered formulation of fractional optimal control problems is stated in Section 3. The detailed implementation of direct methods is presented in Section 4. In Section 5, four numerical examples are provided to show the efficiency and reliability of the proposed methods. Finally, a conclusion is given in Section 6.

Definition and approximation of fractional integral and derivative
In this section, some definitions, together with three approximation formulas for fractional integrals, are reviewed. Moreover, for easy implementation, the matrix forms of these approximation formulas are derived. Definition 2.1 (See, e.g., [16,17]). The left sided fractional integral with order α > 0 of a given function y(x), x ∈ (a, b), is defined as where Γ(·) is Euler's gamma function.
Definition 2.2 (See, e.g., [16,17]). The left sided Caputo fractional derivative of order α ∈ (0, 1) of a func- The first step in developing our numerical methods is to approximate the fractional integral of a given function. For this purpose, we briefly review Grünwald-Letnikov (GL), trapezoidal (TR) and Simpson (SI) formulas for approximating the fractional integral.

The Grünwald-Letnikov formula for approximating a fractional integral
Let the fractional integral of function y be defined on the interval [0, 1]. To approximate the fractional integral of function y, the interval [0, 1] is divided into n equal parts by the following mesh points: x i := ih, i = 0, 1, . . . , n, where h = 1/n. The Grünwald-Letnikov (GL) formula for approximating the fractional integral of function y at x = x i , can be expressed as where G (α) ij are the coefficients of the GL approximation formula, defined by such that h α (see [69,70]). Let Then, all the n + 1 formulas in (3) can be written simultaneously in the following matrix form: where W (α) The matrix W (α) GL is called the fractional GL integration matrix.

The trapezoidal formula for approximating a fractional integral
Let h = 1/n and x i := i h, i = 0, . . . , n. If the given function y is approximated by its piecewise linear interpolant based on the nodes x i , i = 0, . . . , n, and the fractional integral is applied to this approximation, then the trapezoidal (TR) formula is derived as where T (α) A C C E P T E D M A N U S C R I P T (see [69]). If we setā then, we can express the coefficients of the trapezoidal formula (8) as Considering (5), the trapezoidal formula (8) can be expressed in the following matrix form: where W (α) TR is the fractional TR integration matrix defined by

The Simpson formula for approximating a fractional integral
Let n be even and the interval [0, 1] be divided into n parts by the mesh points x i := ih, i = 0, . . . , n. Then the given function y on the intervals [x 2i , x 2i+2 ], i = 0, . . . , n 2 − 2, can be approximated by its corresponding quadratic interpolation polynomial, i.e., y is approximated by a piecewise quadratic polynomial on the whole interval [0, 1]. By replacing y(t) with this piecewise quadratic polynomial approximation in (1), the Simpson formula for approximating the fractional integral of function y on the nodes x i , i = 0, . . . , n, is derived as where S (α) A C C E P T E D M A N U S C R I P T such that the parameters γ Considering (5), the Simpson formula (13) can be expressed in the following matrix form where The Simpson and trapezoidal approximation formulas have been presented in [69]. Due to the fact that these matrix forms have some advantages in the implementation of our method, in this paper we derived the explicit forms of these matrices.

The fractional optimal control problem
In this work, we consider a general formulation for FOCPs, which is described as follows: find the optimal ∈ R p and possibly the terminal time t f that minimize the performance index subject to the fractional dynamic system with the initial and boundary conditions and with the path constraint Here, functions h, g, f , ψ and φ are sufficiently continuously differentiable and defined by the following mappings: Most FOCPs found in the literature, in various disciplines like engineering, control, signal processing and others, belong to the above general class of FOCPs [23,24,25].
For easy application of our numerical method, the time domain [0, t f ] is mapped to the canonical interval [0, 1], by using the affine transformation t → τ t f . Applying this mapping and by noting that the optimal control problem (17) is converted to the following form: It should to be noted that, after applying the mentioned transformation, the symbols of variables in (18) should be changed to new symbols. However, for the sake of simplicity, we shall retain the symbols already used.

The proposed method
We present a direct numerical approach to solve the FOCPs stated in (18). Let the interval [0, 1] be divided into n equal parts by the following mesh points: where h = 1/n. Moreover, for the sake of simplicity, we consider the notations

Discretization of the performance index
The performance index (18a) is discretized by a quadrature formula, such as the trapezoidal, Simpson or other rules. In this way, we have where, in the case of using the trapezoidal rule, and, in the case of using Simpson rule,

Discretization of the fractional dynamic equation
By applying the fractional integral operator 0 I α τ to both sides of the fractional dynamic equation (18b), and by noting that , the above fractional dynamic equations can be converted into the following fractional integral form: By collocating the above equation at τ = τ i , i = 1, . . . , n, we have Considering the initial condition (18c) and notations (19), and by applying GL/TR/SI approximation formulas (3)/(8)/(13), we can approximate the equation (22) as follows: ij are the coefficients of GL/TR/SI approximation formula defined in (4)/(9)/ (14). Similarly, we can discretize the path constraint (18e) as follows: Moreover, the terminal condition (18d) can be discretized as In summary, the FOCP (18) is transcribed into the following nonlinear programming problem (NLP): It should be remembered that the decision variables of the above optimization problem are: and maybe t f . By utilizing an optimization solver, we can solve the optimization problem (26) and find the optimal value of the decision variables, which are the approximations of state and control functions at the points τ i , i = 0, . . . , n.

Reformat of the resulted optimization problem to the classical form
Traditionally, in an NLP, it is preferred that the decision variables are placed in a vector. This helps to analyze and utilize a solver for solving the problem. Here, we reformulate the NLP (26) If we define new functions h and g as then we can express the objective function (26a) as where w is the vector of quadrature weights in (20), i.e., It should be noted that, in case of using trapezoidal quadrature, we have and, for Simpson quadrature, we have The constraints (26b), for i = 0, 1, . . . , n, can be represented in the following matrix form: where X, X 0 and F are p × (n + 1) matrices defined as and W (α) is the fractional GL/TR/SI integration matrix defined in (7)/(12)/ (16). The matrix equation (31) can be reformulated as follows: where vec is the vectorization operator, which converts a matrix into a column vector by stacking the columns of the matrix on the top of each other.
Let I k be the identity matrix of order k and the matrix R be defined as Then, we can express vec (X) based on z as where ⊗ denotes the Kronecker product (or tensor product). Moreover, if 1 n is a row n-vector whose entries are all equal to 1, then we can express vec(X 0 ) based on z as To express vec F W (α) T based on z, we use Theorem 13.26 in [71], which states that for any three matrices A, B, and E for which the matrix product ABE is defined, one has vec(ABE) = (E T ⊗ A)vec(B).

A C C E P T E D M A N U S C R I P T
Now, by using the above equation and by noting that F W (α) T = I p F W (α) T , we conclude that In view of (32), we can write Using the above notation and equations (34)-(36), we can write (33) as To replace the constraints (26c) and (26d) by other constraints based on z, we define functions φ and ψ as Then the constraints (26c) and (26d) can be expressed as The above results are summarized in the following theorem.
Theorem 4.1. The optimization problem (26) can be expressed in the standard form In summary, the solution of the FOCP (17) is reduced to the solution of the NLP (38). In this NLP, w can be w TR or w SI defined in (29) and (30), respectively. Moreover, the fractional integration matrix W (α) can be SI , defined in (7), (12) and (16), respectively. We categorize our methods as follows: GL , then the method is called the "direct GL method".
• If w = w TR and W (α) = W (α) TR , then the method is called the "direct TR method".
• If w = w SI and W (α) = W (α) SI , then the method is called the "direct SI method". For free time problems, the gradient of the objective function can be obtained as For fixed final time problems, the gradient is the same as above, except that the last row is removed.
The Jacobian of constraints (38b), for free final time problems, can be derived as such that In the case of fixed final time problems, the Jacobian of constraints (38b) is obtained as

A C C E P T E D M A N U S C R I P T
The Jacobian of constraint (38c), for fixed final time FOCPs, is derived as while for free final time FOCPs is given by In case of free final time FOCPs, the Jacobian of the inequality constraints (38d) is obtained as where [φ x ] i and [φ u ] i are the r 2 -vectors defined by In case of fixed final time problems, the last column of the Jacobian matrix (44) is removed.
As we can see, the partial derivatives of g, f , ψ and φ with respect to x and u are needed to provide the gradient of the objective function and the Jacobian of the constraints. We can use symbolic computation to supply these partial derivatives.

Numerical examples
We now illustrate the direct methods presented in Section 4 through numerical experiments. We have implemented the direct GL, TR and SI methods using Matlab on a 3.5 GHz Core i7 personal computer with 8 GB of RAM. Moreover, for solving the nonlinear programming problem (38), the solver Ipopt [72] was used, which is based on an interior-point algorithm. In Ipopt, we can adjust the accuracy of solution by the input parameter tolrfun. In our numerical experiments, we set tolrfun=10 −12 . We

Example 1: A FOCP with exact solution
In this example, we consider the following nonlinear FOCP: x(0) = 1, The exact solution for α = 1/2 is where J 0 is the first kind Bessel function of order zero [73]. The domain of this problem is large and the control and state functions have algebraic singularities and an oscillating behavior. Consequently, this example is a suitable test problem for assessing the accuracy and efficiency of the presented methods.
By applying the direct GL and TR methods with n = 100, the obtained control, state and error functions are plotted in Figure 1. We see that, with n = 100, the accuracy of the direct GL method is not satisfactory.
However, the direct TR method generates a solution with reasonable accuracy. The results of the direct GL, TR and SI methods, with various values of n, are reported in Table 1. In this table, to show the impact of using the gradient of the objective function and the Jacobian of the constraints, we compare the CPU times in the two cases "NPD" and "PD". By "NPD" we refer to the CPU time of the methods without providing derivatives  of the objective and constraint functions. In contrast, "PD" refers to the CPU time when the first derivatives of the objective and constraint functions are supplied to the direct methods. In addition, in Table 1, E n (u) and E n (x) are the 2 norm of the error in control and state, i.e.,

A C C E P T E D M A N U S C R I P T
where u i and x i are the obtained control and state functions in τ = τ i by the presented methods with n nodes. Figure 2 visualizes the results in Table 1. In Figure 2 (left), the loglog plot of E n (u) and E n (x) versus n are plotted. Moreover, the regression lines and their slopes are reported too. In Figure 2 (right), the CPU times of the direct GL, TR and SI methods are plotted for the two cases NPD and PD. From Table 1 and Figure 2, we can see that by supplying the gradient of the objective function and the Jacobian of the constraints, the CPU time is significantly reduced. Moreover, we find out that the computational times for direct GL, TR and SI methods are almost the same. However, the direct SI method is more accurate than the direct TR and GL methods. In view of Figure 2

Example 2: A free final time FOCP with path and terminal constraints
In this example, a free final time FOCP with path constraint is considered. Moreover, the state function is forced to lie on a circle at the final time: x(0) = 1, Because of existence of path and terminal constraints, derivation of the optimality condition for this example is difficult. Consequently, solving this problem by indirect methods is cumbersome. We remind that in this paper we use a direct approach, which does not rely on the existence of optimality conditions.
By applying the direct GL, TR and SI methods on this example, the obtained final time t f , the optimal objective value and CPU times for various values of α and n are reported in Table 2. It is worthwhile to note that problem (46) is a free final time problem and has a nonsmooth solution. As a result, obtaining an accurate solution for this problem is not easily accessible. However, from Table 2, the precision and accuracy of the methods, especially the direct SI method, are concluded.

Example 3: A Bang-bang problem
In this example, we apply our new numerical technique to the following bang-bang problem, which is treated in [64]: determine the state x(t) and control u(t) on the interval t ∈ [0, 2] that minimize the performance index , and the initial conditions The exact solution to this problem is given by It is stressed that, in optimal control theory, the field of bang-bang optimal control problems is a classical topic.
For bang-bang problems, the control function switches abruptly between its bounds. Bang-bang optimal control A C C E P T E D M A N U S C R I P T Error of u(t) problems have received considerable attention, due to the difficulties arising in their numerical solutions. See [74] and references therein. Here, we assess our methods with this bang-bang problem.
The obtained control and state functions, by the direct TR method with n = 100, are plotted in Figure 4.
In addition, the errors of the obtained state and control are plotted in this figure too. It is seen that the control and state functions are accurately approximated by the direct TR method. Moreover, we apply the direct TR method with n = 100 on this example for α = 0.1, 0.2, . . . , 0.9, 1.0. The obtained state functions are plotted in Figure 5. We can see that in the all cases the control functions are bang-bang. In addition, the measured CPU times, optimal value of the objective function and switching times are reported in Table 3. To show the precision of the method, the results with n = 400 are also reported in this table. We note that this problem is bang-bang and the control function is discontinuous and the state functions are nonsmooth. Naturally, the accuracy of a numerical method for such problems is lower in comparison with smooth problems. Nevertheless, according to Table 3, the direct TR method provides solutions with reasonable accuracy for this bang-bang problem.
A C C E P T E D M A N U S C R I P T  A HIV optimal control problem is now considered. The problem is stated in [23] as follows:  Table 4. In [23], the authors used an indirect method to solve this FOCP. There, necessary optimality conditions for the problem are firstly derived, then the problem is solved using an iterative algorithm. Here, we solve the problem by the proposed direct methods.
In Figure 6, we plot the obtained control and state functions, by applying the direct TR method with n = 1000 to the problem for α = 0.90, 0.95, 1.00. Our results are in good agreement with those of [23]. In addition, to report the precision of the methods, the values of the performance index obtained by the direct TR and SI methods are given in Table 5. Based on these results, we observe that, without deriving the optimality conditions, the problem can be solved with a reasonable accuracy.

Conclusion
In this paper, three direct methods based on Grünwald-Letnikov, trapezoidal and Simpson approximation formulas are presented for the numerical solution of a general class of fractional optimal control problems. Optimality conditions are not needed in these methods. Thus, they can be applied to any type of fractional optimal control problem.
The proposed direct methods are illustrated in four academic and practical test problems and the results confirm that our methods are reliable. According to our numerical results, for problems with a smooth solution, the accuracy of the direct Simpson method is superior and for problems with discontinuous or a nonsmooth solution, the accuracy is satisfactory. Moreover, by providing the gradient of the objective function and the Jacobian of the constraints, the CPU time is significantly reduced. We conclude that the direct methods here proposed are simple, reliable, reasonably accurate and fast for solving fractional optimal control problems.
As further research works, we can refer to costate estimation and combining mesh generation techniques with the presented method to improve the accuracy and speed of the methods in solving problems with nonsmooth solutions.
A C C E P T E D M A N U S C R I P T