Numerical Solution of Some Types of Fractional Optimal Control Problems

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.

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 Section 2, 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. Let 𝑥∶[𝑎,𝑏]→ℝ be a function, let 𝛼>0 be a real number, and let 𝑛=⌈𝛼⌉, 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Γ(𝛼)^𝑡_𝑎(𝑡−𝜏)^𝛼−1𝑥(𝜏)𝑑𝜏(leftRLFI_𝑡),𝐼^𝛼_𝑏𝑥(𝑡)=1Γ(𝛼)^𝑏_𝑡(𝜏−𝑡)^𝛼−1𝑥(𝜏)𝑑𝜏rightRLFI. (1)

The left (left RLFD) and right (right RLFD) Riemann-Liouville fractional derivatives are given according to _𝑎𝐷^𝛼_𝑡𝑥(𝑡)=1Γ()𝑑𝑛−𝛼𝑑𝑡^𝑛𝑛𝑡_𝑎(𝑡−𝜏)^{𝑛−𝛼−1}𝑥(𝜏)𝑑𝜏(leftRLFD_𝑡),𝐷^𝛼_𝑏𝑥(𝑡)=(−1)^𝑛Γ()𝑑𝑛−𝛼𝑑𝑡^𝑛𝑛𝑏_𝑡(𝜏−𝑡)^{𝑛−𝛼−1}𝑥(𝜏)𝑑𝜏rightRLFD. (2) Moreover, the left (left CFD) and right (right CFD) Caputo fractional derivatives are defined by means of ^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)=1Γ(𝑛−𝛼)^𝑡_𝑎(𝑡−𝜏)^{𝑛−𝛼−1}𝑥^(𝑛)(𝜏)𝑑𝜏(leftCFD_𝐶𝑡),𝐷^𝛼_𝑏𝑥(𝑡)=(−1)^𝑛Γ(𝑛−𝛼)^𝑏_𝑡(𝜏−𝑡)^{𝑛−𝛼−1}𝑥^(𝑛)(𝜏)𝑑𝜏rightCFD. (3) The relation between the right RLFD and the right CFD is as follows [13]: ^𝐶_𝑡𝐷^𝛼_𝑏𝑥(𝑡)=_𝑡𝐷^𝛼_𝑏𝑥(𝑡)−^𝑛−1_𝑘=0𝑥^(𝑘)(𝑏)Γ()()𝑘−𝛼+1𝑏−𝑡^𝑘−𝛼. (4) Further, it holds ^𝐶₀𝐷^𝛼_𝑡𝑐=0,where𝑐isaconstant,_𝐶0𝐷^𝛼_𝑡𝑡^𝑛=0,for𝑛∈ℕ₀⌈𝛼⌉Γ(),𝑛<𝑛+1Γ()𝑡𝑛+1−𝛼^𝑛−𝛼,for𝑛∈ℕ₀,𝑛≥⌈𝛼⌉, (5) where ℕ₀={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]: 𝑇_𝑛+1(𝑧)=2𝑧𝑇_𝑛(𝑧)−𝑇_𝑛−1(𝑧),𝑇₀(𝑧)=1,𝑇₁(𝑧)=𝑧,𝑛=1,2,…. (6) The analytic form of the Chebyshev polynomials 𝑇_𝑛(𝑧) of degree 𝑛 is as follows: 𝑇_𝑛(𝑧)=^⌊𝑛/2⌋_𝑖=0(−1)^𝑖2^{𝑛−2𝑖−1}𝑛(𝑛−𝑖−1)!(𝑖)!()!𝑧𝑛−2𝑖^𝑛−2𝑖, (7) where ⌊𝑛⌋ denotes the biggest integer less than or equal to 𝑛. The orthogonality condition reads ¹₋₁𝑇_𝑖(𝑧)𝑇_𝑗(𝑧)√1−𝑧²𝜋𝑑𝑧=𝜋,for𝑖=𝑗=0;2,for𝑖=𝑗≠0;0,for𝑖≠𝑗. (8) In order to use these polynomials on the interval [0,𝐿], we use the so-called shifted Chebyshev polynomials by introducing the change of variable 𝑧=(2𝑡/𝐿)−1. The shifted Chebyshev polynomials are defined according to 𝑇^∗_𝑛(𝑡)=𝑇_𝑛2𝑡𝐿−1,where𝑇^∗₀(𝑡)=1𝑇^∗₁(𝑡)=2𝑡𝐿−1. (9) Their analytic form is given by 𝑇^∗_𝑛(𝑡)=𝑛^𝑛_𝑘=0(−1)^𝑛−𝑘2^2𝑘(𝑛+𝑘−1)!(2𝑘)!(𝑛−𝑘)!𝐿^𝑘𝑡^𝑘,𝑛=1,2,…. (10) We note that (10) implies that 𝑇^∗_𝑛(0)=(−1)^𝑛, 𝑇^∗_𝑛(𝐿)=1. Further, it is easy to see that the orthogonality condition reads ^𝐿₀𝑇^∗_𝑗(𝑡)𝑇^∗_𝑘(𝑡)𝑤(𝑡)𝑑𝑡=𝛿_𝑗𝑘ℎ_𝑘, (11) with the weight function 𝑤(𝑡)=1/√𝐿𝑡−𝑡²,ℎ_𝑘=(𝑏_𝑘/2)𝜋,𝑏₀=2,𝑏_𝑘=1, for 𝑘≥1.

