Fractional Order Optimal Control Problems with Free Terminal Time

We consider fractional order optimal control problems in which the dynamic control system involves integer and fractional order derivatives and the terminal time is free. Necessary conditions for a state/control/terminal-time triplet to be optimal are obtained. Situations with constraints present at the end time are also considered. Under appropriate assumptions, it is shown that the obtained necessary optimality conditions become sufficient. Numerical methods to solve the problems are presented, and some computational simulations are discussed in detail.


1.
Introduction. Fractional calculus generalizes the standard integral and differential calculus to arbitrary order. If the order of the fractional derivative operator is an integer m, then we recover an m-fold integral when m is negative, and the classical derivative of order m when m is positive. Let x : [a, b] → R be a function, and let α be a real number representing the order of the integral/derivative of x. The most common fractional operators considered in the literature are the Riemann-Liouville fractional integrals (RLFI) and derivatives (RLFD), and the Caputo fractional derivatives (CFD): where n stands for ⌊α⌋ + 1 (n is the smallest integer larger than α), and Γ is the gamma function, that is, If x and x (i) , i = 1, . . . , n − 1, vanish at t = a, then a D α t x(t) = C a D α t x(t), and if they vanish at t = b, then t D α b x(t) = C t D α b x(t). The reader interested in the theory and applications of fractional calculus is referred to the books [11,16]. For the state of the art on fractional order optimal control, see the recent book [15].
For numerical purposes, it is sometimes useful to approximate the fractional operators as sums of integer order (standard) derivatives. The most common procedure is to replace the Riemann-Liouville fractional derivative by the series (cf. [16]). This formula can be easily deduced by applying integration-by-parts to the integral 1 Γ(n − α) The main drawback is that to obtain good accuracy, we need higher-order derivatives of x, which restricts the class of admissible functions that one can consider. Recently, in [5], a new expansion formula was obtained, with the advantage that we only need the first derivative: where V p (t) is defined as the solution of the system V p (t) = (1 − p)(t − a) p−2 x(t), V p (a) = 0, p = 2, 3, . . ., and the coefficients A, B and C are given by the formulas For computational purposes, we truncate the sum and consider the finite expansion where A(α, N ) and B(α, N ) are now defined by We refer to [21] where expansion formulas with higher-order derivatives are obtained, and an error estimation is proven. In [20] analogous results are proven for the Riemann-Liouville fractional integral and in [19] for Hadamard fractional operators. For the right Riemann-Liouville fractional derivative, we have where W p is the solution of the differential equation To approximate the Caputo fractional derivatives, we may use the formulas that relate Caputo and Riemann-Liouville fractional operators, and then use (1) and (2). In this article we are interested in investigating such ideas in the context of fractional order optimal control. An important tool is the integration by parts formula for Caputo fractional derivatives, which is stated in the following theorem.
Theorem 1.1 (cf., e.g., [2]). Let α ∈ (0, 1), and x, y : [a, b] → R be two functions of class C 1 . Then the following integration by parts formula holds: The text is organized as follows. In Section 2 we formulate the optimal control problem under consideration and deduce necessary optimality conditions for it (Theorem 2.1). Another approach consists of using the approximation methods mentioned above, thereby converting the original problem into a classical optimal control problem that can be solved by standard computational techniques (Section 2.2). A generalization of the obtained results is considered in Section 3, where the lower bound of the cost integral is different from the lower bound of the fractional derivative (Theorem 3.1). Under some additional assumptions, the necessary optimality conditions are also sufficient (Theorem 4.1). We end with Section 5 containing numerical computations.
2. Necessary optimality conditions. Let α ∈ (0, 1), a ∈ R, L and f be two differentiable functions with domain [a, +∞) × R 2 , and φ : [a, +∞) × R → R a differentiable function. The fundamental problem is stated in the following way: subject to the control system and the initial boundary condition with (M, N ) = (0, 0) and x a a fixed real number. Our goal is to generalize previous works on fractional optimal control problems by considering the end time T free and the dynamic control system (4) involving integer and fractional order derivatives. For convenience, we consider the one-dimensional case. However, using similar techniques as the ones given here, the results can be easily extended to problems with multiple states and multiple controls. Later we consider the cases T and/or x(T ) fixed. Here, T is a variable number with a < T < ∞. Thus, we are interested not only on the optimal trajectory x and optimal control function u, but also on the corresponding time T for which the functional J attains its minimum value. We assume that the state variable x is differentiable and that the control u is piecewise continuous. When N = 0 we obtain a classical optimal control problem; the case M = 0 with fixed T has already been studied for different types of fractional order derivatives (see, e.g., [1,3,4,8,9,22,23]). In [10] a special type of the proposed problem is also studied for fixed T .
Remark 1. In this paper the terminal time T is a free decision variable and, a priori, no constraints are imposed. For future research, one may wish to consider a class of fractional optimal control problems in which the terminal time is governed by a stopping condition. Such problems were recently investigated, within the classical (integer order) framework, in [13,14].
2.1. Fractional necessary conditions. To deduce necessary optimality conditions that an optimal triplet (x, u, T ) must satisfy, we use a Lagrange multiplier to adjoin the dynamic constraint (4) to the performance functional (3). To start, we define the Hamiltonian function H by where λ is a Lagrange multiplier, so that we can rewrite the initial problem as minimizing Next, we consider variations of the form x + δx, u + δu, T + δT, λ + δλ with δx(a) = 0 by the imposed boundary condition (5). Using the well-known fact that the first variation of J must vanish when evaluated along a minimizer, we get with the partial derivatives of H evaluated at (t, x(t), u(t), λ(t)). Integration by parts gives the relations Thus, we deduce the following formula: Becauseδ x(T ) is arbitrary, in particular one can consider variation functions for whichδ x(T ) = 0. By Taylor's theorem, we arrive at the expression Since the variation functions were chosen arbitrarily, the following theorem is proven.
is a minimizer of (3) under the dynamic constraint (4) and the boundary condition (5), then there exists a function λ for which the triplet (x, u, λ) satisfies: • the Hamiltonian system for all t ∈ [a, T ]; • the stationary condition for all t ∈ [a, T ]; • and the transversality conditions where the Hamiltonian H is defined by (6).

