A Numerical Scheme to Solve Fractional Optimal Control Problems

Abstract

We review recent results obtained to solve fractional order optimal control problems with free terminal time and a dynamic constraint involving integer and fractional order derivatives. Some particular cases are studied in detail. A numerical scheme is given, based on expansion formulas for the fractional derivative. The efficiency of the method is illustrated through examples.

1. Introduction

In a letter dated September 30, 1695 l’Hôpital posed the question to Leibniz: what would be the derivative of order 𝛼=1/2? Leibniz’s response was “an apparent paradox, from which one day useful consequences will be drawn.” In these words fractional calculus was born. In 1730, based on the formula 𝑑^𝑛𝑥^𝑚𝑑𝑥^𝑛=𝑚(𝑚−1)⋯(𝑚−𝑛+1)𝑥^𝑚−𝑛=Γ(𝑚+1)Γ()𝑥𝑚−𝑛+1^𝑚−𝑛,(1) Euler suggested to use this relationship also for real values of 𝑛. Taking 𝑚=1 and 𝑛=1/2, he obtained 𝑑^1/2𝑥𝑑𝑥^1/2=4𝑥𝜋.(2) Since then, many different approaches have been carried out, trying to find proper definitions for what should be a derivative and an integral of real order. Starting with Cauchy’s formula for an 𝑛-fold integral, ^𝑡_𝑎𝑑𝜏₁^𝜏1_𝑎𝑑𝜏₂⋯^𝜏𝑛−1_𝑎𝑥𝜏_𝑛𝑑𝜏_𝑛=1(𝑛−1)!^𝑡_𝑎(𝑡−𝜏)^𝑛−1𝑥(𝜏)𝑑𝜏,(3) Riemann defined fractional integration as _𝑎𝐼^𝛼_𝑡𝑥(𝑡)=1Γ(𝛼)^𝑡_𝑎(𝑡−𝜏)^𝛼−1𝑥(𝜏)𝑑𝜏.(4) This is nowadays the most common definition for fractional integral. We remark that when the order 𝛼 is an integer, then the fractional integral becomes a multiple integral, recovering by this way the classical case.

We begin with some basic definitions and properties about fractional operators [1, 2]. To avoid too many details, we omit here the conditions that ensure the existence of such fractional operators and the assumptions in which the results given below hold. For an introduction to the fractional variational calculus we refer the reader to [3].

Definition 1. Let 𝑥∶[𝑎,𝑏]→ℝ be a function, 𝛼>0 a real, and 𝑛=[𝛼]+1, where [⋅] denotes the integer part function. The left and right Riemann-Liouville fractional integrals are defined, respectively, by _𝑎𝐼^𝛼_𝑡𝑥(𝑡)=1Γ(𝛼)^𝑡_𝑎(𝑡−𝜏)^𝛼−1𝑥(𝜏)𝑑𝜏,(leftRLFI_𝑡)𝐼^𝛼_𝑏𝑥(𝑡)=1Γ(𝛼)^𝑏_𝑡(𝜏−𝑡)^𝛼−1𝑥(𝜏)𝑑𝜏.rightRLFI The left and right Riemann-Liouville fractional derivatives are defined, respectively, by _𝑎𝐷^𝛼_𝑡𝑥(𝑡)=1Γ()𝑑𝑛−𝛼𝑑𝑡^𝑛𝑛𝑡_𝑎(𝑡−𝜏)^{𝑛−𝛼−1}𝑥(𝜏)𝑑𝜏,(leftRLFD_𝑡)𝐷^𝛼_𝑏𝑥(𝑡)=(−1)^𝑛Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛𝑏_𝑡(𝜏−𝑡)^{𝑛−𝛼−1}𝑥(𝜏)𝑑𝜏.rightRLFD The left and right Caputo fractional derivatives are defined, respectively, by ^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)=1Γ(𝑛−𝛼)^𝑡_𝑎(𝑡−𝜏)^{𝑛−𝛼−1}𝑥^(𝑛)(𝜏)𝑑𝜏,(leftCFD_𝐶𝑡)𝐷^𝛼_𝑏𝑥(𝑡)=(−1)^𝑛Γ(𝑛−𝛼)^𝑏_𝑡(𝜏−𝑡)^{𝑛−𝛼−1}𝑥^(𝑛)(𝜏)𝑑𝜏.rightCFD