A function 𝑦∈𝐿²([0,𝐿]) can be expressed in terms of shifted Chebyshev polynomials as 𝑦(𝑡)=^∞_𝑗=0𝑐_𝑛𝑇^∗_𝑛(𝑡), (12) where the coefficients 𝑐_𝑛 are given by 𝑐_𝑛=1ℎ_𝑛^𝐿₀𝑦(𝑡)𝑇^∗_𝑛(𝑡)𝑤(𝑡)𝑑𝑡,𝑛=0,1,…. (13)

3. Necessary Optimality Conditions

Let 𝛼∈(0,1) and let 𝐿,𝑓∶[𝑎,+∞[×ℝ²→ℝ be two differentiable functions. We consider the following FOCP [8]:minimize𝐽(𝑥,𝑢,𝑇)=^𝑇_𝑎𝐿(𝑡,𝑥(𝑡),𝑢(𝑡))𝑑𝑡, (14a) subject to the dynamical system 𝑀₁𝑥̇(𝑡)+𝑀₂^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)((𝑡)(𝑡,𝑥(𝑎)=𝑓𝑡,𝑥,𝑢))=𝑥_𝑎,𝑥(𝑇)=𝑥_𝑇, (14b)(14c) where 𝑀₁,𝑀₂≠0,𝑇,𝑥_𝑎, and 𝑥_𝑇 are fixed real numbers.

Theorem 2 (see [8]). If (𝑥,𝑢,𝑇) is a minimizer of (14a)–(14c), then there exists an adjoint state 𝜆 for which the triple (𝑥,𝑢,𝜆) satisfies the optimality conditions𝑀₁𝑥̇(𝑡)+𝑀₂^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)=𝜕𝐻((𝑡)(𝑡)(𝑡,𝑀𝜕𝜆𝑡,𝑥,𝑢,𝜆))₁𝜆̇(𝑡)−𝑀_2𝑡𝐷^𝛼_𝑇𝜆(𝑡)=−𝜕𝐻((𝑡)(𝑡)(𝑡,𝜕𝑥𝑡,𝑥,𝑢,𝜆))𝜕𝐻((𝑡)(𝑡)(𝑡𝜕𝑢𝑡,𝑥,𝑢,𝜆))=0, (15a)(15b)(15c)for all 𝑡∈[𝑎,𝑇], where the Hamiltonian 𝐻 is defined by 𝐻(𝑡,𝑥,𝑢,𝜆)=𝐿(𝑡,𝑥,𝑢)+𝜆𝑓(𝑡,𝑥,𝑢). (16)