Remark 2.
In standard optimal control, a free terminal time problem can be converted into a fixed final time problem by using the well-known transformation s = t/T (see Example 2). This transformation does not work in the fractional setting. Indeed, in standard optimal control, translating the problem from time t to a new time variable s is straightforward: the chain rule gives dx ds = dx dt dt ds . For Caputo or Riemann-Liouville fractional derivatives, the chain rule is not valid and such conversion is not possible.
Some interesting special cases are obtained when restrictions are imposed on the end time T or on x(T ).
1. If T is fixed and x(T ) is free, then Theorem 2.1 holds with the transversality conditions (9) replaced by 2. If x(T ) is fixed and T is free, then Theorem 2.1 holds with the transversality conditions (9) replaced by 3. If T and x(T ) are both fixed, then Theorem 2.1 holds with no transversality conditions. 4. If the terminal point x(T ) belongs to a fixed curve, i.e., x(T ) = γ(T ) for some differentiable curve γ, then Theorem 2.1 holds with the transversality conditions (9) replaced by 5. If T is fixed and x(T ) ≥ K for some fixed K ∈ R, then Theorem 2.1 holds with the transversality conditions (9) replaced by 6. If x(T ) is fixed and T ≤ K for some fixed K ∈ R, then Theorem 2.1 holds with the transversality conditions (9) replaced by Proof. The first three conditions are obvious. The fourth follows from To prove 5, observe that we have two possible cases. If x(T ) > K, then δx T may take negative and positive values, and so we get On the other hand, if x(T ) = K, then δx T ≥ 0 and so The proof of the last condition is similar.

2.2.
Approximated integer order necessary optimality conditions. Using approximation (1) up to order K, we can transform the original problem (3)-(5) into the following classical problem: where A = A(α, K), B = B(α, K) and C p = C(α, p) are the coefficients in the approximation (1). Now that we are dealing with an integer order problem, we can follow a classical procedure (see, e.g., [12]), by defining the Hamiltonian H by Let λ = (λ 1 , λ 2 , . . . , λ K ) and x = (x, V 2 , . . . , V K ). The necessary optimality conditions result in a two point boundary value problem. Assume that (T * , x * , u * ) is the optimal triplet. In addition to the boundary conditions (10), the transversality conditions imply where tr denotes the transpose. Because V p , p = 2, . . . , K, are auxiliary variables whose values V p (T ), at the final time T , are free, we have 3. A generalization. The aim is now to consider a generalization of the optimal control problem (3)-(5) studied in Section 2. Observe that the initial point t = a is in fact the initial point for two different operators: for the integral in (3) and, secondly, for the left Caputo fractional derivative given by the dynamic constraint (4). We now consider the case where the lower bound of the integral of J is greater than the lower bound of the fractional derivative. The problem is stated as follows: under the constraints where (M, N ) = (0, 0), x A is a fixed real, and a < A.

Remark 3.
We have chosen to consider the initial condition on the initial time A of the cost integral, but the case of initial condition x(a) instead of x(A) can be studied using similar arguments. Our choice seems the most natural: the interval of interest is [A, T ] but the fractional derivative is a nonlocal operator and has "memory" that goes to the past of the interval [A, T ] under consideration.

