# A New Numerical Approach for Solving Fractional Optimal Control Problems with the Caputo–Fabrizio Fractional Operator

Authors: Sara Ghaderi; Sohrab Effati; Aghileh Heydari
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In this study, the necessary optimality conditions are achieved for a class of fractional optimal control problems (FOCPs) involving the Caputo–Fabrizio (CF) fractional derivative. We offer a new numerical method based on the Chebyshev cardinal functions for solving the problem. Our goal was to reduce the original problem into a nonlinear system of algebraic equations. For doing this, we approximate the state and control variables in terms of the Chebyshev cardinal functions. The operational matrices (OMs) of left and right CF fractional integrals are also derived for the Chebyshev cardinal functions. Error estimation and convergence analysis of the method are also introduced. The numerical results illustrate the effectiveness and the rapid convergence of the proposed technique. Due to the high accuracy of the solutions, the suggested method is very useful for numerical techniques in control theory.


---

## Body

## 1. Introduction

Fractional calculus is a generalization of the standard integer order calculus to arbitrary order. It brings more degrees of freedom for differentiation in the modeling of various phenomena that are encountered in many application areas, such as biology, engineering, and various domains of applied sciences [1–3]. However, such flexible differential operators do not have unique definitions. The Grunwald–Letnikov, Riemann–Liouville, and Caputo definitions are examples of commonly used approaches that have been used in various papers.The singularity in the definition of the Caputo derivative causes some difficulties in numerical and analytical methods for fractional calculus. Here, we study the Caputo–Fabrizio (CF) fractional operator that at first comes from the definition of the Caputo fractional derivative [4]. This new derivative has many interesting properties. It can portray substance heterogeneities and configurations with different scales, which cannot be managed with the renowned local theories [5]. CF fractional derivative is considered as a promising differentiation operator to model many problems arising in different fields of science and engineering such as biology, physics, control systems, material science, and fluid dynamics [6–9]. For instance, the Allen–Cahn equation with the CF and Atangana–Baleanu fractional derivative has been compared numerically in [10]. A Fourier transform approach has been done in [11] to obtain the solution of diffusion-convection equation in terms of CF derivative. In addition, the Sumudu transformation was applied in [12] to study the nonlinear Fisher’s reaction-diffusion equation in the CF sense.Motivated by the above discussion, we consider the FOCP with CF fractional derivative in its dynamical system. In recent years, FOCPs begin to attract the attention of many researchers due to the increasing applications in different fields such as engineering and biology [13]. However, to the best of our knowledge, so far there are only a few works devoted in the literature for solving FOCPs involving CF derivative. For instance, in [14] authors have studied a practical example and investigated a numerical method for solving fractional optimal control of a spring-mass-viscodamper system with CF fractional operator. Due to the increasing use of CF fractional derivative in different fields, demands for adding control variables to such systems are also rising. Therefore, investigating precise numerical methods seems crucial for FOCPs. Finding necessary optimality conditions has been always one of the important methods for solving optimal control problems, which is considered as an indirect technique. In this study, we are interested in developing an indirect numerical technique for FOCPs involving CF fractional derivative. For this purpose, for the first time, we investigate the necessary optimality conditions for the FOCPs in terms of CF fractional derivative. In addition, we formulate a new numerical method based on the OM of fractional integrals and some useful properties of the Chebyshev cardinal functions to transform the original problem into nonlinear algebraic equations, which can be solved with mathematical software such as Maple and MATLAB. We analyze the convergence of the proposed numerical approach and obtain the error bound. In one numerical example, we also compare the difference between the Caputo and CF fractional derivative in terms of settling time, which is important in control theory.This work is organized as follows: in the next section, a brief overview of the fractional calculus is given. In Section3, we investigate an indirect approach, and then, the necessary optimality conditions for the FOCP are achieved. In Section 4, our numerical method based on the Chebyshev cardinal function is presented and OM of CF fractional integrals is obtained to transform the FOCP to the system of algebraic equations. In Section 5, some error bounds and the convergence of the method are discussed. Finally, we compare and illustrate some example cases to investigate the impact of CF fractional derivative on the problem.


## 2. Preliminaries

In this part, we state the definitions and theorems of fractional calculus that will be used in the next sections of the paper.


### 2.1. Fractional Calculus

The Caputo definition of fractional derivative is defined as follows:(1)DCa+αft=1Γ1−α∫att−τ−αf′τdτ,0<α<1,t≥a.The kernelt−τ−α in equation (1) causes a singularity at t=τ, which is a drawback of this definition. In 2015, Caputo and Fabrizio defined the following fractional derivative as [15]follows:(2)DCFa+αft=2−αMα21−α∫atexp−α1−αt−τf′τdτ,t≥a,where f is a continuous and differentiable function on C1a,b and Mα is a normalization constant and satisfies the following equation:(3)21−α2−αMα+2α2−αMα=1,and then we have the following equation:(4)Mα=22−α,0<α<1.Equation (2) is named as left CF fractional derivative. The right CF fractional derivative is defined as follows:(5)DCFb−αft=−2−αMα21−α∫tbexp−α1−ατ−tf′τdτ,t≤b.Due to this, we get the following definition of the fractional derivative.Definition 1.
Let0<α<1. The left and right CF fractional derivatives of order 0<α<1 of function f∈H1a,b are given by(6)DCFa+αft=11−α∫atexp−α1−αt−τf′τdτ,t≥a,(7)DCFb−αft=−11−α∫tbexp−α1−ατ−tf′τdτ,t≤b.Definition 2.
The corresponding left and right CF fractional integrals off are defined by(8)ICFa+αft=1−αft+α∫atfsds,(9)ICFb−αft=1−αft+α∫tbfsds.Definition 3.
The left and right new RL fractional derivatives are given by(10)DCFRa+αft=11−αddt∫atexp−α1−αt−τfτdτ,t≥a,DCFRb−αft=−11−αddt∫tbexp−α1−ατ−tfτdτ,t≤b.Remark 1.
The relationship between CF fractional derivatives in Caputo and RL senses is as follows:(11)DCFa+αft=DCFRa+αft−fa1−αexp−α1−αt−a,(12)DCFb−αft=DCFRb−αft−fb1−αexp−α1−αb−t.Now, we state some theorems, which play an important role in the sequel.Theorem 1.
For CF fractional derivative in RL sense, we have the integration by part formula as follows:(13)∫abftDCFRa+αgtdt=∫abgtDCFRb−αftdt.Proof.
See [16].Theorem 2.
For CF fractional derivative in Caputo sense, we have the following equation:(14)∫abftDCFa+αgtdt=∫abgtDCFb−αftdt−ga1−α∫abgtexp−α1−αt−adt+fb1−α∫abftexp−α1−αb−tdt.Proof.
By substituting equations (11) and (12) into (13), we get the following equation:(15)∫abftDCFa+αgt+ga1−αexp−α1−αt−adt=∫abgtDCFb−αft+fb1−αexp−α1−αb−tdt,and then we have the following equation:(16)∫abftDCFa+αgtdt=∫abgtDCFb−αftdt−ga1−α∫abgtexp−α1−αt−adt+fb1−α∫abftexp−α1−αb−tdt.Theorem 3.
Let0<α<1. For CF fractional derivative, we have the following equation:(17)ICFa+αDCFa+αft=ft−fa,ICFb−αDCFb−αgt=gt−gb.Proof.
See [16].For more details about applications of CF fractional derivative, see [17].


