# Numerical Solution of Some Types of Fractional Optimal Control Problems

Authors: Nasser Hassan Sweilam; Tamer Mostafa Al-Ajami; Ronald H. W. Hoppe
Journal: The Scientific World Journal
 (2013)

---

## Abstract

We present two different approaches for the numerical solution of fractional optimal control problems (FOCPs) based on a spectral method using Chebyshev polynomials. The fractional derivative is described in the Caputo sense. The first approach follows the paradigm “optimize first, then discretize” and relies on the approximation of the necessary optimality conditions in terms of the associated Hamiltonian. In the second approach, the state equation is discretized first using the Clenshaw and Curtis scheme for the numerical integration of nonsingular functions followed by the Rayleigh-Ritz method to evaluate both the state and control variables. Two illustrative examples are included to demonstrate the validity and applicability of the suggested approaches.


---

## Body

## 1. Introduction

FOCP refers to the minimization of an objective functional subject to dynamical constraints on the state and the control which have fractional order models. Fractional order models are sometimes more appropriate than conventional integer order models to describe physical systems [1–4]. For example, it has been shown that materials with memory and hereditary effects and dynamical processes including gas diffusion and heat conduction in fractal porous media can be more adequately modeled by fractional order models [5]. Numerical methods for solving FOCPs have been suggested in [6–9].This paper presents two numerical methods for solving some types of FOCPs where fractional derivatives are introduced in the Caputo sense. These numerical methods rely on the spectral method where Chebyshev polynomials are used to approximate the unknown functions. Chebyshev polynomials are widely used in numerical computation [10, 11].For the first numerical method, we follow the approach “optimize first, then discretize” and derive the necessary optimality conditions in terms of the associated Hamiltonian. The necessary optimality conditions give rise to fractional boundary value problems that have left Caputo and right Riemann-Liouville fractional derivatives. We construct an approximation of the right Riemann-Liouville fractional derivatives and solve the fractional boundary value problems by the spectral method. The second method relies on the strategy “discretize first, then optimize.” The Clenshaw and Curtis scheme [12] is used for the discretization of the state equation and the objective functional. The Rayleigh-Ritz method provides the optimality conditions in the discrete regime.The paper is organized as follows: in Section2, some basic notations and preliminaries as well as properties of the shifted Chebyshev polynomials are introduced. Section 3 contains the necessary optimality conditions of the FOCP model. Section 4 is devoted to the approximations of the fractional derivatives. In Section 5, we develop two numerical schemes and present two illustrative examples to demonstrate the validity and applicability of the suggested approaches. Finally, in Section 6, we provide a brief conclusion and some final remarks.


## 2. Basic Notations and Preliminaries

### 2.1. Fractional Derivatives and Integrals