Remark 3. Under some additional assumptions on the objective functional 𝐿 and the right-hand side 𝑓, for example, convexity of 𝐿 and linearity of 𝑓 in 𝑥 and 𝑢, 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,𝐿]; that is, 𝑡_𝑟=𝐿2−𝐿2cos𝜋𝑟𝑁,𝑟=0,1,…,𝑁. (17)

Clenshaw and Curtis [12] introduced an approximation 𝑦_𝑁 of the function 𝑦. We reformulate it to be used with respect to the shifted Chebyshev polynomials as follows: 𝑦_𝑁(𝑡)=^𝑁_𝑛=0^′′𝑎_𝑛𝑇^∗_𝑛(𝑡),𝑎_𝑛=2𝑁^𝑁_𝑟=0^′′𝑦𝑡_𝑟𝑇^∗_𝑛𝑡_𝑟. (18) 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 derivative ^𝐶₀𝐷^𝛼_𝑡𝑦(𝑡) are given.

Theorem 4 (see [17]). An approximation of the fractional derivative of order 𝛼 in the Caputo sense of the function 𝑦 at 𝑡_𝑠 is given by ^𝐶₀𝐷^𝛼_𝑡𝑦_𝑁𝑡_𝑠≅^𝑁_𝑟=0𝑦𝑡_𝑟𝑑^𝛼_𝑠,𝑟,𝛼>0, (19) where𝑑^𝛼_𝑠,𝑟=4𝜃_𝑟𝑁^𝑁_⌈𝛼⌉^𝑁_𝑛=^𝑛_{𝑗=0⌈𝛼⌉𝑘=}𝑛𝜃_𝑛𝑏_𝑗(−1)^𝑛−𝑘(𝑛+𝑘−1)!Γ(𝑘−𝛼+1/2)𝑇^∗_𝑛𝑡_𝑟𝑇^∗_𝑗𝑡_𝑠𝐿^𝛼Γ()Γ,𝑘+1/2)(𝑛−𝑘!Γ𝑘−𝛼−𝑗+1𝑘−𝛼+𝑗+1 (20)where𝑠,𝑟=0,1,…,𝑁, with 𝜃₀=𝜃_𝑁=1/2,𝜃_𝑖=1 for all 𝑖=1,2,…,𝑁−1.

An upper bound for the error in the approximation of the fractional derivative ^𝐶₀𝐷^𝛼_𝑡 of the function 𝑦 is given as follows.

Theorem 5 (see [18]). Let ^𝐶₀𝐷^𝛼_𝑡𝑦_𝑁(𝑡) be the approximation of the fractional derivative ^𝐶₀𝐷^𝛼_𝑡 of the function 𝑦 as given by (19). Then, it holds that ^𝐶₀𝐷^𝛼_𝑡𝑦(𝑡)−^𝐶₀𝐷^𝛼_𝑡𝑦_𝑁(𝑡)₂≤^𝑁_𝑛=0^′′𝑎_𝑛Ω_𝑛𝐺𝑡^𝑘−𝛼;𝑇^∗₀,…,𝑇^∗_𝑁𝐺𝑇^∗₀,…,𝑇^∗_𝑁^1/2, (21) where Ω_𝑛=^𝑛_{⌈𝛼⌉𝑘=}(−1)^𝑛−𝑘()!12𝑛𝑛+𝑘−1×Γ𝑘−𝛼+2×𝑏_𝑗𝐿^𝛼Γ1𝑘+2(𝑛−𝑘)!Γ𝑘−𝛼−𝑗+1×Γ𝑘−𝛼+𝑗+1⁻¹,𝐺𝑥;𝑦₁,𝑦₂,…,𝑦_𝑛=⟨𝑥,𝑥⟩𝑥,𝑦₁⋯𝑥,𝑦_𝑛𝑦₁,𝑥𝑦₁,𝑦₁⋯𝑦₁,𝑦_𝑛𝑦⋮⋮⋱⋮_𝑛,𝑥𝑦_𝑛,𝑦₁⋯𝑦_𝑛,𝑦_𝑛. (22)

4.2. Approximation of the Right RLFD

Let 𝑓 be a sufficiently smooth function in [0,𝑏] and let 𝐽(𝑠;𝑓) be defined as follows: 𝐽𝑠;𝑓=^𝑏_𝑠(𝑡−𝑠)^−𝛼𝑓^′(𝑡)𝑑𝑡,0<𝑠<𝑏. (23) From (3) and (4), we deduce that _𝑠𝐷^𝛼_𝑏𝑓(𝑠)=𝑓(𝑏)Γ()()1−𝛼𝑏−𝑠^−𝛼+𝐽𝑠;𝑓Γ().1−𝛼 (24) We approximate 𝑓(𝑡),0≤𝑡≤𝑏, by a sum of shifted Chebyshev polynomials 𝑇_𝑘(2𝑡/𝑏−1) according to 𝑓(𝑡)≈𝑝_𝑁(𝑡)=^𝑁_𝑘=0^′′𝑎_𝑘𝑇_𝑘2𝑡𝑏,𝑎−1_𝑘=2𝑁^𝑁_𝑗=0^′′𝑓𝑡_𝑗𝑇_𝑘2𝑡_𝑗𝑏−1, (25) where𝑡_𝑗=(𝑏/2)−(𝑏/2)cos(𝜋𝑗/𝑁),𝑗=0,…,𝑁, and obtain 𝐽𝑠;𝑓≈𝐽𝑠;𝑝_𝑁=^𝑏_𝑠𝑝^′_𝑁(𝑡)(𝑡−𝑠)^−𝛼𝑑𝑡. (26)

Lemma 6. Let 𝑝_𝑁 be the polynomial of degree 𝑁 as given by (25). Then, there exists a polynomial 𝐹_𝑁−1 of degree 𝑁−1 such that ^𝑥_𝑠𝑝^′_𝑁(𝑡)−𝑝^′_𝑁(𝑠)(𝑡−𝑠)^−𝛼=𝐹𝑑𝑡_𝑁−1(𝑥)−𝐹_𝑁−1(𝑠)(𝑥−𝑠)^1−𝛼. (27)

Proof. Let 𝑝^′_𝑁(𝑡)−𝑝^′_𝑁(𝑠) be expanded in a Taylor series at 𝑡=𝑠: 𝑝^′_𝑁(𝑡)−𝑝^′_𝑁(𝑠)=^𝑁−1_𝑘=1𝐴_𝑘(𝑠)(𝑡−𝑠)^𝑘. (28) Then, ^𝑥_𝑠𝑝^′_𝑁(𝑡)−𝑝^′_𝑁(𝑠)(𝑡−𝑠)^−𝛼=𝑑𝑡_{𝑁−1𝑘=1}𝐴_𝑘(𝑠)^𝑥_𝑠(𝑡−𝑠)^𝑘−𝛼=)𝑑𝑡(𝑡−𝑠^{1−𝛼𝑁−1}_𝑘=1𝐴_𝑘(𝑠)(𝑡−𝑠)^𝑘𝑘−𝛼+1^𝑥_𝑠. (29) The assertion follows, if we choose 𝐹_𝑁−1(𝑥)=^𝑁−1_𝑘=0𝐴_𝑘(𝑠)(𝑥−𝑠)^𝑘,𝑘−𝛼+1 (30) with an arbitrary constant 𝐴₀(𝑠).

In view of (27), we have 𝐽𝑠;𝑝_𝑁=^𝑏_𝑠𝑝^′_𝑁(𝑡)(𝑡−𝑠)^−𝛼=𝑝𝑑𝑡_′𝑁(𝑠)1−𝛼+𝐹_𝑁−1(𝑏)−𝐹_𝑁−1(𝑠)(𝑏−𝑠)^1−𝛼. (31) Moreover, _𝑠𝐷^𝛼_𝑏𝑓(𝑠) can be approximated by means of _𝑠𝐷^𝛼_𝑏𝑓(𝑠)≈𝑓(𝑏)Γ()()1−𝛼𝑏−𝑠^−𝛼+𝐽𝑠;𝑝_𝑁Γ().1−𝛼 (32) We express 𝐹_𝑁−1(𝑡) 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 𝑥 yields 𝑝^′_𝑁(𝑥)−𝑝^′_𝑁(𝑠)(𝑥−𝑠)^−𝛼=𝐹^′_𝑁−1(𝑥)(𝑥−𝑠)^1−𝛼+𝐹_𝑁−1(𝑥)−𝐹_𝑁−1(𝑠)(1−𝛼)(𝑥−𝑠)^−𝛼, (33) whence 𝑝^′_𝑁(𝑥)−𝑝^′_𝑁(𝑠)=𝐹^′_𝑁−1(𝑥)(𝑥−𝑠)+𝐹_𝑁−1(𝑥)−𝐹_𝑁−1(𝑠)(1−𝛼). (34) To evaluate 𝐹_𝑁−1(𝑠) in (31), we expand 𝐹^′_𝑁−1(𝑥) in terms of the shifted Chebyshev polynomials as 𝐹^′_𝑁−1(𝑥)=^𝑁−2′_𝑘=0𝑏_𝑘𝑇_𝑘2𝑥𝑏−1,0≤𝑥≤𝑏, (35) where the summation symbol with one prime denotes a sum with the first term halved. Integrating both sides of (35) gives 𝐹_𝑁−1(𝑥)−𝐹_𝑁−1(𝑠)=𝑏4^𝑁−1_𝑘=1𝑏_𝑘−1−𝑏_𝑘+1𝑘𝑇_𝑘2𝑥𝑏−1−𝑇_𝑘2𝑠𝑏−1, (36) where 𝑏_𝑁−1=𝑏_𝑁=0. On the other hand, we have (𝑥−𝑠)𝐹^′_𝑁−1(𝑥)=𝑏2𝐹^′_𝑁−1(𝑥)2𝑥𝑏−1−2𝑠𝑏−1. (37) By using the relation 𝑇_𝑘+1(𝑢)+𝑇_𝑘−1(𝑢)=2𝑢𝑇_𝑘(𝑢) and (35), it follows that (𝑥−𝑠)𝐹^′_𝑁−1(𝑥)=𝑏4^𝑁−1′_𝑘=0𝑏_𝑘+1−22𝑠𝑏−1𝑏_𝑘+𝑏_𝑘−1𝑇_𝑘2𝑥𝑏−1, (38) where 𝑏₋₁=𝑏₁. Let 𝑝^′_𝑁(𝑥)=^𝑁−1′_𝑘=0𝑐_𝑘𝑇_𝑘2𝑥𝑏−1. (39) Inserting 𝐹_𝑁−1(𝑥)−𝐹_𝑁−1(𝑠) and (𝑥−𝑠)𝐹^′_𝑁−1(𝑥) as given by (36) and (38) into (34) and taking (39) into account, we get 1−1−𝛼𝑘𝑏_𝑘+1−22𝑠𝑏−1𝑏_𝑘+1+1−𝛼𝑘𝑏_𝑘−1=4𝑏𝑐_𝑘,1≤𝑘. (40)

The Chebyshev coefficients 𝑐_𝑘 of 𝑝^′_𝑁(𝑥) as given by (39) can be evaluated by integrating (39) and comparing it with (25): 𝑐_𝑘−1=𝑐_𝑘+1+4𝑘𝑏𝑎_𝑘,𝑘=𝑁,𝑁−1,…,1, (41) with starting values 𝑐_𝑁=𝑐_𝑁+1=0, where 𝑎_𝑘 are the Chebyshev coefficients of 𝑝_𝑁(𝑥).

5. Numerical Results

In this section, we develop two algorithms (Algorithms A and B) for the numerical solution of FOCPs and apply them to two illustrative examples.

Example 1. We consider the following FOCP from [8]: min𝐽(𝑥,𝑢)=¹₀(𝑡𝑢(𝑡)−(𝛼+2)𝑥(𝑡))²𝑑𝑡, (42a) subject to the dynamical system 𝑥̇(𝑡)+^𝐶₀𝐷^𝛼_𝑡𝑥(𝑡)=𝑢(𝑡)+𝑡² (42b) and the boundary conditions 𝑥(0)(1)=2=0,𝑥Γ().3+𝛼 (42c) The exact solution is given by (𝑥(𝑡),𝑢(𝑡))=2𝑡^𝛼+2Γ(),𝛼+32𝑡^𝛼+1Γ().𝛼+2 (43)

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.

Step 1. Compute the Hamiltonian 𝐻=(𝑡𝑢(𝑡)−(𝛼+2)𝑥(𝑡))²+𝜆𝑢(𝑡)+𝑡². (44)

Step 2. Derive the necessary optimality conditions from Theorem 2: 𝜆̇(𝑡)−_𝑡𝐷^𝛼₁𝜆(𝑡)=−𝜕𝐻((𝑡)−()𝑥(𝑡,𝑥̇(𝑡)+𝜕𝑥=2𝛼+2)(𝑡𝑢𝛼+2))_𝐶0𝐷^𝛼_𝑡𝑥(𝑡)=𝜕𝐻(𝑡)𝜕𝜆=𝑢+𝑡²,0=𝜕𝐻((𝑡)−()𝑥(𝑡𝜕𝑢=2𝑡𝑡𝑢𝛼+2))+𝜆. (45a)(45b)(45c) Use (45c) in (45a) and (45b) to obtain−𝜆̇(𝑡)+_𝑡𝐷^𝛼₁𝜆(𝑡)=(𝛼+2)𝑡𝜆(𝑡),𝑥̇(𝑡)+^𝐶₀𝐷^𝛼_𝑡𝑥(𝑡)𝜆=−2𝑡²+(𝛼+2)𝑡𝑥(𝑡)+𝑡². (46a)(46b)

Step 3. 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 𝑡_𝑠,𝑠=1,2,…,𝑁−1. This gives −^𝑁_𝑟=0𝑑¹_𝑠,𝑟𝜆𝑡_𝑟+𝜆(1)Γ(1−𝛼)1−𝑡_𝑠^−𝛼+𝐽𝑡_𝑠;𝑝_𝑛Γ()=1−𝛼𝛼+2𝑡_𝑠𝜆𝑡_𝑠, (47)𝑠=1,2,…,𝑁−1, where 𝑑¹_𝑠,𝑟 is defined in (20). The system (47) represents 𝑁−1 algebraic equations which can be solved for the unknown coefficients 𝜆(𝑡₁),𝜆(𝑡₂),…,𝜆(𝑡_𝑁−1). Consequently, it remains to compute the two unknowns 𝜆(𝑡₀),𝜆(𝑡_𝑁). This can be done by using any two points 𝑡_𝑎,𝑡_𝑏∈]0,1[ which differ from the Gauss-Lobatto nodes and satisfy (46a). We end up with two equations in two unknowns: 𝜆̇𝑡_𝑎+_𝑡𝐷^𝛼₁𝜆𝑡_𝑎=𝛼+2𝑡_𝑎𝜆𝑡_𝑎,𝜆̇𝑡_𝑏+_𝑡𝐷^𝛼₁𝜆𝑡_𝑏=𝛼+2𝑡_𝑏𝜆𝑡_𝑏. (48)

Step 3.2. In order to solve (46b) by the Chebyshev expansion method, we use (18) to approximate 𝑥. A collocation scheme is defined by substituting (18), (19), and the computed 𝜆 into (46b) and evaluating the results at the shifted Gauss-Lobatto nodes 𝑡_𝑠,𝑠=1,2,…,𝑁−1. This results in ^𝑁_𝑟=0𝑑¹_𝑠,𝑟𝑥𝑡_𝑟+^𝑁_𝑟=0𝑑^𝛼_𝑠,𝑟𝑥𝑡_𝑟𝜆𝑡=−^𝑠2𝑡²_𝑠+𝛼+2𝑡_𝑠𝑥𝑡_𝑠+𝑡²_𝑠,𝑠=1,2,…,𝑁−1, (49) where 𝑑¹_𝑠,𝑟 and 𝑑^𝛼_𝑠,𝑟 are defined in (20). By using the boundary conditions, we have 𝑥(𝑡₀)=0 and 𝑥(𝑡_𝑁)=2/Γ(3+𝛼). The system (49) represents 𝑁−1 algebraic equations which can be solved for the unknown coefficients 𝑥(𝑡₁),𝑥(𝑡₂),…,𝑥(𝑡_𝑁−1).

Figures 1, 2, 3, and 4 display the exact and approximate state 𝑥 and the exact and approximate control 𝑢 for 𝛼=1/2 and 𝑁=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.

Table 1 contains the maximum errors in the state 𝑥 and in the control 𝑢 for 𝑁=2,𝑁=3, and 𝑁=5.

Table 1 

Maximum errors in the state x and in the control u for different values of N.

Algorithm B. The second algorithm follows the “discretize first, then optimize” approach and proceeds according to the following steps.

Step 1. Substitute (42b) into (42a) to obtain min𝐽=¹₀𝑡𝑥̇(𝑡)+^𝐶₀𝐷^𝛼_𝑡𝑥(𝑡)−𝑡²−(𝛼+2)𝑥(𝑡)²𝑑𝑡. (50)

Step 2. Approximate 𝑥 using the Clenshaw and Curtis formula (18) and approximate the Caputo fractional derivative ^𝐶₀𝐷^𝛼_𝑡𝑥 and 𝑥̇ using (19). Then, (50) takes the form min𝐽=¹₀𝑡^𝑁_𝑟=0𝑑¹_𝑡,𝑟𝑥𝑡_𝑟+^𝑁_𝑟=0𝑑^𝛼_𝑡,𝑟𝑥𝑡_𝑟−𝑡²−(𝛼+2)^𝑁_𝑛=0^′′𝑎_𝑛𝑇^∗_𝑛(𝑡)²𝑑𝑡, (51) where 𝑑^𝛼_𝑡,𝑟 is defined as in (20) replacing 𝑡_𝑠 by 𝑡.

Step 3. Use 𝑡=(1/2)(𝜂+1) to transform (51) to 1min𝐽=2¹₋₁12×𝜂+1_𝑁𝑟=0𝑑¹_𝜂,𝑟𝑥𝜂_𝑟+^𝑁_𝑟=0𝑑^𝛼_𝜂,𝑟𝑥𝜂_𝑟−12𝜂+1²−(𝛼+2)^𝑁_𝑛=0^′′𝑎_𝑛𝑇^∗_𝑛𝜂²𝑑𝜂. (52)

Step 4. Use the Clenshaw and Curtis formula [12] ¹₋₁𝐹𝜂2𝑑𝜂≅𝑚^𝑚𝑚_𝑠=0𝑖=0𝜃_𝑠𝐹𝜂_𝑠𝑇2𝑖+1^∗𝑠𝜂_2𝑖−𝑇^∗_𝑠𝜂_2𝑖+2, (53) where 𝜃₀=𝜃_𝑚=12,𝜃_𝑠𝜂=1∀𝑠=1,2,…,𝑚−1,_𝑖()=cos𝜋𝑖𝑚∀𝑖<𝑚,𝜂_𝑖=−1∀𝑖>𝑚, (54) to approximate the integral (52) as a finite sum of shifted Chebyshev polynomials as follows: =1min𝐽𝑚^𝑚𝑚_𝑠=0𝑖=0𝜃_𝑠×12𝑖+12𝜂_𝑠×+1_𝑁𝑟=0𝑑¹_{𝜂𝑠,𝑟}𝑥𝜂_𝑟+^𝑁_𝑟=0𝑑^𝛼_{𝜂𝑠,𝑟}𝑥𝜂_𝑟−12𝜂_𝑠+1²−(𝛼+2)^𝑁_𝑛=0^′′𝑎_𝑛𝑇^∗_𝑛𝜂_𝑠²×𝑇^∗_𝑠𝜂_2𝑖−𝑇^∗_𝑠𝜂_2𝑖+2. (55)

Step 5. According to the Rayleigh-Ritz method, the critical points of the objective functional (42a) are given by 𝜕𝐽𝜕𝑥𝑡₁=0,𝜕𝐽𝜕𝑥𝑡₂=0,…,𝜕𝐽𝜕𝑥𝑡_𝑁=0, (56) which leads to a system of nonlinear algebraic equations. Solve this system by Newton’s method to obtain 𝑥(𝑡₁),𝑥(𝑡₂),…,𝑥(𝑡_𝑁−1) and use the boundary conditions to get 𝑥(𝑡₀),𝑥(𝑡_𝑁). Then, the pair (𝑥,𝑢) which solves the FOCP has the form𝑥(𝑡)=2𝑁^𝑁𝑁_𝑛=0^′′_𝑟=0^′′𝑥𝑡_𝑟𝑇^∗_𝑛𝑡_𝑟𝑇^∗_𝑛(𝑡),𝑢(𝑡)=𝑥̇(𝑡)+^𝐶₀𝐷^𝛼_𝑡𝑥(𝑡)−𝑡². (57a)(57b)

Figures 5, 6, 7, and 8 display the exact and approximate state 𝑥 and the exact and approximate control 𝑢 for 𝛼=1/2 and 𝑁=𝑚=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.

Table 2 contains the maximum errors in the state 𝑥 and in the control 𝑢 for 𝑁=𝑚=2, 𝑁=𝑚=3, and 𝑁=𝑚=5.

Table 2 

Maximum errors in the state x and in the control u for different values of N.

A comparison of Tables 1 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:min𝐽(𝑥,𝑢)=¹₀(𝑢(𝑡)−𝑥(𝑡))²𝑑𝑡, (58a) subject to the dynamical system 𝑥̇(𝑡)+^𝐶₀𝐷^𝛼_𝑡𝑥(𝑡)=𝑢(𝑡)−𝑥(𝑡)+6𝑡^𝛼+2Γ(𝛼+3)+𝑡³ (58b) and the boundary conditions 𝑥(0)(1)=6=0,𝑥Γ().𝛼+4 (58c) The exact solution is given by (𝑥(𝑡),𝑢(𝑡))=6𝑡^𝛼+3Γ(),𝛼+46𝑡^𝛼+3Γ().𝛼+4 (59) We note that, for Example 2, the optimality conditions stated in Theorem 2 are also sufficient (cf. Remark 3).

Table 3 contains a comparison between the maximum error in the state 𝑥 and in the control 𝑢 for Algorithms A and B.

Table 3 

As opposed to Example 1, 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 by 𝑁-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.

Acknowledgments

R. H. W. Hoppe has been supported by the DFG Priority Programs SPP 1253 and SPP 1506, by the NSF Grants DMS-0914788, DMS-1115658, and by the European Science Foundation within the Networking Programme “OPTPDE.”