Legendre spectral-collocation method for solving some types of fractional optimal control problems.

In this paper, the Legendre spectral-collocation method was applied to obtain approximate solutions for some types of fractional optimal control problems (FOCPs). The fractional derivative was described in the Caputo sense. Two different approaches were presented, in the first approach, necessary optimality conditions in terms of the associated Hamiltonian were approximated. In the second approach, the state equation was discretized first using the trapezoidal rule for the numerical integration followed by the Rayleigh-Ritz method to evaluate both the state and control variables. Illustrative examples were included to demonstrate the validity and applicability of the proposed techniques.

Introduction

Differential Equations (DEs) play a major role in mathematical modeling of real-life models in engineering, science and many other fields. Generally speaking the analytical methods are not suitable for large scale problems with complex solution regions. Numerical methods are commonly used to get an approximate solution for the DEs which are non-linear or the derivation of the analytical methods is difficult. Numerical methods for DEs have been explored rapidly with the development of digital computers. Optimal control deals with the problem of finding a control law for a given dynamical system. An optimal control problem is a set of DEs describing the paths of the control variables that minimize a function of state and control variables. A necessary condition for an optimal control problem can be derived using Pontryagin’s maximum principle and a sufficient condition can be obtained using Hamilton–Jacobi–Bellman equation.

Fractional order DEs have gained considerable importance due to their application in various sciences, such as physics, mechanics, chemistry, and engineering. Fractional order models are more appropriate than conventional integer order to describe physical systems [1-4]. For example, it has been illustrated that the so-called fractional Cable equation, which is similar to the traditional Cable equation except that the order of derivative with respect to the space and/or time is fractional, can be more adequately modeled by fractional order models than integer order models [5].

In the recent years, the dynamic behaviors of fractional-order differential systems have received increasing attention. FOCP refers to the minimization of an objective functional subject to dynamic constraints, on state and control variables, which have fractional order models. Some numerical methods for solving some types of FOCPs were recorded [6-10] and the references cited therein.

This paper is a continuation of the authors work in this area of research [9,10]. The main aim of this work was to use the advantage of the Legender spectral-collocation method to study FOCPs, two efficient numerical methods for solving some types of FOCPs are presented where fractional derivatives are introduced in the Caputo sense. These numerical methods depend upon the spectral method where the Legendre polynomials are used to approximate the unknown functions. Legendre polynomials are well known family of orthogonal polynomials on the interval that have many applications [11]. They are widely used because of their good properties in the approximation of functions.

The structure of this paper was arranged in the following way: In Section ‘Preliminaries and notations’, preliminaries, notations and properties of the shifted Legendre polynomials were introduced. In Section ‘Necessary optimality conditions’, necessary optimality conditions of the FOCP model were given. In Section ‘Numerical approximation’, the basic formulation of the proposed approximate formulas of the fractional derivatives was obtained. In Section ‘Error estimates’, error estimates for the approximated fractional derivatives were given. In Section ‘Numerical results’, illustrative examples were included to demonstrate the validity and applicability of the proposed technique. Finally, in Section ‘Conclusions’, this paper ends with a brief conclusion and some remarks.

Preliminaries and notations

Fractional derivatives and integrals

Let be a function, a real number, and , where denotes the smallest integer greater than or equal to . The left (left RLFI) and right (right RLFI) Riemann–Liouville fractional integrals are defined, respectively, by:The left (left RLFD) and right (right RLFD) Riemann–Liouville fractional derivatives are defined, respectively, by:

The left (left CFD) and right (right CFD) Caputo fractional derivatives are defined respectively, by:In the following some basic properties are presented:

The relation between right RLFD and right CFD [12]:

where . 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 its properties see [13,14].

The shifted Legendre polynomials

The well known Legendre polynomials are defined on the interval and can be determined with the aid of the following recurrence formula [15]:The analytic form of the Legendre polynomials of degree n is given bywhere denotes the biggest integer less than or equal to n. Moreover, we have [16]:andIn order to use these polynomials on the interval we use the so-called shifted Legendre polynomials by introducing the change of variable . The shifted Legendre polynomials are defined as follows:The analytic form of the shifted Legendre polynomials of degree n is given by:Note that from Eq. (9), we can see that .

The function which belongs to the space of square integrable in , may be expressed in terms of shifted Legendre polynomials aswhere the coefficients are given by:

Necessary optimality conditions

Let and let be two differentiable functions.

Consider the following FOCP [8]:subject to the dynamic system:where the boundary conditions are as follows:where are fixed real numbers.

[8] If is a minimizer of (11)–(13), then there exists an adjoint state for which the triple satisfies the optimality conditions

for all ,

and the transversality condition:where the Hamiltonian H is defined byIf is fixed, there is no transversality condition.

Under some additional assumptions on the objective functional L and the right-hand side f, e.g., convexity of L and linearity of f in x and u, the optimality conditions (14)–(16) are also sufficient.

Numerical approximation