## 2.1. Fractional Calculus

The Caputo definition of fractional derivative is defined as follows:(1)DCa+αft=1Γ1−α∫att−τ−αf′τdτ,0<α<1,t≥a.The kernelt−τ−α in equation (1) causes a singularity at t=τ, which is a drawback of this definition. In 2015, Caputo and Fabrizio defined the following fractional derivative as [15]follows:(2)DCFa+αft=2−αMα21−α∫atexp−α1−αt−τf′τdτ,t≥a,where f is a continuous and differentiable function on C1a,b and Mα is a normalization constant and satisfies the following equation:(3)21−α2−αMα+2α2−αMα=1,and then we have the following equation:(4)Mα=22−α,0<α<1.Equation (2) is named as left CF fractional derivative. The right CF fractional derivative is defined as follows:(5)DCFb−αft=−2−αMα21−α∫tbexp−α1−ατ−tf′τdτ,t≤b.Due to this, we get the following definition of the fractional derivative.Definition 1.
Let0<α<1. The left and right CF fractional derivatives of order 0<α<1 of function f∈H1a,b are given by(6)DCFa+αft=11−α∫atexp−α1−αt−τf′τdτ,t≥a,(7)DCFb−αft=−11−α∫tbexp−α1−ατ−tf′τdτ,t≤b.Definition 2.
The corresponding left and right CF fractional integrals off are defined by(8)ICFa+αft=1−αft+α∫atfsds,(9)ICFb−αft=1−αft+α∫tbfsds.Definition 3.
The left and right new RL fractional derivatives are given by(10)DCFRa+αft=11−αddt∫atexp−α1−αt−τfτdτ,t≥a,DCFRb−αft=−11−αddt∫tbexp−α1−ατ−tfτdτ,t≤b.Remark 1.
The relationship between CF fractional derivatives in Caputo and RL senses is as follows:(11)DCFa+αft=DCFRa+αft−fa1−αexp−α1−αt−a,(12)DCFb−αft=DCFRb−αft−fb1−αexp−α1−αb−t.Now, we state some theorems, which play an important role in the sequel.Theorem 1.
For CF fractional derivative in RL sense, we have the integration by part formula as follows:(13)∫abftDCFRa+αgtdt=∫abgtDCFRb−αftdt.Proof.
See [16].Theorem 2.
For CF fractional derivative in Caputo sense, we have the following equation:(14)∫abftDCFa+αgtdt=∫abgtDCFb−αftdt−ga1−α∫abgtexp−α1−αt−adt+fb1−α∫abftexp−α1−αb−tdt.Proof.
By substituting equations (11) and (12) into (13), we get the following equation:(15)∫abftDCFa+αgt+ga1−αexp−α1−αt−adt=∫abgtDCFb−αft+fb1−αexp−α1−αb−tdt,and then we have the following equation:(16)∫abftDCFa+αgtdt=∫abgtDCFb−αftdt−ga1−α∫abgtexp−α1−αt−adt+fb1−α∫abftexp−α1−αb−tdt.Theorem 3.
Let0<α<1. For CF fractional derivative, we have the following equation:(17)ICFa+αDCFa+αft=ft−fa,ICFb−αDCFb−αgt=gt−gb.Proof.
See [16].For more details about applications of CF fractional derivative, see [17].


## 3. Indirect Approach for Solving the FOCPs

Our goal is to solve the following FOCP with fixed final time and free final state:(18)minJx.,u.=∫abGt,xt,utdt,and with the CF fractional dynamical system:(19)DCFa+αxt=Ft,xt,ut,0≤α≤1,and the initial condition is as follows:(20)xa=xa,where xt and ut represent the state and control variables, respectively. The expression DCF0+αxt denotes the left CF fractional derivative defined in equation (6) and G and F are scalar functions.


### 3.1. Necessary Optimality Conditions

