Numerical solutions of fractional optimal control with Caputo–Katugampola derivative

*Correspondence: abdelhameed_nagy@yahoo.com 2Department of Mathematics, Faculty of Science, Kuwait University, Safat 13060, Kuwait 3Department of Mathematics, Faculty of Science, Benha University, Benha 13518, Egypt Full list of author information is available at the end of the article Abstract In this paper, we present a numerical technique for solving fractional optimal control problems with a fractional derivative called Caputo–Katugampola derivative. This derivative is a generalization of the Caputo fractional derivative. The proposed technique is based on a spectral method using shifted Chebyshev polynomials of the first kind. The Clenshaw and Curtis scheme for the numerical integration and the Rayleigh–Ritz method are used to estimate the state and control variables. Moreover, the error bound of the fractional derivative operator approximation of Caputo–Katugampola is derived. Illustrative examples are provided to show the validity and applicability of the presented technique.


Introduction
In recent years, many engineering and science problems have arisen in fractional differential equations (FDEs). Torvik and Bagley [6] utilized a fractional derivative to depict the behavior of viscoelastic materials. The fractional logistic model with feedback control has been suggested and analyzed in [17]. Characterization and synthesis of the frequencyband complex noninteger differentiator are studied in [26]. A rational approximation technique is used to approximate the fractional order optimal control problems (FOCPs) in [37].
The main idea of FDEs is the order of the derivative is replaced with a noninteger order. Various definitions of a fractional derivative have been presented in the literature such as Riemann-Liouville, Caputo (see [22,23,25,27,29]). Recently, several definitions of fractional operators and generalized fractional derivatives have been presented [5,30,40]. Katugampola proposed a generalization of the Caputo fractional derivative (see [1,19,20]). Furthermore, Abdeljawad et al. evolved them in various research papers [1,18].
Optimal control problems (OCPs) appear widely in various applications such as aerospace, atmospheric reentry, quantum systems, and space shuttle. Bonnard et al. in [9] used the smooth continuation method in optimal control and performed the analysis of a dissipative two-level quantum system. In [16], the homotopy method combined with the shooting method is used to solve an optimal control problem of the atmospheric reentry of a space shuttle [36].
Nowadays, numerous researchers have widespread OCPs to the fractional case where the dynamic of systems involves a fractional derivative term. Fractional optimal control problems (FOCPs) have been widely investigated in the literature. Agrawal [2] provided a general formulation for a class of FOCPs, and this formulation is similar to the classical OCPs. In [24], the authors used a modified hat function to solve a class of FOCPs. FOCPs with free terminal time are considered in [28]. In [39], time-fractional optimal control problems with Caputo-Fabrizio fractional derivative are presented. In all these papers, FOCPs have been defined with respect to different definitions of fractional derivatives such as Riemann-Liouville and Caputo. However, FOCPs with Caputo-Katugampola derivative have not been investigated yet, and to the best of our knowledge, it is the first time that FOCPs with Caputo-Katugampola derivative are studied in the literature. The Caputo-Katugampola derivative is highly influenced by the value of α and another parameter ρ which is beneficial in graphical simulations associated with real data, see [14]. The motivation of this paper is to investigate the nature of a class of FOCPs with various fractional-order values and an additional parameter ρ. Therefore, in this paper, we deal with the numerical solutions of FOCPs with Caputo-Katugampola derivative, and we derive an estimation of the Caputo-Katugampola definition in terms of Chebyshev polynomials.
In order to find the numerical solutions of FOCPs, the spectral method which tries to approximate the unknown functions by means of orthogonal polynomials, namely Chebyshev polynomials, is used ( [8,21]). Furthermore, the Clenshaw and Curtis technique (see [12,35]) has been utilized to discretize the objective function, and the optimality conditions have been obtained by the Rayleigh-Ritz method.
The main aim of this paper is to study the numerical solution for the FOCP in the following form: subject to where C D α,ρ a + is the Caputo-Katugampola derivative operator, Q 1 , Q 2 = 0, t f , x a , x f are fixed real numbers and χ(t, x(t)) = 0, ∀t ∈ [a, t f ].
The rest of the paper is given as follows: In Sect. 2, the properties of shifted Chebyshev polynomials of the first kind and various definitions of the fractional derivatives and integrals are given. In Sect. 3, the Caputo-Katugampola derivative is approximated in terms of Chebyshev polynomials, and the error bound of the Caputo-Katugampola derivative operator is derived. The structure of the proposed numerical scheme is provided in Sect. 4. Section 5 deals with applying the presented scheme to some FOCPs to show its validity and applicability. The conclusions are stated in Sect. 6.

Basic notations and preliminaries
This section provides some definitions of fractional derivatives and integrals that have been presented in the literature and some properties of the Chebyshev polynomials of the first kind.

Fractional derivatives and integrals
In the following, we recall some definitions regarding fractional integrals and derivatives.    [4,7,18,19]) Suppose that y is an integrable function on [a, b], where 0 < a < b < ∞. For 0 < α < 1 and ρ > 0, the Caputo-Katugampola fractional derivative of order α is defined by Also, we give the following theorem that gives an equivalent form of Caputo-Katugampola fractional derivatives in Definition 2.3, when y ∈ C 1 [a, b].