We remark that if 𝛼=𝑛 in Definition 1, then we have the usual operators: _𝑎𝐼^𝑛_𝑡𝑥(𝑡)=^𝑡_𝑎𝑑𝜏₁^𝜏1_𝑎𝑑𝜏₂⋯^𝜏𝑛−1_𝑎𝑥𝜏_𝑛𝑑𝜏_𝑛,_𝑡𝐼^𝑛_𝑏𝑥(𝑡)=^𝑏_𝑡𝑑𝜏₁^𝑏_𝜏1𝑑𝜏₂⋯^𝑏_𝜏𝑛−1𝑥𝜏_𝑛𝑑𝜏_𝑛,_𝑎𝐷^𝑛_𝑡𝑥(𝑡)=𝑥^(𝑛)(𝑡),_𝑡𝐷^𝑛_𝑏𝑥(𝑡)=(−1)^𝑛𝑥^(𝑛)(𝑡),_𝐶𝑎𝐷^𝑛_𝑡𝑥(𝑡)=𝑥^(𝑛)(𝑡),^𝐶_𝑡𝐷^𝛼_𝑏𝑥(𝑡)=(−1)^𝑛𝑥^(𝑛)(𝑡).(5)

Some basic properties are useful, namely, a relationship between the Riemann-Liouville and the Caputo fractional derivatives and a fractional integration by parts formula.

Theorem 2. The following conditions hold: (1)^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)=_𝑎𝐷^𝛼_𝑡𝑥(𝑡)−∑^𝑛−1_𝑘=0(𝑥^(𝑘)(𝑎)/Γ(𝑘−𝛼+1))(𝑡−𝑎)^𝑘−𝛼,(2)_𝑎𝐼^𝛼_𝑡𝑎𝐼^𝛽_𝑡𝑥(𝑡)=_𝑎𝐼_𝑡^𝛼+𝛽𝑥(𝑡),(3)^𝐶_𝑎𝐷^𝛼_𝑡𝑎𝐼^𝛼_𝑡𝑥(𝑡)=𝑥(𝑡),(4)_𝑎𝐼^𝛼_𝑡^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)=𝑥(𝑡)−∑^𝑛−1_𝑘=0(𝑥^(𝑘)(𝑎)/𝑘!)(𝑡−𝑎)^𝑘,(5)∫^𝑏_𝑎𝑦(𝑡)⋅^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)𝑑𝑡=∫^𝑏_𝑎𝑥(𝑡)⋅_𝑡𝐷^𝛼_𝑏𝑦(𝑡)𝑑𝑡+∑^𝑛−1_𝑗=0[_𝑡𝐷_𝑏^{𝛼+𝑗−𝑛}𝑦(𝑡)⋅_𝑡𝐷_𝑏^{𝑛−1−𝑗}𝑥(𝑡)]^𝑏_𝑎.

For numerical purposes, one of the most common procedures is to replace the fractional operators by a series that involves integer derivatives only. The usual one is given by _𝑎𝐷^𝛼_𝑡𝑥(𝑡)=^∞_𝑛=0𝛼𝑛(𝑡−𝑎)^𝑛−𝛼Γ()𝑥𝑛+1−𝛼^(𝑛)(𝑡),(6) where 𝛼𝑛=(−1)^𝑛−1𝛼Γ(𝑛−𝛼)Γ(1−𝛼)Γ(𝑛+1).(7) Although very simple to use, it is easy to conclude that in order to have a small error when we approximate_𝑎𝐷^𝛼_𝑡𝑥 by a finite sum up to order 𝑁, we need to consider a large value for 𝑁; that is, we need to consider the set of admissible functions to be 𝐶^𝑁[𝑎,𝑏] which is an important restriction of the set of the space of functions. Recently, in [4], a new expansion formula is given, with the advantage that we only need the first derivative: _𝑎𝐷^𝛼_𝑡𝑥(𝑡)=𝐴(𝛼)(𝑡−𝑎)^−𝛼𝑥(𝑡)+𝐵(𝛼)(𝑡−𝑎)^1−𝛼𝑥̇(𝑡)−^∞_𝑝=2𝐶𝛼,𝑝(𝑡−𝑎)^{1−𝑝−𝛼}𝑉_𝑝(𝑡),(8) where 𝑉_𝑝(𝑡) is the solution of the system 𝑉̇_𝑝(𝑡)=1−𝑝(𝑡−𝑎)^𝑝−2𝑥(𝑡),𝑉_𝑝(𝑎)=0,(9) for 𝑝=2,3,…, and 𝐴, 𝐵, and 𝐶 are given by 𝐴(𝛼)=1Γ(1−𝛼)1+^∞_𝑝=2Γ𝑝−1+𝛼Γ(𝛼)!,𝐵(𝛼)=1𝑝−1Γ(2−𝛼)1+^∞_𝑝=1Γ𝑝−1+𝛼Γ(),𝐶=1𝛼−1𝑝!𝛼,𝑝Γ()Γ()Γ2−𝛼𝛼−1𝑝−1+𝛼!.𝑝−1(10)

