A numerical technique for solving multi-dimensional fractional optimal control problems

ABSTRACT In this article, we use the operation matrix (OM) of Riemann–Liouville fractional integral of the shifted Gegenbauer polynomials with the Lagrange multiplier method to provide efficient numerical solutions to the multi-dimensional fractional optimal control problems. The proposed technique transforms the under consideration problems into sets of nonlinear equations which are easy to solve. Numerical tests including numerical comparisons with some existing methods are introduced to demonstrate the accuracy and efficiency of the suggested technique.


Introduction
The mathematical topic of optimal control problems (OCPs) has received much attentions in the last years because it has a lot of important applications in physics, chemistry, engineering, etc. [1,2]. Recently, fractional calculus has demonstrated its accuracy in modelling many applications in various scientific fields [3][4][5] Fractional optimal control problem (FOCP) is an OCP in which the differential equations governing the dynamics of the system contain at least one fractional derivative operator. Due to the various applications of FOCPs in physics and engineering, they received much attentions from many researchers. these applications include the materials with memory and hereditary effects, dynamical processes containing gas diffusion and heat conduction in fractal porous media. Other applications of FOCPs are given in [1,6]. Because of a lot of FOCPs does not have analytical exact solutions, numerous numerical methods are offered to overcome these problems [7][8][9][10][11][12]. In recent years, various OMs for different polynomials like Chebyshev polynomials [13,14], Legendre polynomials [11,15], Jacobi polynomials [10,16,17], Bernstein polynomials [18] and Lagguerre polynomial [19] have been developed to cover the numerical solutions of different types of fractional differential equations [20][21][22][23]. Our technique depends on the OM of RL fractional integral of shifted Gegenbauer polynomial (SGPs). The Gegenbauer polynomials have many useful properties. They generalize Legendre polynomials L j (t) and Chebyshev polynomials of the first T j (t) and second U j (t) kind, and they are special cases of Jacobi polynomials. The most important characteristic of SGPs is achieving rapid rates of convergence [24][25][26][27]. This encourages us to use these polynomials in the suggested technique as basis functions for finding the approximate solutions to a wide class of FOCPs in the following form: subject to with FOCPs is presented. In Section 6, some illustrative examples are offered. Finally concluding remarks end the paper.

Fractional calculus definitions
Definition 2.1: One of the popular definitions of the fractional integral is the Riemann-Liouville (RL), which get from the relation The operator I ν has properties said in [28], we just recall the next property Definition 2.2: D ν is the RL fractional derivative of order ν is defined by where m is the smallest integer greater than ν.

Shifted Gegenbauer polynomials and their properties
The shifted ultraspherical (Gegenbauer) polynomials C S,j (t), of degree j ∈ Z + , with the associated parameter α > −1 2 are a sequence of real polynomials in the finite domain [0, L]. They are a family of orthogonal polynomials which has many properties: The analytical form of the SGP is given by and S,0 (t) = 1, The orthogonal relation of SGPs is: where ω (α) S (t) is the weight function, and it is even function given from the relation This polynomial recover the shifted Chebyshev polynomial of the first kind T S, S,j (t), and the shifted Chebyshev polynomial of the second kind C The square integrable function y(t) in [0, L] is approximated by SGPs as: where the coefficientsỹ j are getting from The approximation of function y(t) in the vector form is defined by where Y T = [ỹ 0 ,ỹ 1 , . . . ,ỹ N ] is the shifted Gegenbauer coefficient vector, and is the shifted Gegenbauer vector. The q times repeated integration of Gegenbauer vector is where P (q) is called the OM of the integration of φ(t).

Fractional SGOM of integration
In this section, the SGOM of RLFI will be derived.
where t ∈ [0, 1] and P (ν) is called OMFI of order ν in the RL sense, it is an (N + 1) × (N + 1) and is written as where ξ i,j,k is given by . (15) Proof: From relation (7) and by using Equations (3) and (4), we can write The function t k+ν can be written as a series of N+1 terms of Gegenbauer polynomial, wherẽ Now, by employing Equations (16)-(18) we obtain: where ξ i,j,k is given in Equation (15).
Writing the last equation in a vector form gives which finishes the proof.

Error estimation
In the following theorem, the error estimation for the approximated functions will be expressed in terms of Gram determinant [29].
Let y(t) be an arbitrary element of H and y * (t) be the unique best approximation of y(t) out of Y, then , (21) where .