In this section, we derive the necessary optimality conditions corresponding to the new FOCP (18)–(20) inspired by the paper [16]. The main advantages of this approach are that we obtain the results only in the sense of left and right CF fractional derivatives.Theorem 4.
Letxt,ut be a minimizer of (18) under the dynamical system (19); then, there exists a function λt for which the triplet xt,ut,λt satisfies the following equation:(21)DCFa+αxt=Ft,xt,ut,(22)DCFb−αλt=∂G∂x+λt∂F∂x,(23)∂G∂u+λt∂F∂u=0,(24)xa=xa,λb=0.Proof.
We follow the traditional approach for finding the necessary conditions. To do this, we use a Lagrange multiplier to adjoin dynamic constraint (19) to performance index (18). Thus, we form the augmented functions as follows:(25)J¯x.,u.=∫abGt,xt,ut+λtFt,xt,ut−DCFa+αxtdt,where λt is the Lagrange multiplier also known as costate or adjoint variable. Taking the first variation of augmented performance index J¯x.,u. in equation (25), we obtain the following equation:(26)δJ¯x.,u.=∫ab∂G∂xδxt+∂G∂uδut+δλtFt,xt,ut−DCFa+αxt+λt∂F∂xδxt+∂F∂uδut−DCFa+αδxtdt.
Using Theorem2, the last integral in equation (26) can be written as follows:(27)∫abλtDCFa+αδxtdt=∫abδxtDCFb−αλxtdt−δxa1−α∫abδxtexp−α1−αt−adt+λb1−α∫abλtexp−α1−αb−tdt.
Becausexa is specified, we have δxa=0, and since xb is not specified, we require λb to be zero. With these assumptions, the identity in equation (27) can be written as follows:(28)∫abλtDCFa+αδxtdt=∫abδxtDCFb−αλtdt.
Using equations (26) and (28), we obtain the following equation:(29)δJ¯x.,u.=∫abδλtFt,xt,ut−DCFa+αxt+δxt∂G∂x+λt∂F∂x−DCFb−αλt+δut∂G∂u+λt∂F∂udt.
Minimization ofJ¯x.,u. (and hence minimization of Jx.,u.) requires that the coefficients of δλt, δxt, and δut in equation (29) are zero. This leads to the system of equations (21)–(23) and completes the proof.


## 3.1. Necessary Optimality Conditions

In this section, we derive the necessary optimality conditions corresponding to the new FOCP (18)–(20) inspired by the paper [16]. The main advantages of this approach are that we obtain the results only in the sense of left and right CF fractional derivatives.Theorem 4.
Letxt,ut be a minimizer of (18) under the dynamical system (19); then, there exists a function λt for which the triplet xt,ut,λt satisfies the following equation:(21)DCFa+αxt=Ft,xt,ut,(22)DCFb−αλt=∂G∂x+λt∂F∂x,(23)∂G∂u+λt∂F∂u=0,(24)xa=xa,λb=0.Proof.
We follow the traditional approach for finding the necessary conditions. To do this, we use a Lagrange multiplier to adjoin dynamic constraint (19) to performance index (18). Thus, we form the augmented functions as follows:(25)J¯x.,u.=∫abGt,xt,ut+λtFt,xt,ut−DCFa+αxtdt,where λt is the Lagrange multiplier also known as costate or adjoint variable. Taking the first variation of augmented performance index J¯x.,u. in equation (25), we obtain the following equation:(26)δJ¯x.,u.=∫ab∂G∂xδxt+∂G∂uδut+δλtFt,xt,ut−DCFa+αxt+λt∂F∂xδxt+∂F∂uδut−DCFa+αδxtdt.
Using Theorem2, the last integral in equation (26) can be written as follows:(27)∫abλtDCFa+αδxtdt=∫abδxtDCFb−αλxtdt−δxa1−α∫abδxtexp−α1−αt−adt+λb1−α∫abλtexp−α1−αb−tdt.
Becausexa is specified, we have δxa=0, and since xb is not specified, we require λb to be zero. With these assumptions, the identity in equation (27) can be written as follows:(28)∫abλtDCFa+αδxtdt=∫abδxtDCFb−αλtdt.
Using equations (26) and (28), we obtain the following equation:(29)δJ¯x.,u.=∫abδλtFt,xt,ut−DCFa+αxt+δxt∂G∂x+λt∂F∂x−DCFb−αλt+δut∂G∂u+λt∂F∂udt.
Minimization ofJ¯x.,u. (and hence minimization of Jx.,u.) requires that the coefficients of δλt, δxt, and δut in equation (29) are zero. This leads to the system of equations (21)–(23) and completes the proof.


## 4. Numerical Treatment of the Problem

### 4.1. Chebyshev Cardinal Functions

Here, we consider the Chebyshev cardinal functions and give some basic of their properties to construct our numerical scheme. The accuracy and efficiency of the Chebyshev cardinal functions have been investigated in several papers such as [18, 19].Definition 4.
The Chebyshev cardinal functions of orderN on the interval −1,1 are defined as follows [18]:(30)ψit=TN+1tTN+1′tit−ti,i=1,2,…,N+1,where TN+1x is the first kind Chebyshev polynomial of order N+1 in −1,1 defined by:(31)TN+1x=cosN+1arccosx,and tj, j=1,2,…,N+1 are the roots of TN+1x, which are as follows:(32)ti=−cos2i−1π2N+1,i=1,2,…,N+1.For the Chebyshev cardinal functions, we have the following important property:(33)ψitj=δij,i,j=1,2,…,N+1,where δij is the Kronecker delta function. The Chebyshev cardinal functions are orthogonal with respect to the weight function wt=1/1−t2 on the interval −1,1, and their orthogonality condition is given by [19]:(34)ψit,ψjtw=∫−11ψitψjtwtdt=πN+1δij.whereDefinition 5.
A functionf∈C−1,1 can be approximated in N+1 terms of the shifted Chebyshev cardinal functions as follows [18]:(35)ft≃∑i=1N+1ftiψit=F⊤Ψt,(36)F=ft1ft2⋯ftN+1⊤,Ψt=ψ1tψ2t⋯ψN+1t⊤.


### 4.2. Operational Matrices of the Caputo–Fabrizio Fractional Integrals