We mention the recent papers [5–7], where similar results are proven for fractional integrals and for other types of fractional operators.

2. Necessary and Sufficient Optimality Conditions

Let 𝛼∈(0,1), and let 𝐿,𝑓∶[𝑎,+∞[×ℝ²→ℝ be two differentiable functions and 𝜙∶[𝑎,+∞[×ℝ→ℝ a differentiable function. The fundamental problem, as studied in [8], is the following: minimize𝐽(𝑥,𝑢,𝑇)=^𝑇_𝑎𝐿(𝑡,𝑥(𝑡),𝑢(𝑡))𝑑𝑡+𝜙(𝑇,𝑥(𝑇))(11) subject to the dynamic control system 𝑀𝑥̇(𝑡)+𝑁^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)=𝑓(𝑡,𝑥(𝑡),𝑢(𝑡))(12) and the initial condition 𝑥(𝑎)=𝑥_𝑎,(13) with (𝑀,𝑁)≠(0,0) and 𝑥_𝑎 being a fixed real number. Thus, we are not only interested in finding the optimal state function 𝑥 and the optimal control 𝑢, but also the optimal time 𝑇.

Theorem 3. If (𝑥,𝑢,𝑇) is a minimizer of (11) under the dynamic constraint (12) and the boundary condition (13), then there exists a function 𝜆 for which the triplet (𝑥,𝑢,𝜆) satisfies (i)the Hamiltonian system 𝑀𝜆̇(𝑡)−𝑁_𝑡𝐷^𝛼_𝑇𝜆(𝑡)=−𝜕𝐻((𝑡)(𝑡)(𝑡,𝑀𝑥̇(𝑡)𝜕𝑥𝑡,𝑥,𝑢,𝜆))+𝑁^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)=𝜕𝐻((𝑡)(𝑡)(𝑡𝜕𝜆𝑡,𝑥,𝑢,𝜆))(14) for all 𝑡∈[𝑎,𝑇]; (ii)the stationary condition 𝜕𝐻((𝑡)(𝑡)(𝑡[];𝜕𝑢𝑡,𝑥,𝑢,𝜆))=0,∀𝑡∈𝑎,𝑇(15)(iii) the transversality conditions 𝐻(𝑡,𝑥(𝑡),𝑢(𝑡),𝜆(𝑡))−𝑁𝜆(𝑡)^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)+𝑁𝑥̇(𝑡)_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)+𝜕𝜙((𝑡𝜕𝑡𝑡,𝑥))_𝑡=𝑇(𝑡)=0,𝑀𝜆+𝑁_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)−𝜕𝜙((𝑡𝜕𝑥𝑡,𝑥))_𝑡=𝑇=0;(16)

where the Hamiltonian 𝐻 is defined by𝐻(𝑡,𝑥,𝑢,𝜆)=𝐿(𝑡,𝑥,𝑢)+𝜆𝑓(𝑡,𝑥,𝑢).(17)

This theorem states the general condition that the optimal solution (𝑥,𝑢,𝑇) must fulfill. Next, depending on extra conditions imposed over the final time 𝑇 or in 𝑥(𝑇), new transversality conditions are obtained.