Shifted Chebyshev polynomials of the first kind
The Chebyshev polynomials of the first kind with degree n, n ≥ 0, are defined as follows: These polynomials satisfy the following recurrence relation: with starting values T 0 (x) = 1 and T 1 (x) = x. The explicit analytic form of T n (x) is given by where n/2 denotes the greatest integer less than or equal to n/2. It is known that the Chebyshev polynomials of the first kind satisfy the following orthogonality relations: Some of useful properties of Chebyshev polynomials [33] are the following: The shifted Chebyshev polynomials of the first kind T * n (t) are defined on [0, τ ] as follows: Again, it is known that the explicit analytic form of the shifted Chebyshev polynomials of the first kind is given by and T * n (t) are orthogonal polynomials on the interval [0, τ ] with respect to the weight function 1 3 Numerical approximations can be expanded in terms of shifted Chebyshev polynomials as follows: with c 0 = 2 and c i = 1 for all i ≥ 1.
Let f ∈ C[0, τ ] and N be a positive integer. Denote where f N (t) is the Chebyshev-Gauss-Lobatto (CGL) interpolation of f (t) on [0, τ ] and the CGL points ( [11,34]) It is important to mention that the double-primed summation means halving the first and the last terms.
To compute an approximate solution for a definite integral of a continuous function g : [-1, 1] → R over the interval [-1, 1], we use the Clenshaw and Curtis formula (see [12,35] where ω j , j = 0, 1, . . . , N , are weights and t j , j = 0, 1, . . . , N , are the roots of (1t 2 ) dT N (t) dt . For N even, the weights are: while the weights for N odd are determined by Next, we study the convergence of the proposed technique.

Theorem 2
The series in Eq. (10) converges to the exact solution.
Proof Consider the exact and the truncated series such that M ≥ N are given as follows: Then we have In what follows, we show that f M (t) is a Cauchy sequence in L 2 w [0, 1] and consequently it converges. From Eq. (13), we have Using Bessel's inequality, M i=N +1 |a i | 2 → 0 as M, N → ∞, therefore f M is a Cauchy, and by the completeness of L 2 w , f M converges to ϑ(t) ∈ L 2 w . Now, we prove that f (t) = ϑ(t),

Approximation of the Caputo-Katugampola fractional derivative
In what follows, we give some fundamental results for the fractional derivative C D α,ρ 0 + f (t). First, we present a lemma which helps us in the proofs that follow.
Proof Using the substitution x = 1 2 (1 + cos θ ), we have Integrating by parts, For ν = 0, Moreover, integration by parts for (15) brings higher derivatives of f and shows more cosine terms. It is noticeable that there is a factor m in the denominator of (15). The second integration by parts gives factors m -1 and m + 1, the third provides factors m -2, m, and m + 2, and so on. For simplicity, we replace all various denominators with m -1 at the second differentiation, m -2 at the third differentiation, and so on until mν at the (ν + 1)st differentiation, hence we fulfilled the desired result. Now, we give an estimation of C D α,ρ 0 + f N (t) in terms of Chebyshev polynomials.

Theorem 3 The Caputo-Katugampola fractional derivatives of the function f N (t) defined in Eq. (11) can be obtained as follows:
where the primed summation means that the last term is halved.
Proof Let us assume that y(t) = t n in (5) and take θ = s ρ t ρ , then we get where B(·, ·) is the beta function. Then, from Eq. (8) and Eq. (11), we have In the following theorem, based on the proof of Theorem 2.1 [38], we show an estimation of the error bound of the fractional derivative operator of Caputo-Katugampola.

Theorem 4 Suppose that f satisfies conditions of Lemma
.
The proof is complete. In what follows, we list the main steps to find the solution of Eqs. (1)-(3). It is important to mention that this technique is generally known as a "direct method". The basic steps are given as follows: Step 1: Use Eq. (2) to write Eq. (1) without the control function u.
Step 2: Approximate both the unknown function x and its Caputo-Katugampola fractional derivative by using Eq. (11) and Eq. (16).
Step 3: In order to use the integration on the interval [-1, 1], we take the transformation t = a + t f -a 2 (η + 1).
Step 4: Approximate the integral whose limits are -1 to 1 by using the Clenshaw and Curtis formula given by Eq. (12).
Step 7: Compute the solutions x of the FOPC by where x(t 0 ) and x(t N ) are calculated from the boundary conditions. Moreover, we use Eq.
(2) to find the control function u.

Numerical results
In this section, we apply the proposed numerical approach in Sect. 4 to illustrative examples.

Example 1
Consider a FOCP as follows: The exact solution is provided by , 2t αρ+1 (αρ + 2) . Table 1 shows the maximum absolute errors for the proposed technique and other published results in the literature at ρ = 1 and α = 0.5 with different values of N . It is obvious that the accuracy of the proposed technique is better than that of other methods. In addition, Tables 2, 3, and 4 show the maximum absolute errors E x and E u of the state and control with various choices of α, ρ, and N , respectively. Moreover, the approximate value of the objective function J is given. Figures 1 and 2 elucidate the approximate of the state and control at N = 3, with different values of ρ, α = 0.5 and α = 0.9, respectively. Figure 3 shows the exact and the approximate solutions of x(t) and u(t) using ρ = 1.5 and α = 0.5. Finally, in Table 5, we compute E x and E u at α = 1 and ρ = 1 with various choices of N . It is clear that the numerical results agree with the analytical solutions and the proposed technique gives accurate numerical results.
The exact solution is provided by In Tables 6, 7, and 8, we obtain the computational results of E x , E u , and J with different values of ρ, α when the value of N increases. The approximate state and control with various choices of ρ and α are displayed in Figs. 4 and 5. Furthermore, Fig. 6 displays the approximate solutions of x(t) and u(t) with their exact solutions. It is evident from   4.5864e-7 5.5041e-5 6.7793e-11

Figure 4
The approximate state x(t) (left) and the approximate control u(t) (right) at various values of ρ and α = 0.5, N = 5

Figure 5
The approximate state x(t) (left) and the approximate control u(t) (right) at various values of ρ and α = 0.9, N = 5