In this section, numerical approximations for the left CFD and the right RLFD using Legendre polynomials are presented.

Let be a function defined on the interval , and N be positive integer. Denote bywhere is an approximation of . If is the interpolation of on the Legendre–Gauss–Lobatto points , then can be determined bywhere for , and are the corresponding quadrature weights [17,18].

In the following, approximation of the fractional derivative is given.

[9] let be approximated by shifted Legendre polynomials as (18) and (19) and also , then

where is given by:

Approximation of right RLFD

Let be a sufficiently smooth function in and be defined as follows:from (2) and (3), we have:let be approximated by shifted Legendre polynomials as (18) and (19)

Then we claim:

Let be a polynomial of degree N given by (18). Then there exists a polynomial of degree such that

Let be expanded in Taylor series at as follows:where .

Then,Then,We have (24) if we choosewith an arbitrary constant . □

From (24) we have:and can be approximated as follows,

Now, we express in (25) by a sum of the Legendre polynomials and show the recurrence relation satisfied by the Legendre coefficients. Differentiating both sides of (24) with respect to x yieldsThen,To evaluate in (25) we expand in terms of the shifted Legendre polynomialsIntegrating both sides of (28) giveswhere . On the other hand, we haveThen, by using the relation and Eq. (28), we have:where . Let

By inserting and given by (29) and (30), respectively, into (27), and from (31), we have:The Legendre coefficients of given by (31) can be evaluated by integrating (31) and comparing it with (18) and (19)with starting values , where are the Legendre coefficients of .

Error estimates

In the following, we give an upper bound for the coefficients of Legendre expansion of a function f on .

If are absolutely continuous on and if for some , then for each ,

We have:

Using the substitution , we have:

Integrating by parts, using Eq. (8),

Again, integrating by parts,

For , to keep the formula simple, we do not keep track of these different denominators but weaken the inequality slightly by replacing them with ,since and .

Further integrations by parts, The result is Eq. (34).  □

Suppose that f satisfies hypotheses of Lemma 4. Let be the truncated Legendre expansion of f. Then for and ,

We have:since Eq. (7). Then, from Lemma 4,

Now, in order to estimate the error of the approximated fractional derivatives, we have to estimate the error of the first derivative of the LGL interpolation as the following.

Suppose that f satisfies hypotheses of Lemma 5. Let be LGL interpolation of f. Assume that and . We have for :

Markov’s inequality asserts thatfor all polynomials of degree at most n with real coefficients [19], so