Theorem 4. Let (𝑥,𝑢) be a minimizer of (11) under the dynamic constraint (12) and the boundary condition (13). (i)If  𝑇 is fixed and 𝑥(𝑇) is free, then Theorem 3 holds with the transversality conditions (16) replaced by 𝑀𝜆(𝑡)+𝑁_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)−𝜕𝜙((𝑡𝜕𝑥𝑡,𝑥))_𝑡=𝑇=0.(18)(ii)If 𝑥(𝑇) is fixed and  𝑇 is free, then Theorem 3 holds with the transversality conditions (16) replaced by 𝐻(𝑡,𝑥(𝑡),𝑢(𝑡),𝜆(𝑡))−𝑁𝜆(𝑡)^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)+𝑁𝑥̇(𝑡)_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)+𝜕𝜙((𝑡𝜕𝑡𝑡,𝑥))_𝑡=𝑇=0.(19)(iii)If 𝑇 and 𝑥(𝑇) are fixed, then Theorem 3 holds with no transversality conditions.(iv)If the terminal point 𝑥(𝑇) belongs to a fixed curve, that is, 𝑥(𝑇)=𝛾(𝑇) for some differentiable curve 𝛾, then Theorem 3 holds with the transversality conditions (16) replaced by 𝐻(𝑡,𝑥(𝑡),𝑢(𝑡),𝜆(𝑡))−𝑁𝜆(𝑡)^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)+𝑁𝑥̇(𝑡)_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)+𝜕𝜙((𝑡−𝛾̇(𝑡(𝑡)𝜕𝑡𝑡,𝑥)))𝑀𝜆+𝑁_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)−𝜕𝜙((𝑡𝜕𝑥𝑡,𝑥))_𝑡=𝑇=0.(20)(v)If 𝑇 is fixed and 𝑥(𝑇)≥𝐾 for some fixed 𝐾∈ℝ, then Theorem 3 holds with the transversality conditions (16) replaced by 𝑀𝜆(𝑡)+𝑁_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)−𝜕𝜙((𝑡𝜕𝑥𝑡,𝑥))_𝑡=𝑇(𝑥(𝑇)(𝑡)≤0,−𝐾)𝑀𝜆+𝑁_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)−𝜕𝜙((𝑡𝜕𝑥𝑡,𝑥))_𝑡=𝑇=0.(21)

Numerically, by using approximation (8) up to order 𝐾, we can transform the original problem into the following classical optimal control problem: minimize𝐽̃(𝑥,𝑢,𝑇)=^𝑇_𝑎𝐿(𝑡,𝑥(𝑡),𝑢(𝑡))𝑑𝑡+𝜙(𝑇,𝑥(𝑇))(22) subject to the dynamic constraints 𝑥̇(𝑡)=𝑓(𝑡,𝑥(𝑡),𝑢(𝑡))−𝑁𝐴(𝑡−𝑎)^−𝛼𝑥(𝑡)𝑀+𝑁𝐵(𝑡−𝑎)^1−𝛼+∑^𝐾_𝑝=2𝑁𝐶_𝑝(𝑡−𝑎)^{1−𝑝−𝛼}𝑉_𝑝(𝑡)𝑀+𝑁𝐵(𝑡−𝑎)^1−𝛼,𝑉̇_𝑝(𝑡)=1−𝑝(𝑡−𝑎)^𝑝−2𝑥(𝑡),𝑝=2,…,𝐾(23) and the initial conditions 𝑥(𝑎)=𝑥_𝑎,𝑉_𝑝(𝑎)=0,𝑝=2,…,𝐾.(24)

Theorem 3 can be generalized in the following way. Observe that we have two initial points for the problem, one for the fractional derivative and a second one for the integral of the functional. We now consider a more general approach, where the initial time for the integral is greater than the initial time of the fractional derivative. We impose a boundary condition on 𝑡=𝐴, but similar conditions could be obtained if we considered conditions at 𝑡=𝑎 instead. The problem is formulated as follows. Let 𝛼∈(0,1), and let 𝐿,𝑓∶[𝑎,+∞[×ℝ²→ℝ be two differentiable functions, 𝜙∶[𝑎,+∞[×ℝ→ℝ a differentiable function, and 𝐴>𝑎 a real. We wish to minimize𝐽(𝑥,𝑢,𝑇)=^𝑇_𝐴𝐿(𝑡,𝑥(𝑡),𝑢(𝑡))𝑑𝑡+𝜙(𝑇,𝑥(𝑇))(25) subject to 𝑀𝑥̇(𝑡)+𝑁^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)((𝑡)(𝑡,𝑥(𝐴)=𝑓𝑡,𝑥,𝑢))(26)=𝑥_𝐴,(27) with (𝑀,𝑁)≠(0,0) and 𝑥_𝐴 being a fixed real number.

Theorem 5. If (𝑥,𝑢,𝑇) is a minimizer of (25) under the dynamic constraint (26) and the boundary condition (27), then there exists a function 𝜆 for which the triplet (𝑥,𝑢,𝜆) satisfies(i)the Hamiltonian system 𝑀𝜆̇(𝑡)−𝑁_𝑡𝐷^𝛼_𝑇𝜆(𝑡)=−𝜕𝐻((𝑡)(𝑡)(𝑡,𝑀𝑥̇(𝑡)𝜕𝑥𝑡,𝑥,𝑢,𝜆))+𝑁^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)=𝜕𝐻((𝑡)(𝑡)(𝑡𝜕𝜆𝑡,𝑥,𝑢,𝜆))(28) for all 𝑡∈[𝐴,𝑇] and _𝑡𝐷^𝛼_𝑇𝜆(𝑡)−_𝑡𝐷^𝛼_𝐴𝜆(𝑡)=0(29) for all 𝑡∈[𝑎,𝐴];(ii)the stationary condition 𝜕𝐻((𝑡)(𝑡)(𝑡[];𝜕𝑢𝑡,𝑥,𝑢,𝜆))=0,∀𝑡∈𝐴,𝑇(30)(iii) the transversality conditions 𝐻(𝑡,𝑥(𝑡),𝑢(𝑡),𝜆(𝑡))−𝑁𝜆(𝑡)^𝐶_𝑎𝐷^𝛼_𝑡𝑥(𝑡)+𝑁𝑥̇(𝑡)_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)+𝜕𝜙((𝑡𝜕𝑡𝑡,𝑥))_𝑡=𝑇(𝑡)=0,𝑀𝜆+𝑁_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)−𝜕𝜙((𝑡𝜕𝑥𝑡,𝑥))_𝑡=𝑇=0,_𝑡𝐼_𝑇^1−𝛼𝜆(𝑡)−_𝑡𝐼_𝐴^1−𝛼𝜆(𝑡)_𝑡=𝑎=0;(31)