In this part of the study, we obtain the OM of the left and right CF fractional integrals.Lemma 1.
LetΨt be the vector of Chebyshev cardinal functions and Θt=1,t,…,tN⊤. Then, we have the following equation [18]:(37)Ψt=V−1Θt,where N+1×N+1 matrix V is Vandermonde’s matrix and(38)V=11…11t1t2…tNtN+1t12t22…tN2tN+12⋮⋮…⋮⋮t1Nt2N…tNNtN+1N,where tj,j=1,2,…,N+1, is defined in equation (32).Remark 2.
LetV be Vandermonde’s matrix defined in Lemma 1. Then, V−1=B can be specified as follows [18]:(39)Bij=−1j−1∑m1<m2<…<mN−j+2≤N+1mk≠itm1tm2⋯tN−j+2∏k=1k≠jN+1tk−tj.Corollary 1.
Using Lemma1 and Remark 2, the Chebyshev cardinal functions ψit, i=1,2,…,N+1, can be represented by:(40)ψit=∑l=1N+1Biltl−1,where Bil, i,l=1,2,…,N+1, is the components of V−1=B.Theorem 5.
The left CF fractional integral of the vectorΨt can be approximated as follows:(41)ICF−1+αΨt≃∑j=1N+1P1jψjt,∑j=1N+1P2jψjt,…,∑j=1N+1PN+1jψjt=PαΨt,where Pα is the N+1×N+1 left CF fractional operational matrix of integration and(42)Pij=1−α∑l=1N+1Biltjl−1+α∑l=1N+1Biltjll−∑l=1N+1Bil−1ll,i,j=1,2,…,N+1.Proof.
Using equation (8) and Corollary 1, we calculate ICF−1+αψit, i=1,…,N+1, as follows:(43)ICF−1+αψit=1−α∑l=1N+1Biltl−1+α∑l=1N+1Bil∫−1tsl−1ds=1−α∑l=1N+1Biltl−1+α∑l=1N+1Biltll−∑l=1N+1Bil−1ll,and so, we have the following equation:(44)ICF−1+αψit≃∑j=1N+1Pijψjt,where(45)Pij=1−α∑l=1N+1Biltjl−1+α∑l=1N+1Biltjll−∑l=1N+1Bil−1ll.Theorem 6.
For the right CF fractional integral, we can write the following equation:(46)ICF1−αΨt≃∑j=1N+1Q1jψjt,∑j=1N+1Q2jψjt,…,∑j=1N+1QN+1jψjt=QαΨt,where Qα is the N+1×N+1 right CF fractional operational matrix of integration and(47)Qij=1−α∑l=1N+1Biltjl−1+α∑l=1N+1Bill−∑l=1N+1Biltjll,i,j=1,2,…,N+1.Proof.
Using the definition of right CF fractional integral and similar to Theorem5, the proof is straightforward.


### 4.3. Description of the Proposed Method

The state variablext, the control variable ut, and costate variable λt can be approximated in terms of the Chebyshev cardinal functions with unknown coefficients as follows:(48)xt≃X⊤Ψt,(49)ut≃U⊤Ψt,(50)λt≃B⊤Ψt,where(51)X⊤=c1,c2,…,cN+1,U⊤=u1,u2,…,uN+1,B⊤=b1,b2,…,bN+1,are unknown vectors, which should be determined, and Ψt is the vector of Chebyshev cardinal functions.We can approximate functionwt=t by the shifted Chebyshev cardinal functions as follows:(52)wt≃t1,t2,…,tN+1⊤Ψt,where ti, i=1,2,…,N+1, is defined in equation (32).Corollary 2.
Suppose thatW⊤Ψt, X⊤Ψt, and U⊤Ψt are the expansions of wt=t, xt, and ut by the Chebyshev cardinal functions, respectively. Then, for any continuous function ϕ:ℝ×ℝ×ℝ⟶ℝ, we have the following equation:(53)ϕt,xt,ut≃φ⊤Ψt,where(54)φ=ϕt1,xt1,ut1,ϕt2,xt2,ut2,…,ϕtN+1,xtN+1,utN+1⊤,orφ=ϕw1,x1,u1,ϕw2,x2,u2,…,ϕwN+1,xN+1,uN+1⊤.Proof.
By considering equation (35), the proof is clear.Now, consider Theorem4. Using Corollary 2, we can also write the following equation:(55)Ft,xt,ut≃ℱ⊤Ψt,(56)∂G∂x+λt∂F∂x=δ1t,xt,ut,λt≃Δ1⊤Ψt,(57)∂G∂u+λt∂F∂u=δ2t,xt,ut,λt≃Δ2⊤Ψt,where(58)ℱ⊤=Ft1,xt1,ut1,Ft2,xt2,ut2,…,FtN+1,xtN+1,utN+1,Δ1⊤=δ1t1,xt1,ut1,λt1,…,δ1tN+1,xtN+1,utN+1,λtN+1,Δ2⊤=δ2t1,xt1,ut1,λt1,…,δ2tN+1,xtN+1,utN+1,λtN+1.By applying the left and right CF fractional integrals of equations (21) and (22), respectively, we have the following equation:(59)xt−x−1=ICF−1+αFt,xt,ut,λt−λ1=ICF1−α∂G∂x+λt∂F∂x.Now, using equations (48)–(50) and (55)–(57) we have the following equation:(60)X⊤−r⊤−ℱ⊤PαΨt=0,(61)B⊤−s⊤−Δ1⊤QαΨt=0,(62)Δ2⊤Ψt=0,where Pα and Qα are the OM of left and right CF fractional integrals, respectively, and r⊤≃x−1x−1⋯x−1 and s⊤≃λ1λ1⋯λ1. Then, the dynamical system (21)–(23) is converted to the following nonlinear system of algebraic equations:(63)X⊤−r⊤−ℱ⊤Pα=0,B⊤−s⊤−Δ1⊤Qα=0,Δ2⊤=0.Remark 3.
It is worth mentioning that the CF fractional derivative does not satisfy the chain rule. In this case, to transform the intervala,b of our problem into −1,1 (since the Chebyshev cardinal functions were considered on the interval −1,1) the changing variable on DCFa+αxt must be done using the following theorem.Theorem 7.
Lett∈a,b and θ=2t/b−a−b+a/b−a∈−1,1. Then, for CF fractional derivative we have the following equation:(64)DCFa+αxt=λDCF−1+βyθ,where λ=2/21−α+αb−a and β=αb−a/21−α+αb−a∈0,1.Proof.
First, we definext=xb−a/2θ+b+a/2=yθ. Now, using the definition of the left CF fractional derivative (6) we can write the following equation:(65)DCFa+αxt=11−α∫atexp−α1−αt−τdxdtτdτ=11−α∫ab−a/2θ+b+a/2exp−α1−αb−a2θ+b+a2−τ2b−adydθ2τb−a−b+ab−adτ,and by changing variable ξ=2τ/b−a−b+a/b−a, we can write the following equation:(66)DCFa+αxt=11−α∫−1θexp−αb−a21−αθ−ξdydθξdξ=221−α+αb−aDCF−1+αb−a/21−α+αb−ayθ.In the following, we use Theorem7 to transform an arbitrary interval to the interval −1,1.