Remark 4.
In the theory of fractional differential equations, the initial condition is given at t = a. To the best of our knowledge there is no general theory about unicity of solutions for problems like (12), where the fractional derivative involves x(t) for a < t < A and the initial condition is given at t = A. Unicity of solution is, however, possible. Consider, for example, C 0 D α t x(t) = t 2 . Applying the fractional integral to both sides of equality we get x(t) = x(0) + 2t 2+α /Γ(3 + α) so, knowing a value for x(t), not necessarily at t = 0, one can determine x(0) and by doing so x(t). A different approach than the one considered here is to provide an initialization function for t ∈ [a, A]. This initial memory approach was studied for fractional continuous-time linear control systems in [17] and [18], respectively for Caputo and Riemann-Liouville derivatives.
The method to obtain the required necessary optimality conditions follows the same procedure as the one discussed before. The first variation gives where the Hamiltonian H is as in (6). Now, if we integrate by parts, we get Substituting these relations into the first variation of J, we conclude that Repeating the calculations as before, we prove the following optimality conditions. Theorem 3.1. If the triplet (x, u, T ) is an optimal solution to problem (11)-(12), then there exists a function λ for which the following conditions hold: with the Hamiltonian H given by (6).  Then (x, u) is an optimal solution to problem (3)-(5).

Numerical treatment and examples.
Here we apply the necessary conditions of Section 2 to solve some test problems. Solving an optimal control problem, analytically, is an optimistic goal and is impossible except for simple cases. Therefore, we apply numerical and computational methods to solve our problems. In each case we try to solve the problem either by applying fractional necessary conditions or by approximating the problem by a classical one and then solving the approximate problem.

Fixed final time.
We first solve a simple problem with fixed final time. In this case the exact solution, i.e., the optimal control and the corresponding optimal trajectory, is known, and hence we can compare it with the approximations obtained by our numerical method.
Let us apply the fractional necessary conditions to the above problem. The Hamiltonian is H = (tu − (α + 2)x) 2 + λu + λt 2 . The stationary condition (8) implies that for t = 0 Finally, (7) gives . At this point, we encounter a fractional boundary value problem that needs to be solved in order to reach the optimal solution. A handful of methods can be found in the literature to solve this problem. Nevertheless, we use approximations (1) and (2), up to order N , that have been introduced in [5] and used in [10,21]. With our choice of approximation, the fractional problem is transformed into a classical (integer order) boundary value problem: The solutions are depicted in Figure 1 for N = 2, N = 3 and α = 1/2. Since the exact solution for this problem is known, for each N we compute the approximation error by using the maximum norm. Assume that x(t i ) are the approximated values on the discrete time horizon a = t 0 , t 1 , . . . , t n . Then the error is given by Another approach is to approximate the original problem by using (1) for the fractional derivative. Following the procedure discussed in Section 2, the problem of Example 1 is approximated bỹ subject to the control system     (1) and (2) up to order N in the fractional necessary optimality conditions. and boundary conditions , V p (0) = 0, p = 2, 3, . . . , N.
The Hamiltonian system for this classical optimal control problem is Using the stationary condition ∂H ∂u = 0, we have for t = 0.

Finally, the Hamiltonian becomes
where and The Hamiltonian systemẋ = ∂H ∂λ ,λ = − ∂H ∂x , gives Example 2. Find an optimal triplet (x(·), u(·), T ) that minimizes subject to the control systeṁ and boundary conditions An exact solution to this problem is not known and we apply the two numerical procedures already used with respect to the fixed final time problem in Example 1.
We begin by using the fractional necessary optimality conditions that, after approximating the fractional terms, result in     The only difference here with respect to Example 1 is that there is an extra unknown, the terminal time T . The boundary condition for this new unknown is chosen appropriately from the transversality conditions discussed in Corollary 1, i.e., [H(t, x, u, λ) − λ(t) C a D α t x(t) +ẋ(t) t I 1−α T λ(t)] t=T = 0, where H is given as in (13). Since we require λ to be continuous, t I 1−α T λ(t)| t=T = 0 (cf. [16, pag. 46]) and so λ(T ) = 0. One possible way to proceed consists in translating the problem into the interval [0, 1] by the change of variable t = T s [6]. In this setting, either we add T to the problem as a new state variable with dynamicsṪ (s) = 0, or we treat it as a parameter. We use the latter, to get the   This parametric boundary value problem is solved for N = 2 and α = 0.5 with Matlab's bvp4c function. The result is shown in Figure 3 (dashed lines).
We also solve Example 2 with α = 1/2 by directly transforming it into an integer order optimal control problem with free final time. As is well known in the classical theory of optimal control, the Hamiltonian must vanish at the terminal point when the final time is free, i.e., one has H| t=T = 0 with H given by (14) [12]. For N = 2, the necessary optimality conditions give the following two point boundary value problem: where φ 0 (t) and φ 1 (t) are given by (15) and φ 2 (t) and φ 3 (t) by (16) with p = N = 2.
The trajectory x and corresponding u are shown in Figure 3 (dash-dotted lines).