where the Hamiltonian 𝐻 is defined by𝐻(𝑡,𝑥,𝑢,𝜆)=𝐿(𝑡,𝑥,𝑢)+𝜆𝑓(𝑡,𝑥,𝑢).(32)

We remark that when 𝐴=𝑎, Theorem 5 reduces to Theorem 3.

Under some additional conditions, namely, convexity conditions over 𝐿, 𝑓, and 𝜙, Theorem 3 provides also sufficient conditions to ensure optimal solutions. The result is given in the next theorem.

Theorem 6. Let (𝑥,𝑢,𝜆) be a triplet satisfying the necessary conditions of Theorem 3. Moreover, assume that (1)𝐿 and 𝑓 are convex on 𝑥 and 𝑢, and 𝜙 is convex in 𝑥; (2)𝑇 is fixed; (3)𝜆(𝑡)≥0 for all 𝑡∈[𝑎,𝑇] or 𝑓 is linear in 𝑥 and 𝑢. Then (𝑥,𝑢) is an optimal solution to problem (11)–(13).

3. Numerical Treatment

So far, we have provided a theoretical approach to fractional optimal control problems, which involves solving fractional differential equations. As it is known, solving such equations is in most cases impossible to do, and numerical methods are used to find approximated solutions for the problem (see, e.g., [9, 10]). We describe next, briefly, how formula (8) is deduced and generalized for arbitrary size expansions.

Let 𝑥∈𝐶²[𝑎,𝑏]. Using integration by parts two times, we deduce that _𝑎𝐷^𝛼_𝑡𝑥(𝑡)=𝑥(𝑎)Γ()()1−𝛼𝑡−𝑎^−𝛼+𝑥̇(𝑎)Γ()()2−𝛼𝑡−𝑎^1−𝛼+(𝑡−𝑎)^1−𝛼Γ(2−𝛼)^𝑡_𝑎1−𝜏−𝑎𝑡−𝑎^1−𝛼𝑥̈(𝜏)𝑑𝜏.(33) By the binomial formula, we can rewrite the fractional derivative as _𝑎𝐷^𝛼_𝑡𝑥(𝑡)=𝑥(𝑎)Γ()()1−𝛼𝑡−𝑎^−𝛼+(𝑡−𝑎)^1−𝛼Γ()𝑥̇(𝑎)+()2−𝛼𝑡−𝑎^1−𝛼Γ(2−𝛼)^𝑡_𝑎^∞_𝑝=0Γ𝑝−1+𝛼Γ(𝛼−1)𝑝!𝜏−𝑎𝑡−𝑎^𝑝×𝑥̈(𝜏)𝑑𝜏.(34) Further integration by parts gives _𝑎𝐷^𝛼_𝑡𝑥(𝑡)=𝐴(𝛼)(𝑡−𝑎)^−𝛼𝑥(𝑡)+𝐵(𝛼)(𝑡−𝑎)^1−𝛼𝑥̇(𝑡)−^∞_𝑝=2𝐶𝛼,𝑝(𝑡−𝑎)^{1−𝑝−𝛼}𝑉_𝑝(𝑡),(35) where 𝑉_𝑝(𝑡), 𝐴(𝛼), 𝐵(𝛼), and 𝐶(𝛼,𝑝) are given in (9)-(10). Following similar calculations, we can deduce the next theorem.