Convergence analysis
Consider the error, E I ν of the operational matrix of RLFI as where is an error vector. From Equation (17), we had approximated t k+ν as N j=0t j C α S,j (t). From the above theorem, we have By using Equation (19), the upper bound of the operational matrix of integration will be The following theorem illustrates that by increasing the number of SGPs the error tends to zero. where For the proof [30].

Application of SGOM for FOCPs
In this section, the SGPs are used with the OM of the RLFI for solving the following FOCPs (1) and (2)

Shifted Gegenbauer approximation
Firstly, approximating D (ν) x j (t) and u(t) by, SGPs, and where X and U are unknown coefficients matrices which are written as ⎛ respectively. By using the relation into Equation (25), we have By using SGOM relation together with Equation (25), we get From Equations (26) and (27), we get where R and Q j are known vectors defined as ⎛ By using Equations (25), (28) and (29) into Equation (1), we get Employing Equations (25), (28) and (29), the dynamic constraints (2) can be approximated as where H j is an N × N matrix. The elements of the matrix H j are getting from the following relation where l = 0, 1, . . . , N and k = 0, 1, . . . , N.
By using Equation (33) into Equation (32), we obtain Secondly, assume that where is the unknown Lagrange multiplier vector, ⎛ By applying the necessary conditions of optimality to Equation (36), we have By using Newton iterative method, this system of NEs can be solved for the unknown coefficients of the vectors X, U and .

Illustrative problems
One-Dimensional Problems Problem 6.1: Consider the following FOCP [21], with the dynamic constraint At ν = 1 the exact solution is By applying the suggested technique to Problem 6.1, the resultant numerical results for the state x(t) and the control u(t) variables are displayed through Figure 1(a,b), respectively, at ν = 0.75, 0.85, 0.95, 1 with the exact solutions for N = 8. We noted that the obtained solutions cover the classical results when the value of the fractional order tends to unity. Also in Tables 1 and 2, the absolute errors of the state x(t) and the control u(t) variables for Problem 6.1 are calculated at various choices of N. It is observed that the efficiency of the proposed method is increased by increasing N. Problem 6.2: Consider the following FOCP [21,32] In Figure 2(a,b), the obtained results of the variables x(t) and u(t) of Problem 6.2 are plotted for different values of ν. In Table 3, comparisons of our numerical results for the minimum values of J of Problem 6.2 with different values of ν at N = 8 together with the results obtained in [13] and [33] are tabulated. Obviously, our estimated results are in a good agreement with the results in [13,33].

Problem 6.3: Consider the following FOCP [14]
Min J =       Table 4 comparisons of obtained results for the minimum values of J of Problem 6.3 at different choices of N together with the results obtained in [14] are tabulated. It is noted that the approximated results obtained by the suggested technique are more accurate than the results in [14].
Two-Dimensional Problems Problem 6.4: Consider the following FOCS [34] min J = 1 2  At ν = 1 the exact solution is      Figure 4(a-c) illustrate the behaviour of state variables x 1 (t), x 2 (t) and control variable u(t), respectively, for N = 8 and ν = 0.5, 0.75, 0.85, 0.95 and 1 with the exact solutions. Tables 5-7 display the comparison of the absolute errors by using our mechanism and the method mention in [35] of the variables x 1 (t), x 2 (t) and u(t) for Problem 6.4 at ν = 1 and N = 3, 4, 6. These results show that the suggested method is more accurate than the results of [35]. This problem was solved in [34] by a different technique. The results shown in Figure 4(a-c) are in a good agreement with the results established in [34]. It is remarkable that we achieved satisfactory numerical results with at last 8 numbers of the SGP while in [34], number of approximations starts in 8 and increases up to 128 are used to obtain satisfactory results. So we can deduce that our numerical technique is less computational than that in [34].

Conclusions
In this paper, a new numerical mechanism has been derived to find the approximate solutions of the multi-dimensional FOCPs, this numerical mechanism depends on the SGOM of RLFI. The SGOM of fractional integration reduced the FOCP into an equivalent integral problem. The properties of the SGPs together with the LMM are used to transform the equivalent functional integral equation problem into an algebraic system of equations, which is easily to solve. The applicability, accuracy and rapidity by using few terms of the SGPs of the proposed mechanism are illustrated by numerical applications and the numerical comparisons with some existing methods in the literatures.