## 4.1. Chebyshev Cardinal Functions

Here, we consider the Chebyshev cardinal functions and give some basic of their properties to construct our numerical scheme. The accuracy and efficiency of the Chebyshev cardinal functions have been investigated in several papers such as [18, 19].Definition 4.
The Chebyshev cardinal functions of orderN on the interval −1,1 are defined as follows [18]:(30)ψit=TN+1tTN+1′tit−ti,i=1,2,…,N+1,where TN+1x is the first kind Chebyshev polynomial of order N+1 in −1,1 defined by:(31)TN+1x=cosN+1arccosx,and tj, j=1,2,…,N+1 are the roots of TN+1x, which are as follows:(32)ti=−cos2i−1π2N+1,i=1,2,…,N+1.For the Chebyshev cardinal functions, we have the following important property:(33)ψitj=δij,i,j=1,2,…,N+1,where δij is the Kronecker delta function. The Chebyshev cardinal functions are orthogonal with respect to the weight function wt=1/1−t2 on the interval −1,1, and their orthogonality condition is given by [19]:(34)ψit,ψjtw=∫−11ψitψjtwtdt=πN+1δij.whereDefinition 5.
A functionf∈C−1,1 can be approximated in N+1 terms of the shifted Chebyshev cardinal functions as follows [18]:(35)ft≃∑i=1N+1ftiψit=F⊤Ψt,(36)F=ft1ft2⋯ftN+1⊤,Ψt=ψ1tψ2t⋯ψN+1t⊤.


## 4.2. Operational Matrices of the Caputo–Fabrizio Fractional Integrals

In this part of the study, we obtain the OM of the left and right CF fractional integrals.Lemma 1.
LetΨt be the vector of Chebyshev cardinal functions and Θt=1,t,…,tN⊤. Then, we have the following equation [18]:(37)Ψt=V−1Θt,where N+1×N+1 matrix V is Vandermonde’s matrix and(38)V=11…11t1t2…tNtN+1t12t22…tN2tN+12⋮⋮…⋮⋮t1Nt2N…tNNtN+1N,where tj,j=1,2,…,N+1, is defined in equation (32).Remark 2.
LetV be Vandermonde’s matrix defined in Lemma 1. Then, V−1=B can be specified as follows [18]:(39)Bij=−1j−1∑m1<m2<…<mN−j+2≤N+1mk≠itm1tm2⋯tN−j+2∏k=1k≠jN+1tk−tj.Corollary 1.
Using Lemma1 and Remark 2, the Chebyshev cardinal functions ψit, i=1,2,…,N+1, can be represented by:(40)ψit=∑l=1N+1Biltl−1,where Bil, i,l=1,2,…,N+1, is the components of V−1=B.Theorem 5.
The left CF fractional integral of the vectorΨt can be approximated as follows:(41)ICF−1+αΨt≃∑j=1N+1P1jψjt,∑j=1N+1P2jψjt,…,∑j=1N+1PN+1jψjt=PαΨt,where Pα is the N+1×N+1 left CF fractional operational matrix of integration and(42)Pij=1−α∑l=1N+1Biltjl−1+α∑l=1N+1Biltjll−∑l=1N+1Bil−1ll,i,j=1,2,…,N+1.Proof.
Using equation (8) and Corollary 1, we calculate ICF−1+αψit, i=1,…,N+1, as follows:(43)ICF−1+αψit=1−α∑l=1N+1Biltl−1+α∑l=1N+1Bil∫−1tsl−1ds=1−α∑l=1N+1Biltl−1+α∑l=1N+1Biltll−∑l=1N+1Bil−1ll,and so, we have the following equation:(44)ICF−1+αψit≃∑j=1N+1Pijψjt,where(45)Pij=1−α∑l=1N+1Biltjl−1+α∑l=1N+1Biltjll−∑l=1N+1Bil−1ll.Theorem 6.
For the right CF fractional integral, we can write the following equation:(46)ICF1−αΨt≃∑j=1N+1Q1jψjt,∑j=1N+1Q2jψjt,…,∑j=1N+1QN+1jψjt=QαΨt,where Qα is the N+1×N+1 right CF fractional operational matrix of integration and(47)Qij=1−α∑l=1N+1Biltjl−1+α∑l=1N+1Bill−∑l=1N+1Biltjll,i,j=1,2,…,N+1.Proof.
Using the definition of right CF fractional integral and similar to Theorem5, the proof is straightforward.


## 4.3. Description of the Proposed Method