Theorem 7. Fix 𝑛∈ℕ and let 𝑥∈𝐶^𝑛[𝑎,𝑏]. Then, _𝑎𝐷^𝛼_𝑡𝑥(𝑡)=1Γ()()1−𝛼𝑡−𝑎^−𝛼𝑥(𝑡)+^𝑛−1_𝑖=1𝐴(𝛼,𝑖)(𝑡−𝑎)^𝑖−𝛼𝑥^(𝑖)(𝑡)+^∞_𝑝=𝑛−Γ𝑝−𝑛+1+𝛼Γ()Γ()!()−𝛼1+𝛼𝑝−𝑛+1𝑡−𝑎^−𝛼𝑥(𝑡)+𝐵𝛼,𝑝(𝑡−𝑎)^{𝑛−1−𝑝−𝛼}𝑉_𝑝(𝑡),(36) where 𝐴()=1𝛼,𝑖Γ(𝑖+1−𝛼)1+^∞_{𝑝=𝑛−𝑖}Γ𝑝−𝑛+1+𝛼Γ()!,𝐵=Γ𝛼−𝑖𝑝−𝑛+𝑖+1𝛼,𝑝𝑝−𝑛+1+𝛼Γ()Γ()!,𝑉−𝛼1+𝛼𝑝−𝑛+1_𝑝(𝑡)=𝑝−𝑛+1^𝑡_𝑎(𝜏−𝑎)^𝑝−𝑛𝑥(𝜏)𝑑𝜏.(37)

The idea is to replace the fractional derivative with such expansions and to consider finite sums only. When we use the approximation _𝑎𝐷^𝛼_𝑡𝑥(𝑡)≈^𝑛−1_𝑖=0𝐴(𝛼,𝑖,𝑁)(𝑡−𝑎)^𝑖−𝛼𝑥^(𝑖)(𝑡)+^𝑁_𝑝=𝑛𝐵𝛼,𝑝(𝑡−𝑎)^{𝑛−1−𝑝−𝛼}𝑉_𝑝(𝑡),(38) the error is bounded by 𝐸_𝑡𝑟(𝑡)≤𝐿_𝑛𝑒^{(𝑛−1−𝛼)2+𝑛−1−𝛼}Γ(𝑛−𝛼)(𝑛−1−𝛼)𝑁^{𝑛−1−𝛼}(𝑡−𝑎)^𝑛−𝛼,(39) where 𝐿_𝑛=max_{[]𝜏∈𝑎,𝑡}𝑥^(𝑛)(𝜏).(40)

To see the accuracy of the method, we exemplify it by considering some functions and compare the exact expression of the fractional derivative with the approximated one. To start, consider 𝑥₁(𝑡)=𝑡⁴ and 𝑥₂(𝑡)=𝑒^2𝑡 and expansions with 𝑛=2 and different values for 𝑁. The result is exemplified in Figures 1(a) and 1(b).

(a) 𝑥1(𝑡)=𝑡4

(b) 𝑥2(𝑡)=𝑒2𝑡

Figure 1 

Analytic versus numerical approximation for a fixed 𝑛 (𝑛=2).

A different approach is to consider a fixed 𝑁 and vary the size of the expansion, that is, to consider different values for 𝑛. For the same functions 𝑥₁ and 𝑥₂, with 𝑁=6, the results are shown in Figures 2(a) and 2(b).

(a) 𝑥1(𝑡)=𝑡4

(b) 𝑥2(𝑡)=𝑒2𝑡

Figure 2 

Analytic versus numerical approximation for a fixed 𝑁 (𝑁=6).

4. Examples