ButSince in [16]:, where is the best approximation of f and is the Lebesgue constant for which the following estimate holds, on [20], and from Eq. (7) and Lemma 4, we haveIt is well known that the truncated Chebyshev expansion is very close to the best polynomial approximation [21]. Therefore, from [22] (we reformulate the Chebyshev error bound on ,Hence,

Numerical results

In this section, we develop two algorithms (Algorithms 1 and 2) for the numerical solution of FOCPs and apply them to two illustrative examples. For the first Algorithm, 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. We solve the fractional boundary value problems by the spectral method. The second Algorithm relies on the strategy “discretize first, then optimize”. The Rayleigh–Ritz method provides the optimality conditions in the discrete regime.

We consider the following FOCP from [8,10]:subject to the dynamical systemand the boundary conditionsThe exact solution is given by

The first algorithm for the solution of (36)–(38) follows the “optimize first, then discretize” approach. It is based on the necessary optimality conditions from Theorem 1 and implements the following steps:

Step 1: Compute the Hamiltonian

Step 2: Derive the necessary optimality conditions from Theorem 1:Use (43) in (41) and (42) to obtain

Step 3: By using Legendre expansion, get an approximate solution of the coupled system (44) and (45) under the boundary conditions (38):

Step 3a: In order to solve (44) by the Legendre expansion method, use (18) and (19) to approximate . A collocation scheme is defined by substituting (18), (19), (20) and (26) into (44) and evaluating the results at the shifted Legendre–Gauss–Lobatto nodes . This gives:, where is defined in (21). The system (46) represents algebraic equations which can be solved for the unknown coefficients .

Consequently, it remains to compute the two unknowns . This can be done by using any two points which differ from the Legendre–Gauss–Lobatto nodes and satisfy (44). We end up with two equations in two unknowns:

Step 3b: In order to solve (45) by the Legendre expansion method, we use (18) and (19) to approximate the state x. A collocation scheme is defined by substituting (18)–(20) and then computed into (45) and evaluating the results at the shifted Legendre–Gauss–Lobatto nodes . This results in system of algebraic equations which can be solved for the unknown coefficients . By using the boundary conditions, we have and . Figs. 1a,1b,1c and 1d display the exact and approximate state x and the exact and approximate control u for and , 3. Table 1 contains the maximum errors in the state x and in the control u for and .

Table 1

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

	N=2	N=3	N=5
Max. error in x	3.1055E-2	4.0702E-3	3.5526E-4
Max. error in u	2.0410E-1	4.5860E-2	9.1353E-3

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

Step 1: Substitute (37) into (36) to obtain

Step 2: Approximate x using the Legendre expansion (18) and (19) and approximate the Caputo fractional derivative and using (20) on the Legendre–Gauss–Lobatto nodes. Then, (47) takes the formwhere is defined as in (21).

Step 3: DefineUsing the composite trapezoidal integration technique,

Step 4: The extremal values of functionals of the general form (6.1), according to Rayleigh–Ritz method giveso, after using the boundary conditions, we obtain a system of algebraic equations.

Step 5: Solve the algebraic system by using the Newton–Raphson method to obtain and using the boundary conditions to get , then the function which extremes FOCPs has the following form:Figs. 1e,1f,1g and 1h display the exact and approximate state x and the exact and approximate control u for , and . Table 2 contains the maximum errors in the state x and in the control u for and .

Table 2

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

	N=2	N=3	N=5
Max. error in x	2.7313E-2	2.2570E-3	1.6006E-4
Max. error in u	2.5699E-1	4.4538E-2	8.2254E-3

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

We consider the following linear-quadratic optimal control problem [10]:subject to the dynamical systemand the boundary conditions

The exact solution is given by

We note that for Example 2 the optimality conditions stated in Theorem 1 are also sufficient (cf. Remark 1).

Table 3 contains a comparison between the maximum error in the state x and in the control u for Algorithms 1 and 2.

Table 3

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

N=3	Alg. 1	Alg. 2
Max. error in x	8.8025E-3	5.1966E-3
Max. error in u	8.8025E-3	4.3260E-2
N=5	Alg. 1	Alg. 2
Max. error in x	1.0903E-4	4.5321E-5
Max. error in u	1.0903E-4	6.3134E-4

The next two examples are modifications of the problems presented in [23,24].

Consider the following time invariant problem:subject to the dynamical systemand the boundary conditionswhereFor this problem we have the exact solution in the case of as follows [24]:

Figs. 2a and 2b display Algorithm 1 approximate solutions of and for and , 0.9, 0.99, and exact solution for .

Figs. 2c and 2d display Algorithm 1 approximate solutions of and for , 5 and , and exact solution for .

Figs. 2e and 2f display Algorithm 2 approximate solutions of and for and , 0.9, 0.99 and exact solution for .

Figs. 2g and 2h display Algorithm 2 approximate solutions of and for , 5 and and exact solution for .

Figs. 2b, 2d, 2f and 2h illustrate that the approximate control converges better to the exact solution in Algorithm 1 than Algorithm 2.

Table 4 contains a comparison between approximate J in Algorithms 1 and 2 for “ with different values of ” and “ with ” where the exact is “ for ”.

Table 4

Approximate J for Algorithms 1 and 2.

N=3	J, Alg. 1	J, Alg. 2
α=0.8	0.193035	0.185312
α=0.9	0.193929	0.196629
α=0.99	0.195687	0.212169
N=5	J, Alg. 1	J, Alg. 2
α=0.9	0.187676	0.19636

Consider the following time variant problem:subject to the dynamical system,and the initial condition,Algorithm 1 has a modification to step 3a and step 3b where we have and and we use any two point which differ from LGL nodes and satisfy the necessary equation like (44) or (44) to determine and . Also in Algorithm 2 , there is a modification to step 5 where we solve the non-linear algebraic system of equations to obtain and use the initial condition to get .

Figs. 3a and 3b display Algorithm 1 approximate solutions of and for and .

Figs. 3c and 3d display Algorithm 2 approximate solutions of and for and .

Table 5 contains a comparison between approximate J in Algorithms 1 and 2 for different values of and .

Table 5

Approximate J for Algorithms 1 and 2.

N=3	J, Alg. 1	J, Alg. 2
α=0.8	0.488123	0.481819
α=0.9	0.487306	0.487719
α=0.99	0.484141	0.497106

Conclusions

In this work, Legendre spectral-collocation method is used to study some types of fractional optimal control problems. Two efficient algorithms for the numerical solution of a wide class of fractional optimal control problems are presented. In the first algorithm we 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 drive an approximation of right Riemann–Liouville fractional derivatives and solve these fractional boundary value problems using the spectral method. In the second algorithm, the state equation is adjoined to the objective functional which discretized and then the composite trapezoidal integration technique and the Rayleigh–Ritz method are used to evaluate both the state and control variables. In both algorithms, the solution is approximated by N-term truncated Legendre series. Numerical results show that the two algorithms converge as the number of terms increase. For the first example, it is noted that Algorithm 2 is more accurate than Algorithm 1 but in the second one Algorithm 1 is better in finding the control variable. Also Examples 3 and 4 show that Algorithm 1 is preferable than Algorithm 2. In general, the two algorithms are efficient and give the optimum solution.

Conflict of interest

The authors have declared no conflict of interest

Compliance with Ethics Requirements

This article does not contain any studies with human or animal subjects.