The state variablext, the control variable ut, and costate variable λt can be approximated in terms of the Chebyshev cardinal functions with unknown coefficients as follows:(48)xt≃X⊤Ψt,(49)ut≃U⊤Ψt,(50)λt≃B⊤Ψt,where(51)X⊤=c1,c2,…,cN+1,U⊤=u1,u2,…,uN+1,B⊤=b1,b2,…,bN+1,are unknown vectors, which should be determined, and Ψt is the vector of Chebyshev cardinal functions.We can approximate functionwt=t by the shifted Chebyshev cardinal functions as follows:(52)wt≃t1,t2,…,tN+1⊤Ψt,where ti, i=1,2,…,N+1, is defined in equation (32).Corollary 2.
Suppose thatW⊤Ψt, X⊤Ψt, and U⊤Ψt are the expansions of wt=t, xt, and ut by the Chebyshev cardinal functions, respectively. Then, for any continuous function ϕ:ℝ×ℝ×ℝ⟶ℝ, we have the following equation:(53)ϕt,xt,ut≃φ⊤Ψt,where(54)φ=ϕt1,xt1,ut1,ϕt2,xt2,ut2,…,ϕtN+1,xtN+1,utN+1⊤,orφ=ϕw1,x1,u1,ϕw2,x2,u2,…,ϕwN+1,xN+1,uN+1⊤.Proof.
By considering equation (35), the proof is clear.Now, consider Theorem4. Using Corollary 2, we can also write the following equation:(55)Ft,xt,ut≃ℱ⊤Ψt,(56)∂G∂x+λt∂F∂x=δ1t,xt,ut,λt≃Δ1⊤Ψt,(57)∂G∂u+λt∂F∂u=δ2t,xt,ut,λt≃Δ2⊤Ψt,where(58)ℱ⊤=Ft1,xt1,ut1,Ft2,xt2,ut2,…,FtN+1,xtN+1,utN+1,Δ1⊤=δ1t1,xt1,ut1,λt1,…,δ1tN+1,xtN+1,utN+1,λtN+1,Δ2⊤=δ2t1,xt1,ut1,λt1,…,δ2tN+1,xtN+1,utN+1,λtN+1.By applying the left and right CF fractional integrals of equations (21) and (22), respectively, we have the following equation:(59)xt−x−1=ICF−1+αFt,xt,ut,λt−λ1=ICF1−α∂G∂x+λt∂F∂x.Now, using equations (48)–(50) and (55)–(57) we have the following equation:(60)X⊤−r⊤−ℱ⊤PαΨt=0,(61)B⊤−s⊤−Δ1⊤QαΨt=0,(62)Δ2⊤Ψt=0,where Pα and Qα are the OM of left and right CF fractional integrals, respectively, and r⊤≃x−1x−1⋯x−1 and s⊤≃λ1λ1⋯λ1. Then, the dynamical system (21)–(23) is converted to the following nonlinear system of algebraic equations:(63)X⊤−r⊤−ℱ⊤Pα=0,B⊤−s⊤−Δ1⊤Qα=0,Δ2⊤=0.Remark 3.
It is worth mentioning that the CF fractional derivative does not satisfy the chain rule. In this case, to transform the intervala,b of our problem into −1,1 (since the Chebyshev cardinal functions were considered on the interval −1,1) the changing variable on DCFa+αxt must be done using the following theorem.Theorem 7.
Lett∈a,b and θ=2t/b−a−b+a/b−a∈−1,1. Then, for CF fractional derivative we have the following equation:(64)DCFa+αxt=λDCF−1+βyθ,where λ=2/21−α+αb−a and β=αb−a/21−α+αb−a∈0,1.Proof.
First, we definext=xb−a/2θ+b+a/2=yθ. Now, using the definition of the left CF fractional derivative (6) we can write the following equation:(65)DCFa+αxt=11−α∫atexp−α1−αt−τdxdtτdτ=11−α∫ab−a/2θ+b+a/2exp−α1−αb−a2θ+b+a2−τ2b−adydθ2τb−a−b+ab−adτ,and by changing variable ξ=2τ/b−a−b+a/b−a, we can write the following equation:(66)DCFa+αxt=11−α∫−1θexp−αb−a21−αθ−ξdydθξdξ=221−α+αb−aDCF−1+αb−a/21−α+αb−ayθ.In the following, we use Theorem7 to transform an arbitrary interval to the interval −1,1.


## 5. Error Estimates and Convergence

In this section, we estimate the error for the suggested approximation method. To do so, we need to state some definitions and theorems in the Sobolev space. We set:(67)Hwm−1,1=f∈Lw2−1,1:for0≤l≤m,fl∈Lw2−1,1,the Sobolev space of integer order m in the interval −1,1 [20]. Now, the Sobolev norm of integer order m≥0 for function f is given by:(68)fHwm−1,1=∑i=0mfiLw2−1,121/2.It is convenient to introduce the semi-norm:(69)fHwm;N−1,1=∑i=minm,N+1mfiLw2−1,121/2.Note thatfHwm;N−1,1≤fHwm−1,1, and when m≤N+1, we have the following equation:(70)fHwm;N−1,1=fmLw2−1,12.Theorem 8.
Letℱ∈Hwm;N−1,1, wt=1/1−t2, and ℱN=∑i=1N+1ℱtiψit, where ψit and ti for i=1,2,…,N+1 are defined in equations (30) and (32), respectively. Then, the truncation error ℱ−ℱN, satisfies the inequality [20]:(71)ℱ−ℱNLw2−1,1≤CN−mℱHwm;N−1,1,where C depends on m.In the sequel, we obtain an error estimate for the left and right CF fractional integrals.Theorem 9.
Let0<α<1, and suppose ℱ∈Hwm;N−1,1, then:(72)ICF−1+αℱ−ICF−1+αℱNLw2−1,1≤1+αCN−mℱHwm;N−1,1.Proof.
We have the following equation:(73)ICF−1+αℱ−ICF−1+αℱNLw2−1,1=1−αℱt−ℱNt+α∫−11ℱt−ℱNtdtLw2−1,1≤1−αℱt−ℱNtLw2−1,1+α∫−11ℱt−ℱNtLw2−1,1dt≤1−αCN−mℱHwm;N−1,1+4αCN−mℱHwm;N−1,1=1+3αCN−mℱHwm;N−1,1.Similarly, we can obtain an error bound for the right CF fractional integral.Theorem 10.
Let0<α<1, and suppose ℱ∈Hwm;N−1,1, then:(74)ICF1−αℱ−ICF1−αℱNLw2−1,1≤1+3αCN−mℱHwm;N−1,1.Proof.
In a similar approach with Theorem9, the proof is straightforward.Theorem8 ensures that as N⟶∞, the approximated solutions X⊤Ψt, U⊤Ψt, and B⊤Ψt tend to be xt, ut, and λt, respectively. Now, by considering Theorems 9 and 10 one can easily observe that increasing in the number of Chebyshev cardinal functions, ICF1−αℱ−ICF1−αℱNLw2−1,1 and ICF−1+αℱ−ICF−1+αℱNLw2−1,1 go to zero, and equations (60)–(62) converge to the system of equations (21)–(23). This means that the approximate solutions of our numerical technique tend to be the exact solutions as N⟶∞.


## 6. Numerical Simulations