We will see that applying the numerical method given in the previous section, we are able to solve fractional optimal control problems applying known techniques from the classical optimal control theory. First, consider the following optimal control problem: 𝐽(𝑥,𝑢)=¹₀(𝑡𝑢(𝑡)−(𝛼+2)𝑥(𝑡))²𝑑𝑡⟶min(41) subject to the control system 𝑥̇(𝑡)+^𝐶₀𝐷^𝛼_𝑡𝑥(𝑡)=𝑢(𝑡)+𝑡²(42) and the boundary conditions 𝑥(0)(1)=2=0,𝑥Γ().3+𝛼(43) The solution is given by (𝑥(𝑡),𝑢(𝑡))=2𝑡^𝛼+2Γ(),𝛼+32𝑡^𝛼+1Γ().𝛼+2(44) Using the necessary conditions given in Theorem 4, we arrive at 𝑥̇(𝑡)+^𝐶₀𝐷^𝛼_𝑡𝑥(𝑡)𝜆=−2𝑡²+𝛼+2𝑡𝑥(𝑡)+𝑡²,−𝜆̇(𝑡)+_𝑡𝐷^𝛼₁𝜆(𝑡)=𝛼+2𝑡𝜆(𝑡),𝑥(0)𝑥(1)=2=0,Γ(),3+𝛼(45) which is a fractional boundary value problem. We approximate this problem by approximation in (8) up to order 𝑁: 𝑥̇(𝑡)=𝛼+2𝑡−𝐴𝑡^−𝛼𝑥(𝑡)+^𝑁_𝑝=2𝐶_𝑝𝑡^{1−𝑝−𝛼}𝑉_𝑝(𝑡)−𝜆(𝑡)2𝑡²+𝑡²×11+𝐵𝑡^1−𝛼,𝑉̇_𝑝(𝑡)=1−𝑝𝑡^𝑝−2𝑥(𝑡)𝜆̇(𝑡)=𝐴(),𝑝=2,…,𝑁,1−𝑡^−𝛼−𝛼+2𝑡𝜆(𝑡)−^𝑁_𝑝=2𝐶_𝑝(1−𝑡)^{1−𝑝−𝛼}𝑊_𝑝(𝑡)×11+𝐵(1−𝑡)^1−𝛼,𝑊̇_𝑝(𝑡)=−1−𝑝(1−𝑡)^𝑝−2𝜆(𝑡),𝑝=2,…,𝑁,(46) subject to the boundary conditions 𝑥(0)(1)=2=0,𝑥Γ(),𝑉3+𝛼_𝑝(0)𝑊=0,𝑝=2,…,𝑁,_𝑝(1)=0,𝑝=2,…,𝑁.(47) The solutions are depicted in Figure 3, for 𝑁=2, 𝑁=3, and 𝛼=1/2, with the error being given by 𝐸=max_𝑖(|𝑥(𝑡_𝑖)−𝑥(𝑡_𝑖)|).

(a) 𝑥(𝑡);𝑁=2

(b) 𝑢(𝑡);𝑁=2

(c) 𝑥(𝑡);𝑁=3

(d) 𝑢(𝑡);𝑁=3

Figure 3 

Exact solution (solid line) for the fractional optimal control problem (41)–(43) with 𝛼=1/2 versus numerical solutions (dashed lines) obtained approximating the optimality conditions given by Theorem 4.

Another approach is to approximate the original problem by using approximation from (8) directly, getting 𝐽̃(𝑥,𝑢)=¹₀(𝑡𝑢−(𝛼+2)𝑥)²𝑑𝑡⟶min(48) subject to the control system 𝑥̇(𝑡)1+𝐵(𝛼,𝑁)𝑡^1−𝛼+𝐴(𝛼,𝑁)𝑡^−𝛼𝑥(𝑡)−^𝑁_𝑝=2𝐶𝛼,𝑝𝑡^{1−𝑝−𝛼}𝑉_𝑝(𝑡)=𝑢(𝑡)+𝑡²,𝑉̇_𝑝(𝑡)=1−𝑝𝑡^𝑝−2𝑥(𝑡)(49) and boundary conditions 𝑥(0)(1)=2=0,𝑥Γ(),𝑉3+𝛼_𝑝(0)=0,𝑝=2,3,…,𝑁.(50) The (classical) necessary optimality conditions become 𝑥̇(𝑡)=2𝜙₀(𝑡)𝜆₁(𝑡)+𝜙₁(𝑡)𝑥(𝑡)+^𝑁_𝑝=2𝜙_𝑝(𝑡)𝑉_𝑝(𝑡)+𝜙_𝑁+1(𝑡),𝑉̇_𝑝=1−𝑝𝑡^𝑝−2𝑥(𝑡)𝜆̇,𝑝=2,…,𝑁,₁=−𝜙₁(𝑡)𝜆₁(𝑡)+^𝑁_𝑝=2𝑝−1𝑡^𝑝−2𝜆_𝑝,𝜆̇_𝑝=−𝜙_𝑝(𝑡)𝜆₁(𝑡),𝑝=2,…,𝑁(51) subject to the boundary conditions 𝑥(0)𝑉=0,_𝑝(0)𝑥(1)=2=0,𝑝=2,…,𝑁;Γ(),𝜆3+𝛼_𝑝(1)=0,𝑝=2,…,𝑁.(52) The solutions are depicted in Figure 4 for 𝑁=2, 𝑁=3, and 𝛼=1/2.