Definition 1.
Letx:[a,b]→ℝ be a function, let α>0 be a real number, and let n=⌈α⌉, where ⌈α⌉ denotes the smallest integer greater than or equal to α. The left (left RLFI) and right (right RLFI) Riemann-Liouville fractional integrals are defined by
(1)Iatαx(t)=1Γ(α)∫at‍(t-τ)α-1x(τ)dτ(leftRLFI),Itbαx(t)=1Γ(α)∫tb‍(τ-t)α-1x(τ)dτ(rightRLFI).The left (left RLFD) and right (right RLFD) Riemann-Liouville fractional derivatives are given according to(2)Datαx(t)=1Γ(n-α)dndtn∫at‍(t-τ)n-α-1x(τ)dτ(leftRLFD),Dtbαx(t)=(-1)nΓ(n-α)dndtn∫tb‍(τ-t)n-α-1x(τ)dτ(rightRLFD).
Moreover, the left (left CFD) and right (right CFD) Caputo fractional derivatives are defined by means of
(3)DaCtαx(t)=1Γ(n-α)∫at‍(t-τ)n-α-1x(n)(τ)dτ(leftCFD),DtCbαx(t)=(-1)nΓ(n-α)∫tb‍(τ-t)n-α-1x(n)(τ)dτ(rightCFD).
The relation between the right RLFD and the right CFD is as follows [13]:
(4)DtCbαx(t)=Dtbαx(t)-∑k=0n-1‍x(k)(b)Γ(k-α+1)(b-t)k-α.
Further, it holds
(5)D0Ctαc=0,wherecisaconstant,D0Ctαtn={0,forn∈ℕ0,n<⌈α⌉Γ(n+1)Γ(n+1-α)tn-α,forn∈ℕ0,n≥⌈α⌉,
where ℕ0={0,1,2,…}. We recall that, for α∈ℕ, the Caputo differential operator coincides with the usual differential operator of integer order. For more details on the fractional derivatives definitions and their properties, we refer the reader to [3, 8, 14, 15].


### 2.2. Shifted Chebyshev Polynomials

The well-known Chebyshev polynomials are defined on the interval[-1,1] and can be determined by the following recurrence formula [16]:
(6)Tn+1(z)=2zTn(z)-Tn-1(z),T0(z)=1,T1(z)=z,n=1,2,….
The analytic form of the Chebyshev polynomials Tn(z) of degree n is as follows:
(7)Tn(z)=∑i=0⌊n/2⌋‍(-1)i2n-2i-1n(n-i-1)!(i)!(n-2i)!zn-2i,
where ⌊n⌋ denotes the biggest integer less than or equal to n. The orthogonality condition reads
(8)∫-11‍Ti(z)Tj(z)1-z2dz={π,fori=j=0;π2,fori=j≠0;0,fori≠j.
In order to use these polynomials on the interval [0,L], we use the so-called shifted Chebyshev polynomials by introducing the change of variable z=(2t/L)-1. The shifted Chebyshev polynomials are defined according to
(9)Tn*(t)=Tn(2tL-1),whereT0*(t)=1T1*(t)=2tL-1.
Their analytic form is given by
(10)Tn*(t)=n∑k=0n‍(-1)n-k22k(n+k-1)!(2k)!(n-k)!Lktk,n=1,2,….
We note that (10) implies that Tn*(0)=(-1)n, Tn*(L)=1. Further, it is easy to see that the orthogonality condition reads
(11)∫0L‍Tj*(t)Tk*(t)w(t)dt=δjkhk,
with the weight function w(t)=1/Lt-t2,hk=(bk/2)π,b0=2,bk=1, for k≥1.A functiony∈L2([0,L]) can be expressed in terms of shifted Chebyshev polynomials as
(12)y(t)=∑j=0∞‍cnTn*(t),
where the coefficients cn are given by
(13)cn=1hn∫0L‍y(t)Tn*(t)w(t)dt,n=0,1,….


## 2.1. Fractional Derivatives and Integrals

Definition 1.
Letx:[a,b]→ℝ be a function, let α>0 be a real number, and let n=⌈α⌉, where ⌈α⌉ denotes the smallest integer greater than or equal to α. The left (left RLFI) and right (right RLFI) Riemann-Liouville fractional integrals are defined by
(1)Iatαx(t)=1Γ(α)∫at‍(t-τ)α-1x(τ)dτ(leftRLFI),Itbαx(t)=1Γ(α)∫tb‍(τ-t)α-1x(τ)dτ(rightRLFI).The left (left RLFD) and right (right RLFD) Riemann-Liouville fractional derivatives are given according to(2)Datαx(t)=1Γ(n-α)dndtn∫at‍(t-τ)n-α-1x(τ)dτ(leftRLFD),Dtbαx(t)=(-1)nΓ(n-α)dndtn∫tb‍(τ-t)n-α-1x(τ)dτ(rightRLFD).
Moreover, the left (left CFD) and right (right CFD) Caputo fractional derivatives are defined by means of
(3)DaCtαx(t)=1Γ(n-α)∫at‍(t-τ)n-α-1x(n)(τ)dτ(leftCFD),DtCbαx(t)=(-1)nΓ(n-α)∫tb‍(τ-t)n-α-1x(n)(τ)dτ(rightCFD).
The relation between the right RLFD and the right CFD is as follows [13]:
(4)DtCbαx(t)=Dtbαx(t)-∑k=0n-1‍x(k)(b)Γ(k-α+1)(b-t)k-α.
Further, it holds
(5)D0Ctαc=0,wherecisaconstant,D0Ctαtn={0,forn∈ℕ0,n<⌈α⌉Γ(n+1)Γ(n+1-α)tn-α,forn∈ℕ0,n≥⌈α⌉,
where ℕ0={0,1,2,…}. We recall that, for α∈ℕ, the Caputo differential operator coincides with the usual differential operator of integer order. For more details on the fractional derivatives definitions and their properties, we refer the reader to [3, 8, 14, 15].


## 2.2. Shifted Chebyshev Polynomials

The well-known Chebyshev polynomials are defined on the interval[-1,1] and can be determined by the following recurrence formula [16]:
(6)Tn+1(z)=2zTn(z)-Tn-1(z),T0(z)=1,T1(z)=z,n=1,2,….
The analytic form of the Chebyshev polynomials Tn(z) of degree n is as follows:
(7)Tn(z)=∑i=0⌊n/2⌋‍(-1)i2n-2i-1n(n-i-1)!(i)!(n-2i)!zn-2i,
where ⌊n⌋ denotes the biggest integer less than or equal to n. The orthogonality condition reads
(8)∫-11‍Ti(z)Tj(z)1-z2dz={π,fori=j=0;π2,fori=j≠0;0,fori≠j.
In order to use these polynomials on the interval [0,L], we use the so-called shifted Chebyshev polynomials by introducing the change of variable z=(2t/L)-1. The shifted Chebyshev polynomials are defined according to
(9)Tn*(t)=Tn(2tL-1),whereT0*(t)=1T1*(t)=2tL-1.
Their analytic form is given by
(10)Tn*(t)=n∑k=0n‍(-1)n-k22k(n+k-1)!(2k)!(n-k)!Lktk,n=1,2,….
We note that (10) implies that Tn*(0)=(-1)n, Tn*(L)=1. Further, it is easy to see that the orthogonality condition reads
(11)∫0L‍Tj*(t)Tk*(t)w(t)dt=δjkhk,
with the weight function w(t)=1/Lt-t2,hk=(bk/2)π,b0=2,bk=1, for k≥1.A functiony∈L2([0,L]) can be expressed in terms of shifted Chebyshev polynomials as
(12)y(t)=∑j=0∞‍cnTn*(t),
where the coefficients cn are given by
(13)cn=1hn∫0L‍y(t)Tn*(t)w(t)dt,n=0,1,….


## 3. Necessary Optimality Conditions

Letα∈(0,1) and let L,f:[a,+∞[×ℝ2→ℝ be two differentiable functions. We consider the following FOCP [8]:(14a)minimizeJ(x,u,T)=∫aT‍L(t,x(t),u(t))dt,
subject to the dynamical system
(14b)M1x˙(t)+M2DaCtαx(t)=f(t,x(t),u(t)),(14c)x(a)=xa,x(T)=xT, where M1,M2≠0,T,xa, and xT are fixed real numbers.Theorem 2 (see [8]).
If(x,u,T) is a minimizer of (14a)–(14c), then there exists an adjoint state λ for which the triple (x,u,λ) satisfies the optimality conditions(15a)M1x˙(t)+M2DaCtαx(t)=∂H∂λ(t,x(t),u(t),λ(t)),(15b)M1λ˙(t)-M2DtTαλ(t)=-∂H∂x(t,x(t),u(t),λ(t)),(15c)∂H∂u(t,x(t),u(t),λ(t))=0,for all t∈[a,T], where the Hamiltonian H is defined by
(16)H(t,x,u,λ)=L(t,x,u)+λf(t,x,u).Remark 3.
Under some additional assumptions on the objective functionalL and the right-hand side f, for example, convexity of L and linearity of f in x and u, the optimality conditions (15a)–(15c) are also sufficient.


## 4. Numerical Approximations

In this section, we provide numerical approximations of the left CFD and the right RLFD using Chebyshev polynomials. We choose the grid points to be the Chebyshev-Gauss-Lobatto points associated with the interval[0,L]; that is,
(17)tr=L2-L2cos(πrN),r=0,1,…,N.Clenshaw and Curtis [12] introduced an approximation yN of the function y. We reformulate it to be used with respect to the shifted Chebyshev polynomials as follows:
(18)yN(t)=∑n=0N‍′′anTn*(t),an=2N∑r=0N‍′′y(tr)Tn*(tr).
Here, the summation symbol with double primes denotes a sum with both first and last terms halved.


### 4.1. Approximation of the Left CFD

In the sequel, some basic results for the approximation of the fractional derivativeD0Ctαy(t) are given.Theorem 4 (see [17]).
An approximation of the fractional derivative of orderα in the Caputo sense of the function y at ts is given by
(19)D0CtαyN(ts)≅∑r=0N‍y(tr)ds,rα,α>0,
where(20)ds,rα=4θrN∑n=⌈α⌉N‍∑j=0N‍∑k=⌈α⌉n‍nθnbj(-1)n-k(n+k-1)!Γ(k-α+1/2)Tn*(tr)Tj*(ts)LαΓ(k+1/2)(n-k)!Γ(k-α-j+1)Γ(k-α+j+1),wheres,r=0,1,…,N, with θ0=θN=1/2,θi=1 for all i=1,2,…,N-1.An upper bound for the error in the approximation of the fractional derivativeD0Ctα of the function y is given as follows.Theorem 5 (see [18]).
LetD0CtαyN(t) be the approximation of the fractional derivative D0Ctα of the function y as given by (19). Then, it holds that
(21)∥D0Ctαy(t)-D0CtαyN(t)∥2≤∑n=0N‍′′anΩn(G(tk-α;T0*,…,TN*)G(T0*,…,TN*))1/2,
where
(22)Ωn=∑k=⌈α⌉n‍(Γ(k-α+j+1)12)-1(12(-1)n-k2n(n+k-1)!0000.000×Γ(k-α+12))000000×(bjLαΓ(k+12)(n-k)!Γ(k-α-j+1)0000000000×Γ(k-α+j+1)12)-1),G(x;y1,y2,…,yn)=|〈x,x〉〈x,y1〉⋯〈x,yn〉〈y1,x〉〈y1,y1〉⋯〈y1,yn〉⋮⋮⋱⋮〈yn,x〉〈yn,y1〉⋯〈yn,yn〉|.


### 4.2. Approximation of the Right RLFD

Letf be a sufficiently smooth function in [0,b] and let J(s;f) be defined as follows:
(23)J(s;f)=∫sb‍(t-s)-αf′(t)dt,0<s<b.
From (3) and (4), we deduce that
(24)Dsbαf(s)=f(b)Γ(1-α)(b-s)-α+J(s;f)Γ(1-α).
We approximate f(t),0≤t≤b, by a sum of shifted Chebyshev polynomials Tk(2t/b-1) according to
(25)f(t)≈pN(t)=∑k=0N‍′′akTk(2tb-1),ak=2N∑j=0N‍′′f(tj)Tk(2tjb-1),
wheretj=(b/2)-(b/2)cos(πj/N),j=0,…,N, and obtain
(26)J(s;f)≈J(s;pN)=∫sb‍pN′(t)(t-s)-αdt.Lemma 6.
LetpN be the polynomial of degree N as given by (25). Then, there exists a polynomial FN-1 of degree N-1 such that
(27)∫sx‍[pN′(t)-pN′(s)](t-s)-αdt=[FN-1(x)-FN-1(s)](x-s)1-α.Proof.
LetpN′(t)-pN′(s) be expanded in a Taylor series at t=s:
(28)pN′(t)-pN′(s)=∑k=1N-1‍Ak(s)(t-s)k.
Then,
(29)∫sx‍[pN′(t)-pN′(s)](t-s)-αdt=∑k=1N-1‍Ak(s)∫sx‍(t-s)k-αdt=[(t-s)1-α∑k=1N-1‍Ak(s)(t-s)kk-α+1]sx.
The assertion follows, if we choose
(30)FN-1(x)=∑k=0N-1‍Ak(s)(x-s)kk-α+1,
with an arbitrary constant A0(s).In view of (27), we have
(31)J(s;pN)=∫sb‍pN′(t)(t-s)-αdt=[pN′(s)1-α+FN-1(b)-FN-1(s)](b-s)1-α.
Moreover, Dsbαf(s) can be approximated by means of
(32)Dsbαf(s)≈f(b)Γ(1-α)(b-s)-α+J(s;pN)Γ(1-α).
We express FN-1(t) in (31) by a sum of Chebyshev polynomials and provide the recurrence relation satisfied by the Chebyshev coefficients. Differentiating both sides of (27) with respect to x yields
(33){pN′(x)-pN′(s)}(x-s)-α=FN-1′(x)(x-s)1-α+{FN-1(x)-FN-1(s)}(1-α)(x-s)-α,
whence
(34)pN′(x)-pN′(s)=FN-1′(x)(x-s)+{FN-1(x)-FN-1(s)}(1-α).
To evaluate FN-1(s) in (31), we expand FN-1′(x) in terms of the shifted Chebyshev polynomials as
(35)FN-1′(x)=∑k=0N-2‍′bkTk(2xb-1),0≤x≤b,
where the summation symbol with one prime denotes a sum with the first term halved. Integrating both sides of (35) gives
(36)FN-1(x)-FN-1(s)=b4∑k=1N-1‍bk-1-bk+1k{Tk(2xb-1)-Tk(2sb-1)},
where bN-1=bN=0. On the other hand, we have
(37)(x-s)FN-1′(x)=b2FN-1′(x){(2xb-1)-(2sb-1)}.
By using the relation Tk+1(u)+Tk-1(u)=2uTk(u) and (35), it follows that
(38)(x-s)FN-1′(x)=b4∑k=0N-1‍′{bk+1-2(2sb-1)bk+bk-1}Tk(2xb-1),
where b-1=b1. Let
(39)pN′(x)=∑k=0N-1‍′ckTk(2xb-1).
Inserting FN-1(x)-FN-1(s) and (x-s)FN-1′(x) as given by (36) and (38) into (34) and taking (39) into account, we get
(40){1-1-αk}bk+1-2(2sb-1)bk+{1+1-αk}bk-1=4bck,1≤k.The Chebyshev coefficientsck of pN′(x) as given by (39) can be evaluated by integrating (39) and comparing it with (25):
(41)ck-1=ck+1+4kbak,k=N,N-1,…,1,
with starting values cN=cN+1=0, where ak are the Chebyshev coefficients of pN(x).


## 4.1. Approximation of the Left CFD

In the sequel, some basic results for the approximation of the fractional derivativeD0Ctαy(t) are given.Theorem 4 (see [17]).
An approximation of the fractional derivative of orderα in the Caputo sense of the function y at ts is given by
(19)D0CtαyN(ts)≅∑r=0N‍y(tr)ds,rα,α>0,
where(20)ds,rα=4θrN∑n=⌈α⌉N‍∑j=0N‍∑k=⌈α⌉n‍nθnbj(-1)n-k(n+k-1)!Γ(k-α+1/2)Tn*(tr)Tj*(ts)LαΓ(k+1/2)(n-k)!Γ(k-α-j+1)Γ(k-α+j+1),wheres,r=0,1,…,N, with θ0=θN=1/2,θi=1 for all i=1,2,…,N-1.An upper bound for the error in the approximation of the fractional derivativeD0Ctα of the function y is given as follows.Theorem 5 (see [18]).
LetD0CtαyN(t) be the approximation of the fractional derivative D0Ctα of the function y as given by (19). Then, it holds that
(21)∥D0Ctαy(t)-D0CtαyN(t)∥2≤∑n=0N‍′′anΩn(G(tk-α;T0*,…,TN*)G(T0*,…,TN*))1/2,
where
(22)Ωn=∑k=⌈α⌉n‍(Γ(k-α+j+1)12)-1(12(-1)n-k2n(n+k-1)!0000.000×Γ(k-α+12))000000×(bjLαΓ(k+12)(n-k)!Γ(k-α-j+1)0000000000×Γ(k-α+j+1)12)-1),G(x;y1,y2,…,yn)=|〈x,x〉〈x,y1〉⋯〈x,yn〉〈y1,x〉〈y1,y1〉⋯〈y1,yn〉⋮⋮⋱⋮〈yn,x〉〈yn,y1〉⋯〈yn,yn〉|.


## 4.2. Approximation of the Right RLFD

Letf be a sufficiently smooth function in [0,b] and let J(s;f) be defined as follows:
(23)J(s;f)=∫sb‍(t-s)-αf′(t)dt,0<s<b.
From (3) and (4), we deduce that
(24)Dsbαf(s)=f(b)Γ(1-α)(b-s)-α+J(s;f)Γ(1-α).
We approximate f(t),0≤t≤b, by a sum of shifted Chebyshev polynomials Tk(2t/b-1) according to
(25)f(t)≈pN(t)=∑k=0N‍′′akTk(2tb-1),ak=2N∑j=0N‍′′f(tj)Tk(2tjb-1),
wheretj=(b/2)-(b/2)cos(πj/N),j=0,…,N, and obtain
(26)J(s;f)≈J(s;pN)=∫sb‍pN′(t)(t-s)-αdt.Lemma 6.
LetpN be the polynomial of degree N as given by (25). Then, there exists a polynomial FN-1 of degree N-1 such that
(27)∫sx‍[pN′(t)-pN′(s)](t-s)-αdt=[FN-1(x)-FN-1(s)](x-s)1-α.Proof.
LetpN′(t)-pN′(s) be expanded in a Taylor series at t=s:
(28)pN′(t)-pN′(s)=∑k=1N-1‍Ak(s)(t-s)k.
Then,
(29)∫sx‍[pN′(t)-pN′(s)](t-s)-αdt=∑k=1N-1‍Ak(s)∫sx‍(t-s)k-αdt=[(t-s)1-α∑k=1N-1‍Ak(s)(t-s)kk-α+1]sx.
The assertion follows, if we choose
(30)FN-1(x)=∑k=0N-1‍Ak(s)(x-s)kk-α+1,
with an arbitrary constant A0(s).In view of (27), we have
(31)J(s;pN)=∫sb‍pN′(t)(t-s)-αdt=[pN′(s)1-α+FN-1(b)-FN-1(s)](b-s)1-α.
Moreover, Dsbαf(s) can be approximated by means of
(32)Dsbαf(s)≈f(b)Γ(1-α)(b-s)-α+J(s;pN)Γ(1-α).
We express FN-1(t) in (31) by a sum of Chebyshev polynomials and provide the recurrence relation satisfied by the Chebyshev coefficients. Differentiating both sides of (27) with respect to x yields
(33){pN′(x)-pN′(s)}(x-s)-α=FN-1′(x)(x-s)1-α+{FN-1(x)-FN-1(s)}(1-α)(x-s)-α,
whence
(34)pN′(x)-pN′(s)=FN-1′(x)(x-s)+{FN-1(x)-FN-1(s)}(1-α).
To evaluate FN-1(s) in (31), we expand FN-1′(x) in terms of the shifted Chebyshev polynomials as
(35)FN-1′(x)=∑k=0N-2‍′bkTk(2xb-1),0≤x≤b,
where the summation symbol with one prime denotes a sum with the first term halved. Integrating both sides of (35) gives
(36)FN-1(x)-FN-1(s)=b4∑k=1N-1‍bk-1-bk+1k{Tk(2xb-1)-Tk(2sb-1)},
where bN-1=bN=0. On the other hand, we have
(37)(x-s)FN-1′(x)=b2FN-1′(x){(2xb-1)-(2sb-1)}.
By using the relation Tk+1(u)+Tk-1(u)=2uTk(u) and (35), it follows that
(38)(x-s)FN-1′(x)=b4∑k=0N-1‍′{bk+1-2(2sb-1)bk+bk-1}Tk(2xb-1),
where b-1=b1. Let
(39)pN′(x)=∑k=0N-1‍′ckTk(2xb-1).
Inserting FN-1(x)-FN-1(s) and (x-s)FN-1′(x) as given by (36) and (38) into (34) and taking (39) into account, we get
(40){1-1-αk}bk+1-2(2sb-1)bk+{1+1-αk}bk-1=4bck,1≤k.The Chebyshev coefficientsck of pN′(x) as given by (39) can be evaluated by integrating (39) and comparing it with (25):
(41)ck-1=ck+1+4kbak,k=N,N-1,…,1,
with starting values cN=cN+1=0, where ak are the Chebyshev coefficients of pN(x).


## 5. Numerical Results

In this section, we develop two algorithms (AlgorithmsA and B) for the numerical solution of FOCPs and apply them to two illustrative examples.Example 1.
We consider the following FOCP from [8]: (42a)minJ(x,u)=∫01‍(tu(t)-(α+2)x(t))2dt,
subject to the dynamical system
(42b)x˙(t)+D0Ctαx(t)=u(t)+t2
and the boundary conditions
(42c)x(0)=0,x(1)=2Γ(3+α). The exact solution is given by
(43)(x-(t),u-(t))=(2tα+2Γ(α+3),2tα+1Γ(α+2)).Algorithm A. The first algorithm for the solution of (42a)–(42c) follows the “optimize first, then discretize” approach. It is based on the necessary optimality conditions from Theorem 2 and implements the following steps.Step1. Compute the Hamiltonian
(44)H=(tu(t)-(α+2)x(t))2+λ(u(t)+t2).Step2. Derive the necessary optimality conditions from Theorem 2: (45a)λ˙(t)-tD1αλ(t)=-∂H∂x=2(α+2)(tu(t)-(α+2)x(t)),(45b)x˙(t)+0CDtαx(t)=∂H∂λ=u(t)+t2,(45c)0=∂H∂u=2t(tu(t)-(α+2)x(t))+λ. Use (45c) in (45a) and (45b) to obtain(46a)-λ˙(t)+Dt1αλ(t)=(α+2)tλ(t),(46b)x˙(t)+0CDtαx(t)=-λ2t2+(α+2)tx(t)+t2.Step3. By using Chebyshev expansion, get an approximate solution of the coupled system (46a), (46b) under the boundary conditions (42c).Step
3.1
. In order to solve (46a) by the Chebyshev expansion method, use (18) to approximate λ. A collocation scheme is defined by substituting (18), (19), and (32) into (46a) and evaluating the results at the shifted Gauss-Lobatto nodes ts,s=1,2,…,N-1. This gives
(47)-∑r=0N‍ds,r1λ(tr)+λ(1)Γ(1-α)(1-ts)-α+J(ts;pn)Γ(1-α)=α+2tsλ(ts),s=1,2,…,N-1, where ds,r1 is defined in (20). The system (47) represents N-1 algebraic equations which can be solved for the unknown coefficients λ(t1),λ(t2),…,λ(tN-1). Consequently, it remains to compute the two unknowns λ(t0),λ(tN). This can be done by using any two points ta,tb∈]0,1[ which differ from the Gauss-Lobatto nodes and satisfy (46a). We end up with two equations in two unknowns:
(48)λ˙(ta)+Dt1αλ(ta)=α+2taλ(ta),λ˙(tb)+Dt1αλ(tb)=α+2tbλ(tb).Step3.2. In order to solve (46b) by the Chebyshev expansion method, we use (18) to approximate x. A collocation scheme is defined by substituting (18), (19), and the computed λ into (46b) and evaluating the results at the shifted Gauss-Lobatto nodes ts,s=1,2,…,N-1. This results in
(49)∑r=0N‍ds,r1x(tr)+∑r=0N‍ds,rαx(tr)=-λ(ts)2ts2+α+2tsx(ts)+ts2,s=1,2,…,N-1,
where ds,r1 and ds,rα are defined in (20). By using the boundary conditions, we have x(t0)=0 and x(tN)=2/Γ(3+α). The system (49) represents N-1 algebraic equations which can be solved for the unknown coefficients x(t1),x(t2),…,x(tN-1).Figures1, 2, 3, and 4 display the exact and approximate state x and the exact and approximate control u for α=1/2 and N=2,3.Figure 1
Exact and approximate state.Figure 2
Exact and approximate control.Figure 3
Exact and approximate state.Figure 4
Exact and approximate control.Table1 contains the maximum errors in the state x and in the control u for N=2,N=3, and N=5.Table 1
Maximum errors in the statex and in the control u for different values of N.
N
=
2
N
=
3
N
=
5
Max. error inx
3
.
03292
E
-
2
3
.
4641
E
-
3
2
.
6415
E
-
4
Max. error inu
2.12592
E
-
1
4.1878
E
-
2
7.7493
E
-
3Algorithm B. The second algorithm follows the “discretize first, then optimize” approach and proceeds according to the following steps.Step1. Substitute (42b) into (42a) to obtain
(50)minJ=∫01‍(t[x˙(t)+0CDtαx(t)-t2]0000-(α+2)x(t)[x˙(t)+0CDtαx(t)-t2])2dt.Step2. Approximate x using the Clenshaw and Curtis formula (18) and approximate the Caputo fractional derivative D0Ctαx and x˙ using (19). Then, (50) takes the form
(51)minJ=∫01‍(∑n=0N‍′′t[∑r=0N‍dt,r1x(tr)+∑r=0N‍dt,rαx(tr)-t2]00000-(α+2)∑n=0N‍′′anTn*(t))2dt,
where dt,rα is defined as in (20) replacing ts by t.Step3. Use t=(1/2)(η+1) to transform (51) to
(52)minJ=12∫-11(∑n=0N‍′′12(η+1)‍0000000×[∑r=0N‍dη,r1x(ηr)+∑r=0N‍dη,rαx(ηr)0000000000-(12(η+1))2∑r=0N‍]0000000-(α+2)∑n=0N‍′′anTn*(η))2dη.Step4. Use the Clenshaw and Curtis formula [12]
(53)∫-11‍F(η)dη≅2m∑s=0m‍∑i=0m‍θsF(ηs)2i+1[Ts*(η2i)-Ts*(η2i+2)],
where
(54)θ0=θm=12,θs=1∀s=1,2,…,m-1,ηi=cos[(πi)m]∀i<m,ηi=-1∀i>m,
to approximate the integral (52) as a finite sum of shifted Chebyshev polynomials as follows:
(55)minJ=1m∑s=0m‍∑i=0m‍θs2i+100000000000×(∑n=0N‍′′12(ηs+1)000000000000000×[∑r=0N‍dηs,r1x(ηr)+∑r=0N‍dηs,rαx(ηr)0000000000000000000-(12(ηs+1))2∑r=0N‍]000000000000000-(α+2)∑n=0N‍′′anTn*(ηs))2×[Ts*(η2i)-Ts*(η2i+2)].Step5. According to the Rayleigh-Ritz method, the critical points of the objective functional (42a) are given by
(56)∂J∂x(t1)=0,∂J∂x(t2)=0,…,∂J∂x(tN)=0,
which leads to a system of nonlinear algebraic equations. Solve this system by Newton’s method to obtain x(t1),x(t2),…,x(tN-1) and use the boundary conditions to get x(t0),x(tN). Then, the pair (x,u) which solves the FOCP has the form(57a)x(t)=2N∑n=0N‍′′∑r=0N‍′′x(tr)Tn*(tr)Tn*(t),(57b)u(t)=x˙(t)+D0Ctαx(t)-t2.Figures5, 6, 7, and 8 display the exact and approximate state x and the exact and approximate control u for α=1/2 and N=m=2,3.Figure 5
Exact and approximate state.Figure 6
Exact and approximate control.Figure 7
Exact and approximate state.Figure 8
Exact and approximate control.Table2 contains the maximum errors in the state x and in the control u for N=m=2, N=m=3, and N=m=5.Table 2
Maximum errors in the statex and in the control u for different values of N.
N
=
m
=
2
N
=
m
=
3
N
=
m
=
5
Max. error inx
3.03292
E
-
2
3.4641
E
-
3
2.6416
E
-
4
Max. error inu
2.69495
E
-
1
4.8393
E
-
2
8.0532
E
-
3A comparison of Tables1 and 2 reveals that both algorithms yield comparable numerical results which are more accurate than those obtained by the algorithm used in [8].Example 2.
We consider the following linear-quadratic optimal control problem:(58a)minJ(x,u)=∫01‍(u(t)-x(t))2dt,
subject to the dynamical system
(58b)x˙(t)+D0Ctαx(t)=u(t)-x(t)+6tα+2Γ(α+3)+t3
and the boundary conditions
(58c)x(0)=0,x(1)=6Γ(α+4). The exact solution is given by
(59)(x-(t),u-(t))=(6tα+3Γ(α+4),6tα+3Γ(α+4)).
We note that, for Example 2, the optimality conditions stated in Theorem 2 are also sufficient (cf. Remark 3).Table3 contains a comparison between the maximum error in the state x and in the control u for Algorithms A and B.Table 3
Alg. A,N=3
Alg. B,N=m=3
max. error inx
7.6404
E
-
3
1.1943
E
-
2
max. error inu
7.6404
E
-
3
1.6339
E
-
1
Alg. A,N=5
Alg. B,N=m=5
max. error inx
7.8604
E
-
5
1.0304
E
-
4
max. error inu
7.8604
E
-
5
1.0600
E
-
3As opposed to Example1, in this case, Algorithm A performs substantially better than Algorithm B.


## 6. Conclusions

In this paper, we have presented two algorithms for the numerical solution of a wide class of fractional optimal control problems, one based on the “optimize first, then discretize” approach and the other one on the “discretize first, then optimize” strategy. In both algorithms, the solution is approximated byN-term truncated Chebyshev series. Numerical results for two illustrative examples show that the algorithms converge as the number of terms is increased and that the first algorithm is more accurate than the second one.



---

*Source: 306237-2013-12-09.xml*