In this section, three examples are presented to verify and demonstrate the effectiveness and accuracy of the proposed method. All numeric examples are performed via MAPLE 18 to 50 decimal digits on an X64-based PC with Intel (R) CPU core i5, 3.10 GHz with 4.0 GB of RAM.Example 1.
Consider the following FOCP:(75)minJx.,u.=12∫01xt−t2+ut+2e−t+t−22dt,with the CF fractional dynamical system:(76)DCF0+0.5xt=xt+ut,and the initial condition:(77)x0=0.First, by changing variablet=θ/2+1/2, we have the following equation:(78)xθ2+12=yθ,uθ2+12=vθ.Using Theorem7, problems (75)–(95) are converted to the following problem:(79)minJx.,u.=14∫−11yθ−θ2−122+vθ+2e−θ/2−1/2+θ2−322dθ,with the CF fractional dynamical system:(80)DCF−1+1/3yθ=34yθ+vθ,and the initial condition:(81)y−1=0.By applying Theorem4, we have the following equation:(82)DCF−1+1/3yθ=34yθ+vθ,(83)DCF1−1/3λθ=12yθ−θ4−14+34λθ,(84)12vθ+e−θ/2−1/2+θ4−34+34λθ=0,(85)y−1=0,λ1=0.Substitutingvθ=−2e−θ/2−1/2−θ/2+3/2−3/2λθ from equation (84) into (82), we get the following equation:(86)DCF−1+αyθ=34yθ−2e−θ/2−1/2−θ2+32−32λθ,DCF1−αλθ=12yθ−θ4−14+34λθ,y−1=0,λ1=0.Taking the left and right CF fractional integrals, we have the following equation:(87)yθ−y−1=ICF−1+34yθ−32e−θ/2−1/2−3θ8+98−98λθ,(88)λθ−λ1=ICF1−12yθ−θ4−14+34λθ.Now, let:(89)yθ≃C⊤Ψθ,λθ≃B⊤Ψθ,and we also consider:(90)−32e−θ/2−1/2−3θ8+98≃d⊤Ψθ,−θ4−14≃z⊤Ψθ,y−1≃r⊤Ψθ,λ1≃s⊤Ψθ,where coefficient vectors d⊤ and z⊤ can be easily determined by equation (35) and r⊤≃y−1,y−1,…,y−1 and s⊤≃λ1,λ1,…,λ1. Then, equations (87) and (88) reduced into the following algebraic system of equations:(91)C⊤−r⊤−34C⊤Pα−d⊤Pα+98B⊤Pα=0,B⊤−s⊤−12C⊤Qα−z⊤Pα−34B⊤Qα=0.where Pα and Qα are the OM of the left and right CF fractional integrals, respectively. The exact solution for this problem is as follows:(92)xt=t,ut=−2exp−t−t+2.By applying the proposed method and takingN=10, the numerical solutions are shown in Figure 1. In Table 1, the absolute errors of state variable, that is, ex=xt−X⊤Ψt and control variable eu=ut−U⊤Ψt, are reported.Figure 1
Approximate solutions ofx. and u. with N=10 in Example 1. (a) Approximate solution of x.. (b) Approximate solution of u..
(a)(b)Table 1
Absolute errors ofx. and u. with different values of N in Example 1.
tN=5N=10N=15exeuexeuexeu0.12.5171E − 083.2674E − 088.4294E − 0171.2004E − 0163.8349E − 0253.2583E − 0250.24.0582E − 084.4677E − 083.3695E − 0162.7646E − 0164.2936E − 0254.5611E − 0250.32.9863E − 083.2526E − 082.9668E − 0178.3497E − 0173.3422E − 0252.9513E − 0250.48.4724E − 091.1096E − 083.4480E − 0163.7533E − 0161.0585E − 0255.6399E − 0260.51.2954E − 089.4916E − 091.9166E − 0162.0914E − 0164.3513E − 0254.6132E − 0250.62.7922E − 082.3190E − 081.6926E − 0161.4635E − 0162.3562E − 0266.3898E − 0270.73.3397E − 082.7344E − 083.2840E − 0162.9101E − 0164.1856E − 0253.8407E − 0250.82.9131E − 082.2005E − 081.4893E − 0161.0268E − 0163.9179E − 0268.0698E − 0260.91.6980E − 089.2559E − 091.7652E − 0162.1763E − 0164.3259E − 0254.5553E − 0251.02.2699E − 0107.4743E − 093.4333E − 0163.6796E − 0162.8840E − 0263.5987E − 026Example 2.
To compare the proposed method, we consider the following FOCP:(93)minJx.,u.=12∫01x2t+u2tdt,with the CF fractional dynamical system:(94)DCF0+αxt=−xt+ut,and the initial condition:(95)x0=1.The exact solution forα=1 for this problem is as follows [21]:(96)xt=cosh2t+βsinh2t,ut=1+2βcosh2t+2+βsinh2t,where(97)β=−cosh2+2sinh22cosh2+sinh2≃−0.98.By applying the proposed method and takingN=10, the numerical solutions are shown in Figure 2. In Table 2 and Figure 3, the absolute errors of state variable and control variable for α=1 are reported.Figure 2
Approximate solutions ofx. and u. with N=10 for α=1 in Example 2. (a) Approximate solution of x.. (b) Approximate solution of u..
(a)(b)Table 2
Absolute errors ofx. and u. with N=10 in Example 2.
tPresented method withN=10Method of [22] with N=100Method of [23] with N=80exeuexex0.11.718486E − 139.432034E − 145.21E − 063.78968E − 100.21.704877E − 139.417102E − 149.95E − 067.40843E − 100.32.481259E − 131.342161E − 131.438E − 051.06101E − 90.42.133112E − 131.121161E − 131.865E − 051.30656E − 90.51.870411E − 145.911187E − 152.291E − 051.43634E − 90.61.94182E − 131.0692881E − 132.728E − 051.40106E − 90.72.539250E − 131.353746E − 133.191E − 051.14363E − 90.81.936446E − 131.017559E − 133.691E − 055.99654E − 100.91.873796E − 139.83377E − 144.241E − 053.01964E − 101.02.574603E − 131.375833E − 134.856E − 051.63814E − 9Figure 3
Absolute error ofx. and u. with N=10 for α=1 in Example 2. (a) Absolute error of x.. (b) Absolute error of u..
(a)(b)Example 3.
Consider the following FOCP:(98)minJx.,u.=12∫02xt−et2+ut+eαt/α−12dt,subject to the CF fractional dynamical system:(99)DCF0+αxt=xt+ut,0<α<1,and the initial condition:(100)x0=1.For this problem, the exact solutions for0<α<1 are as follows:(101)xt=et,ut=−eαt/α−1.The exact value of cost function isJ=0. Table 3 reports our numerical findings for different values of α and N. The reported results in Table 3 verify that the proposed algorithm converges extremely fast only with a low N. Table 4 shows absolute errors for xt and ut for different values of α and N=15. Simulation curves of ut for α=0.7,0.8,0.9 and xt for α=0.9 with N=15 are plotted in Figure 4.Table 3
Different values ofJ in Example 3.
NJα=0.9α=0.8α=0.752.8096E − 054.9830E − 071.5689E − 08108.1339E − 0102.6408E − 0143.1062E − 018152.0076E − 0153.8304E − 0232.3146E − 029Table 4
Absolute errors ofx. and u. with N=15 in Example 3.
tα=0.9α=0.8α=0.7exeuexeuexeu0.21.0856E − 081.2969E − 086.5329E − 0131.1555E − 0123.5886E − 0166.9671E − 0160.44.9735E − 098.6296E − 096.9455E − 0134.7502E − 0138.1750E − 0161.5201E − 0170.63.0405E − 087.4628E − 092.9706E − 0123.6994E − 0131.9337E − 0155.9303E − 0170.86.8330E − 081.3170E − 088.3804E − 0122.4279E − 0125.8962E − 0152.6579E − 0151.04.8541E − 091.0470E − 097.3583E − 0121.4388E − 0125.5705E − 0152.1262E − 0151.28.4030E − 086.0811E − 099.6432E − 0147.6243E − 0131.9651E − 0166.6001E − 0161.45.9926E − 081.0096E − 087.8574E − 0122.1531E − 0125.6618E − 0152.4884E − 0151.63.6498E − 083.2235E − 096.0358E − 0121.5266E − 0124.7136E − 0152.0706E − 0151.81.9460E − 081.9179E − 091.9826E − 0123.3184E − 0141.2500E − 0171.1647E − 0162.05.5613E − 088.9307E − 097.8508E − 0122.1842E − 0125.8332E − 0152.6252E − 015Figure 4
Approximate solutions ofx. and u. with N=15 in Example 3. (a) Approximate solution of x. for α = 0.9. (b) Approximate solution of u. for α = 0.7, 0.8, and 0.9.
(a)(b)Example 4.
Consider the following FOCP [24]:(102)minJx.,u.=12∫052x2t+u2tdt,subjecttoDCF0+αxt=sinxt+t2ut,0<α<1,and the initial condition:(103)x0=π.Applying the suggested numerical method forN=25, simulation curves of xt and ut for α=0.7,0.8,0.9 are depicted in Figure 5. Table 5 provides a summary of the cost functional values derived by the proposed method for different values of α and N. The reported results in Table 5 indicate that the proposed algorithm converges extremely fast. In control theory, the settling time of a dynamical system is considered an important property. Note that in Figure 5, decreasing the fractional order α leads to decreasing the settling time of the response, and thus, the parameter α can be considered as the performance requirements of a controlled system. Figure 6 compares the approximate results of xt with the Caputo and CF fractional operators. These figures show that applying the CF fractional derivative instead of Caputo provides better results in terms of settling time for all values of 0<α<1.Figure 5
Approximate solutions ofx. and u. with N=25 in Example 4. (a) Approximate solution of u. for α = 0.7, 0.8, and 0.9. (b) Approximate solution of u. for α = 0.7, 0.8, and 0.9.
(a)(b)Table 5
Different values ofJ in Example 4.
NJα=0.9α=0.8α=0.7515.276599214.869237914.56744401013.848876413.482111813.14821791513.777815313.368277812.95346312013.769081913.350187312.92824752513.769091113.350609712.9338765Figure 6
Simulation curves ofxt for α=0.7 and α=0.5 in Example 4 within the Caputo and Caputo–Fabrizio fractional derivatives. (a) Simulation curves of xt for α = 0.7. (b) Simulation curves of xt for α = 0.5.
(a)(b)To end this section, we answer a natural question: what are the advantages of the proposed numerical method? To answer this, we summarize what we have observed from numerical experiments and theoretical results as follows:(i)
As can be seen from the numerical examples, the proposed method with a few Chebyshev cardinal functions has high accurate results.(ii)
The use of Chebyshev cardinal functions for approximation of functions is interesting and simple. These functions are also notable for their rapid convergence for the expansion of functions, and this capability makes them very useful for solving various problems in fractional calculus.(iii)
Instead of solving an optimization problem in direct methods for fractional optimal control problem, we deal with a system of FDEs. We also take the advantages of operational matrices to reduce the complexity of computations.(iv)
Calculating the operational matrices of integral can be also used in other approaches when the basis functions are non-smooth.


## 7. Conclusion

A new numerical approach is constructed for approximating the state and control variables of a class of FOCP involving CF fractional derivative. We investigated the necessary optimality conditions in terms of only left and right CF fractional derivatives. The proposed numerical technique is based on the Chebyshev cardinal functions. We utilized the OM of the left and right CF fractional integrals, to reduce the original problem into a system of algebraic equations. The OM of the left and right CF fractional integrals for the Chebyshev cardinal functions is obtained for the first time. The new OM obtained in this study can be used for other types of FOCPs. In Theorem7, we also obtained a formula for the CF fractional derivative when transferring an interval into −1,1. We have noticed that this point seems to be neglected in some papers on fractional calculus. Our numerical simulations show that the proposed method is efficient and has a significantly high rate of convergence.



---

*Source: 6680319-2022-03-21.xml*