(a) 𝑥(𝑡);𝑁=2

(b) 𝑢(𝑡);𝑁=2

(c) 𝑥(𝑡);𝑁=3

(d) 𝑢(𝑡);𝑁=3

Figure 4 

Exact solution (solid line) for the fractional optimal control problem (41)–(43) with 𝛼=1/2 versus numerical solutions (dashed lines) obtained approximating the original problem by a classical one.

For our next example, we consider the final time 𝑇 free and thus a variable in the problem. We wish to find an optimal triplet (𝑥(⋅),𝑢(⋅),𝑇) that minimizes 𝐽(𝑥,𝑢,𝑇)=^𝑇₀(𝑡𝑢−(𝛼+2)𝑥)²𝑑𝑡(53) subject to the control system 𝑥̇(𝑡)+^𝐶₀𝐷^𝛼_𝑡𝑥(𝑡)=𝑢(𝑡)+𝑡²(54) and boundary conditions 𝑥(0)=0,𝑥(𝑇)=1.(55) In this case, an exact solution to this problem is not known.

The fractional necessary optimality conditions, after approximating the fractional terms, result in 𝑥̇(𝑡)=𝛼+2𝑡−𝐴𝑡^−𝛼𝑥(𝑡)+^𝑁_𝑝=2𝐶_𝑝𝑡^{1−𝑝−𝛼}𝑉_𝑝(𝑡)−𝜆(𝑡)2𝑡²+𝑡²×11+𝐵𝑡^1−𝛼,𝑉̇_𝑝(𝑡)=1−𝑝𝑡^𝑝−2𝑥(𝑡)𝜆̇(𝑡)=𝐴(),𝑝=2,…,𝑁,1−𝑡^−𝛼−𝛼+2𝑡𝜆(𝑡)−^𝑁_𝑝=2𝐶_𝑝(1−𝑡)^{1−𝑝−𝛼}𝑊_𝑝(𝑡)×11+𝐵(1−𝑡)^1−𝛼,𝑊̇_𝑝(𝑡)=−1−𝑝(1−𝑡)^𝑝−2𝜆(𝑡),𝑝=2,…,𝑁(56) subject to the boundary conditions 𝑥(0)(𝑇)𝑉=0,𝑥=1,_𝑝(0)𝑊=0,𝑝=2,…,𝑁,_𝑝(𝑇)=0,𝑝=2,…,𝑁.(57) Another way is transforming the problem into an integer order optimal control problem with free final time. The necessary optimality conditions are 𝑥̇(𝑡)=2𝜙₀(𝑡)𝜆₁(𝑡)+𝜙₁(𝑡)𝑥(𝑡)+𝜙₂(𝑡)𝑉₂(𝑡)+𝜙₃(𝑡),𝑉̇₂(𝑡),𝜆̇=−𝑥₁=−𝜙₁(𝑡)𝜆₁(𝑡)(𝑡),𝜆̇+𝑥₂=−𝜙₂(𝑡)𝜆₁(𝑡).(58) The results obtained are shown in Figure 5.

(a) 𝑥(𝑡);𝑁=2

(b) 𝑢(𝑡);𝑁=2

Figure 5 

Numerical solutions to the free final time problem (53)–(55), using fractional necessary optimality conditions (dashed lines) and approximation of the problem to an integer order optimal control problem (dash-dotted lines).

Acknowledgments

This work was supported by FEDER funds through COMPETE—Operational Programme Factors of Competitiveness (“Programa Operacional Factores de Competitividade”), Portuguese funds through the Center for Research and Development in Mathematics and Applications (University of Aveiro), and the Portuguese Foundation for Science and Technology (“FCT—Fundação para a Ciência e a Tecnologia”), within project PEst-C/MAT/UI4106/2011 with Compete no. FCOMP-01-0124-FEDER-022690. Pooseh was also supported by the FCT Ph.D. fellowship SFRH/BD/33761/2009. This work is part of the first author’s Ph.D., which was carried out at the University of Aveiro under the Doctoral Program in Mathematics and Applications (PDMA) of Universities of Aveiro and Minho.