# A Numerical Scheme to Solve Fractional Optimal Control Problems

Authors: Shakoor Pooseh; Ricardo Almeida; Delfim F. M. Torres
Journal: Conference Papers in Mathematics
 (2013)

---

## 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.


---

## Body

## 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)dnxmdxn=m(m-1)⋯(m-n+1)xm-n=Γ(m+1)Γ(m-n+1)xm-n,
Euler suggested to use this relationship also for real values of n. Taking m=1 and n=1/2, he obtained
(2)d1/2xdx1/2=4xπ.
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 n-fold integral,
(3)∫at‍dτ1∫aτ1‍dτ2⋯∫aτn-1‍x(τn)dτn=1(n-1)!∫at‍(t-τ)n-1x(τ)dτ,
Riemann defined fractional integration as
(4)aItαx(t)=1Γ(α)∫at‍(t-τ)α-1x(τ)dτ.
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.
Letx:[a,b]→ℝ be a function, α>0 a real, and n=[α]+1, where [·] denotes the integer part function. The left and right Riemann-Liouville fractional integrals are defined, respectively, by
(leftRLFI)aItαx(t)=1Γ(α)∫at‍(t-τ)α-1x(τ)dτ,(rightRLFI)tIbαx(t)=1Γ(α)∫tb‍(τ-t)α-1x(τ)dτ.
The left and right Riemann-Liouville fractional derivatives are defined, respectively, by
(leftRLFD)aDtαx(t)=1Γ(n-α)dndtn∫at‍(t-τ)n-α-1x(τ)dτ,(rightRLFD)tDbαx(t)=(-1)nΓ(n-α)dndxn∫tb‍(τ-t)n-α-1x(τ)dτ.
The left and right Caputo fractional derivatives are defined, respectively, by
(leftCFD)aCDtαx(t)=1Γ(n-α)∫at‍(t-τ)n-α-1x(n)(τ)dτ,(rightCFD)tCDbαx(t)=(-1)nΓ(n-α)∫tb‍(τ-t)n-α-1x(n)(τ)dτ.We remark that ifα=n in Definition 1, then we have the usual operators:
(5)Iatnx(t)=∫at‍dτ1∫aτ1‍dτ2⋯∫aτn-1‍x(τn)dτn,tIbnx(t)=∫tb‍dτ1∫τ1b‍dτ2⋯∫τn-1b‍x(τn)dτn,aDtnx(t)=x(n)(t),tDbnx(t)=(-1)nx(n)(t),aCDtnx(t)=x(n)(t),tCDbαx(t)=(-1)nx(n)(t).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)
aCDtαx(t)=aDtαx(t)-∑k=0n-1‍(x(k)(a)/Γ(k-α+1))(t-a)k-α,(2)
aItαaItβx(t)=aItα+βx(t),(3)
aCDtαaItαx(t)=x(t),(4)
aItαaCDtαx(t)=x(t)-∑k=0n-1‍(x(k)(a)/k!)(t-a)k,(5)
∫ab‍y(t)·aCDtαx(t)dt=∫ab‍x(t)·tDbαy(t)dt+∑j=0n-1‍[tDbα+j-ny(t)·tDbn-1-jx(t)]ab.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(6)aDtαx(t)=∑n=0∞‍(αn)(t-a)n-αΓ(n+1-α)x(n)(t),
where
(7)(αn)=(-1)n-1αΓ(n-α)Γ(1-α)Γ(n+1).
Although very simple to use, it is easy to conclude that in order to have a small error when we approximateaDtαx by a finite sum up to order N, we need to consider a large value for N; that is, we need to consider the set of admissible functions to be CN[a,b] 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:
(8)aDtαx(t)=A(α)(t-a)-αx(t)+B(α)(t-a)1-αx˙(t)-∑p=2∞‍C(α,p)(t-a)1-p-αVp(t),
where Vp(t) is the solution of the system
(9)V˙p(t)=(1-p)(t-a)p-2x(t),Vp(a)=0,
for p=2,3,…, and A, B, and C are given by
(10)A(α)=1Γ(1-α)[1+∑p=2∞‍Γ(p-1+α)Γ(α)(p-1)!],B(α)=1Γ(2-α)[1+∑p=1∞‍Γ(p-1+α)Γ(α-1)p!],C(α,p)=1Γ(2-α)Γ(α-1)Γ(p-1+α)(p-1)!.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 L,f:[a,+∞[×ℝ2→ℝ be two differentiable functions and ϕ:[a,+∞[×ℝ→ℝ a differentiable function. The fundamental problem, as studied in [8], is the following:
(11)minimizeJ(x,u,T)=∫aT‍L(t,x(t),u(t))dt+ϕ(T,x(T))
subject to the dynamic control system
(12)Mx˙(t)+NaCDtαx(t)=f(t,x(t),u(t))
and the initial condition
(13)x(a)=xa,
with (M,N)≠(0,0) and xa being a fixed real number. Thus, we are not only interested in finding the optimal state function x and the optimal control u, but also the optimal time T.Theorem 3.
If(x,u,T) is a minimizer of (11) under the dynamic constraint (12) and the boundary condition (13), then there exists a function λ for which the triplet (x,u,λ) satisfies (i)
the Hamiltonian system(14)Mλ˙(t)-NtDTαλ(t)=-∂H∂x(t,x(t),u(t),λ(t)),Mx˙(t)+NaCDtαx(t)=∂H∂λ(t,x(t),u(t),λ(t))
for all t∈[a,T];(ii)
the stationary condition(15)∂H∂u(t,x(t),u(t),λ(t))=0,∀t∈[a,T];(iii)
the transversality conditions(16)[H(t,x(t),u(t),λ(t))-Nλ(t)aCDtαx(t)+Nx˙(t)tIT1-αλ(t)+∂ϕ∂t(t,x(t))]t=T=0,[Mλ(t)+NtIT1-αλ(t)-∂ϕ∂x(t,x(t))]t=T=0;where the Hamiltonian
H
is defined by
(17)
H
(
t
,
x
,
u
,
λ
)
=
L
(
t
,
x
,
u
)
+
λ
f
(
t
,
x
,
u
)
.This theorem states the general condition that the optimal solution(x,u,T) must fulfill. Next, depending on extra conditions imposed over the final time T or in x(T), new transversality conditions are obtained.Theorem 4.
Let(x,u) be a minimizer of (11) under the dynamic constraint (12) and the boundary condition (13). (i)
IfT is fixed and x(T) is free, then Theorem 3 holds with the transversality conditions (16) replaced by
(18)[Mλ(t)+NtIT1-αλ(t)-∂ϕ∂x(t,x(t))]t=T=0.(ii)
Ifx(T) is fixed and  T is free, then Theorem 3 holds with the transversality conditions (16) replaced by
(19)[H(t,x(t),u(t),λ(t))-Nλ(t)aCDtαx(t)+Nx˙(t)tIT1-αλ(t)+∂ϕ∂t(t,x(t))]t=T=0.(iii)
IfT and x(T) are fixed, then Theorem 3 holds with no transversality conditions.(iv)
If the terminal pointx(T) belongs to a fixed curve, that is, x(T)=γ(T) for some differentiable curve γ, then Theorem 3 holds with the transversality conditions (16) replaced by
(20)[∂ϕ∂xH(t,x(t),u(t),λ(t))-Nλ(t)aCDtαx(t)+Nx˙(t)tIT1-αλ(t)+∂ϕ∂t(t,x(t))-γ˙(t)(Mλ(t)+NtIT1-αλ(t)-∂ϕ∂x(t,x(t)))]t=T=0.(v)
IfT is fixed and x(T)≥K for some fixed K∈ℝ, then Theorem 3 holds with the transversality conditions (16) replaced by
(21)[Mλ(t)+NtIT1-αλ(t)-∂ϕ∂x(t,x(t))]t=T≤0,(x(T)-K)[Mλ(t)+NtIT1-αλ(t)-∂ϕ∂x(t,x(t))]t=T=0.Numerically, by using approximation (8) up to order K, we can transform the original problem into the following classical optimal control problem:
(22)minimizeJ~(x,u,T)=∫aT‍L(t,x(t),u(t))dt+ϕ(T,x(T))
subject to the dynamic constraints
(23)x˙(t)=f(t,x(t),u(t))-NA(t-a)-αx(t)M+NB(t-a)1-α+∑p=2KNCp(t-a)1-p-αVp(t)M+NB(t-a)1-α,V˙p(t)=(1-p)(t-a)p-2x(t),p=2,…,K
and the initial conditions
(24)x(a)=xa,Vp(a)=0,p=2,…,K.Theorem3 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 t=A, but similar conditions could be obtained if we considered conditions at t=a instead. The problem is formulated as follows. Let α∈(0,1), and let L,f:[a,+∞[×ℝ2→ℝ be two differentiable functions, ϕ:[a,+∞[×ℝ→ℝ a differentiable function, and A>a a real. We wish to
(25)minimizeJ(x,u,T)=∫AT‍L(t,x(t),u(t))dt+ϕ(T,x(T))
subject to
(26)Mx˙(t)+NaCDtαx(t)=f(t,x(t),u(t)),(27)x(A)=xA,
with (M,N)≠(0,0) and xA being a fixed real number.Theorem 5.
If(x,u,T) is a minimizer of (25) under the dynamic constraint (26) and the boundary condition (27), then there exists a function λ for which the triplet (x,u,λ) satisfies(i)
the Hamiltonian system(28)Mλ˙(t)-NtDTαλ(t)=-∂H∂x(t,x(t),u(t),λ(t)),Mx˙(t)+NaCDtαx(t)=∂H∂λ(t,x(t),u(t),λ(t))
for all t∈[A,T] and
(29)tDTαλ(t)-tDAαλ(t)=0
for all t∈[a,A];(ii)
the stationary condition(30)∂H∂u(t,x(t),u(t),λ(t))=0,∀t∈[A,T];(iii)
the transversality conditions(31)[H(t,x(t),u(t),λ(t))-Nλ(t)aCDtαx(t)+Nx˙(t)tIT1-αλ(t)+∂ϕ∂t(t,x(t))]t=T=0,[Mλ(t)+NtIT1-αλ(t)-∂ϕ∂x(t,x(t))]t=T=0,[tIT1-αλ(t)-tIA1-αλ(t)]t=a=0;where the Hamiltonian
H
is defined by
(32)
H
(
t
,
x
,
u
,
λ
)
=
L
(
t
,
x
,
u
)
+
λ
f
(
t
,
x
,
u
)
.We remark that whenA=a, Theorem 5 reduces to Theorem 3.Under some additional conditions, namely, convexity conditions overL, f, and ϕ, Theorem 3 provides also sufficient conditions to ensure optimal solutions. The result is given in the next theorem.Theorem 6.
Let(x¯,u¯,λ¯) be a triplet satisfying the necessary conditions of Theorem 3. Moreover, assume that  (1)
L and f are convex on x and u, and ϕ is convex in x;(2)
T is fixed;(3)
λ¯(t)≥0 for all t∈[a,T] or f is linear in x and u. Then (x¯,u¯) 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.Letx∈C2[a,b]. Using integration by parts two times, we deduce that
(33)aDtαx(t)=x(a)Γ(1-α)(t-a)-α+x˙(a)Γ(2-α)(t-a)1-α+(t-a)1-αΓ(2-α)∫at‍(1-τ-at-a)1-αx¨(τ)dτ.
By the binomial formula, we can rewrite the fractional derivative as
(34)aDtαx(t)=x(a)Γ(1-α)(t-a)-α+(t-a)1-αΓ(2-α)x˙(a)+(t-a)1-αΓ(2-α)∫at‍(∑p=0∞‍Γ(p-1+α)Γ(α-1)p!(τ-at-a)p)×x¨(τ)dτ.
Further integration by parts gives
(35)aDtαx(t)=A(α)(t-a)-αx(t)+B(α)(t-a)1-αx˙(t)-∑p=2∞‍C(α,p)(t-a)1-p-αVp(t),
where Vp(t), A(α), B(α), and C(α,p) are given in (9)-(10). Following similar calculations, we can deduce the next theorem.Theorem 7.
Fixn∈ℕ and let x∈Cn[a,b]. Then,
(36)aDtαx(t)=1Γ(1-α)(t-a)-αx(t)+∑i=1n-1‍A(α,i)(t-a)i-αx(i)(t)+∑p=n∞‍[-Γ(p-n+1+α)Γ(-α)Γ(1+α)(p-n+1)!(t-a)-αx(t)+B(α,p)(t-a)n-1-p-αVp(t)-Γ(p-n+1+α)Γ(-α)Γ(1+α)(p-n+1)\!],
where
(37)A(α,i)=1Γ(i+1-α)[1+∑p=n-i∞‍Γ(p-n+1+α)Γ(α-i)(p-n+i+1)!],B(α,p)=Γ(p-n+1+α)Γ(-α)Γ(1+α)(p-n+1)!,Vp(t)=(p-n+1)∫at‍(τ-a)p-nx(τ)dτ.The idea is to replace the fractional derivative with such expansions and to consider finite sums only. When we use the approximation(38)aDtαx(t)≈∑i=0n-1‍A(α,i,N)(t-a)i-αx(i)(t)+∑p=nN‍B(α,p)(t-a)n-1-p-αVp(t),
the error is bounded by
(39)|Etr(t)|≤Lne(n-1-α)2+n-1-αΓ(n-α)(n-1-α)Nn-1-α(t-a)n-α,
where
(40)Ln=maxτ∈[a,t]|x(n)(τ)|.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, considerx1(t)=t4 and x2(t)=e2t and expansions with n=2 and different values for N. The result is exemplified in Figures 1(a) and 1(b).Analytic versus numerical approximation for a fixedn (n=2).
(a)
x
1
(
t
)
=
t
4
(b)
x
2
(
t
)
=
e
2
tA different approach is to consider a fixedN and vary the size of the expansion, that is, to consider different values for n. For the same functions x1 and x2, with N=6, the results are shown in Figures 2(a) and 2(b).Analytic versus numerical approximation for a fixedN (N=6).
(a)
x
1
(
t
)
=
t
4
(b)
x
2
(
t
)
=
e
2
t


## 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:(41)J(x,u)=∫01‍(tu(t)-(α+2)x(t))2dt⟶min
subject to the control system
(42)x˙(t)+0CDtαx(t)=u(t)+t2
and the boundary conditions
(43)x(0)=0,x(1)=2Γ(3+α).
The solution is given by
(44)(x¯(t),u¯(t))=(2tα+2Γ(α+3),2tα+1Γ(α+2)).
Using the necessary conditions given in Theorem 4, we arrive at
(45)x˙(t)+0CDtαx(t)=-λ2t2+α+2tx(t)+t2,-λ˙(t)+tD1αλ(t)=α+2tλ(t),x(0)=0,x(1)=2Γ(3+α),
which is a fractional boundary value problem. We approximate this problem by approximation in (8) up to order N:
(46)x˙(t)=[(α+2t-At-α)x(t)+∑p=2N‍Cpt1-p-αVp(t)-λ(t)2t2+t2]×11+Bt1-α,V˙p(t)=(1-p)tp-2x(t),p=2,…,N,λ˙(t)=[(A(1-t)-α-α+2t)λ(t)-∑p=2N‍Cp(1-t)1-p-αWp(t)]×11+B(1-t)1-α,W˙p(t)=-(1-p)(1-t)p-2λ(t),p=2,…,N,
subject to the boundary conditions
(47)x(0)=0,x(1)=2Γ(3+α),Vp(0)=0,p=2,…,N,Wp(1)=0,p=2,…,N.
The solutions are depicted in Figure 3, for N=2, N=3, and α=1/2, with the error being given by E=maxi(|x(ti)-x¯(ti)|).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.
(a)
x
(
t
)
;
N
=
2
(b)
u
(
t
)
;
N
=
2
(c)
x
(
t
)
;
N
=
3
(d)
u
(
t
)
;
N
=
3Another approach is to approximate the original problem by using approximation from (8) directly, getting
(48)J~(x,u)=∫01‍(tu-(α+2)x)2dt⟶min
subject to the control system
(49)x˙(t)[1+B(α,N)t1-α]+A(α,N)t-αx(t)-∑p=2N‍C(α,p)t1-p-αVp(t)=u(t)+t2,V˙p(t)=(1-p)tp-2x(t)
and boundary conditions
(50)x(0)=0,x(1)=2Γ(3+α),Vp(0)=0,p=2,3,…,N.
The (classical) necessary optimality conditions become
(51)x˙(t)=2ϕ0(t)λ1(t)+ϕ1(t)x(t)+∑p=2N‍ϕp(t)Vp(t)+ϕN+1(t),V˙p=(1-p)tp-2x(t),p=2,…,N,λ˙1=-ϕ1(t)λ1(t)+∑p=2N‍(p-1)tp-2λp,λ˙p=-ϕp(t)λ1(t),p=2,…,N
subject to the boundary conditions
(52)x(0)=0,Vp(0)=0,p=2,…,N;x(1)=2Γ(3+α),λp(1)=0,p=2,…,N.
The solutions are depicted in Figure 4 for N=2, N=3, and α=1/2.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.
(a)
x
(
t
)
;
N
=
2
(b)
u
(
t
)
;
N
=
2
(c)
x
(
t
)
;
N
=
3
(d)
u
(
t
)
;
N
=
3For our next example, we consider the final timeT free and thus a variable in the problem. We wish to find an optimal triplet (x(·),u(·),T) that minimizes
(53)J(x,u,T)=∫0T‍(tu-(α+2)x)2dt
subject to the control system
(54)x˙(t)+0CDtαx(t)=u(t)+t2
and boundary conditions
(55)x(0)=0,x(T)=1.
In this case, an exact solution to this problem is not known.The fractional necessary optimality conditions, after approximating the fractional terms, result in(56)x˙(t)=[(α+2t-At-α)x(t)+∑p=2N‍Cpt1-p-αVp(t)-λ(t)2t2+t2]×11+Bt1-α,V˙p(t)=(1-p)tp-2x(t),p=2,…,N,λ˙(t)=[(A(1-t)-α-α+2t)λ(t)-∑p=2N‍Cp(1-t)1-p-αWp(t)]×11+B(1-t)1-α,W˙p(t)=-(1-p)(1-t)p-2λ(t),p=2,…,N
subject to the boundary conditions
(57)x(0)=0,x(T)=1,Vp(0)=0,p=2,…,N,Wp(T)=0,p=2,…,N.
Another way is transforming the problem into an integer order optimal control problem with free final time. The necessary optimality conditions are
(58)x˙(t)=2ϕ0(t)λ1(t)+ϕ1(t)x(t)+ϕ2(t)V2(t)+ϕ3(t),V˙2=-x(t),λ˙1=-ϕ1(t)λ1(t)+x(t),λ˙2=-ϕ2(t)λ1(t).
The results obtained are shown in 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).
(a)
x
(
t
)
;
N
=
2
(b)
u
(
t
)
;
N
=
2



---

*Source: 165298-2013-07-17.xml*