A numerical approach for solving fractional optimal control problems using modified hat functions

TL;DR

This paper presents a novel numerical method using modified hat functions to solve fractional optimal control problems, simplifying computations and avoiding numerical integration, with proven error bounds and applicability to constrained problems.

Contribution

The paper introduces a new numerical approach based on modified hat functions for fractional optimal control problems, including error analysis and extension to inequality constraints.

Findings

Effective in solving fractional optimal control problems

Reduces to solving nonlinear algebraic equations

Demonstrates high accuracy and ease of implementation

Abstract

We introduce a numerical method, based on modified hat functions, for solving a class of fractional optimal control problems. In our scheme, the control and the fractional derivative of the state function are considered as linear combinations of the modified hat functions. The fractional derivative is considered in the Caputo sense while the Riemann-Liouville integral operator is used to give approximations for the state function and some of its derivatives. To this aim, we use the fractional order integration operational matrix of the modified hat functions and some properties of the Caputo derivative and Riemann-Liouville integral operators. Using results of the considered basis functions, solving the fractional optimal control problem is reduced to the solution of a system of nonlinear algebraic equations. An error bound is proved for the approximate optimal value of the performance index obtained by the proposed method. The method is then generalized for solving a class of fractional optimal control problems with inequality constraints. The most important advantages of our method are easy implementation, simple operations, and elimination of numerical integration. Some illustrative examples are considered to demonstrate the effectiveness and accuracy of the proposed technique.

Figures

Tables

Table 1. Table 1 : Example  6.1 : error and convergence order for the state and control functions.

		$x\left(t\right)$			$u\left(t\right)$
$n$		$E_n\left(x\right)$	${ϵ}_n\left(x\right)$		$E_n\left(u\right)$	${ϵ}_n\left(u\right)$
$8$		$1.23e+0$	$2.34$		$3.10e+0$	$3.63$
$16$		$2.43e-1$	$3.09$		$2.51e-1$	$3.56$
$32$		$2.86e-2$	$3.41$		$2.13e-2$	$2.44$
$64$		$2.68e-3$	$3.50$		$3.92e-3$	$3.37$
$128$		$2.36e-4$	$3.52$		$3.79e-4$	$3.57$
$256$		$2.06e-5$	—		$3.18e-5$	—

Table 2. Table 2 : Example  6.1 : comparison of the numerical results obtained by the present method versus methods of [ 11 ] .

	Direct GL method [11]			Direct TR method [11]			Direct SI method [11]			Present method
$n$	$E_n\left(x\right)$	$E_n\left(u\right)$		$E_n\left(x\right)$	$E_n\left(u\right)$		$E_n\left(x\right)$	$E_n\left(u\right)$		$E_n\left(x\right)$	$E_n\left(u\right)$
$100$	$1.11e-1$	$1.68e-1$		$1.48e-2$	$2.07e-2$		$5.60e-4$	$8.99e-4$		$5.63e-4$	$9.03e-4$
$200$	$5.71e-2$	$9.19e-2$		$3.71e-3$	$5.21e-3$		$4.91e-5$	$7.66e-5$		$4.92e-5$	$7.68e-5$
$300$	$3.94e-2$	$6.37e-2$		$1.65e-3$	$2.32e-3$		$1.18e-5$	$1.80e-5$		$1.18e-5$	$1.80e-5$

Table 3. Table 3 : Numerical results for Example  6.2 .

Methods		$E_n\left(x\right)$	${ϵ}_n\left(x\right)$	$J$	CPU Time
Present	$n=2$	$8.07e-4$	$4.02$	$-0.3063957$	$0.015$
method	$n=4$	$4.99e-5$	$4.01$	$-0.3068248$	$0.015$
	$n=8$	$3.09e-6$	$4.01$	$-0.3068511$	$0.015$
	$n=16$	$1.92e-7$	$4.00$	$-0.3068527$	$0.032$
	$n=32$	$1.20e-8$	—	$-0.3068528$	$0.032$
Method of [21]	$M=3,$ $N=1$			$-0.30683$	$40.091$
Method of [20]	$M=3,$ $N=1$			$-0.30684$	$53.133$
Method of [19]	$M=3$			$-0.30685$

Table 4. Table 4 : Numerical results for Example  6.3 .

Methods		$E_n\left(x\right)$	${ϵ}_n\left(x\right)$	$E_n\left(u\right)$	${ϵ}_n\left(u\right)$	$J$	CPU Time
Present	$n=4$	$7.10e-4$	$3.39$	$2.98e-4$	$3.03$	$9.64314e-7$	$0.000$
method	$n=8$	$6.75e-5$	$3.33$	$3.65e-5$	$3.15$	$1.00418e-8$	$0.015$
	$n=16$	$6.69e-6$	$3.28$	$4.10e-6$	$3.18$	$1.06677e-10$	$0.046$
	$n=32$	$6.91e-7$	$3.22$	$4.52e-7$	$3.17$	$1.19487e-12$	$0.109$
	$n=64$	$7.42e-8$	$3.18$	$5.03e-8$	$3.15$	$1.41601e-14$	$0.765$
	$n=128$	$8.20e-9$	$3.15$	$5.66e-9$	$3.14$	$1.75827e-16$	$4.922$
	$n=256$	$9.24e-10$	—	$6.44e-10$	—	$2.25012e-18$	$21.282$
Method of [22]	$m=4$					$5.42028e-7$	$0.141$
	$m=8$					$8.22283e-10$	$0.296$
Method of [23]	$m=4$					$5.16864e-7$	$0.078$
	$m=8$					$4.23025e-11$	$0.171$

Equations

$$\min J={\int }_0^{t_f}f(t,x(t),u(t))dt$$

$$D^{\alpha }x(t)=g\left(t,x(t),D^{{\alpha }_1}x(t),D^{{\alpha }_2}x(t),\dots ,D^{{\alpha }_k}x(t),u(t)\right),m-1<\alpha \le m,$$

$$x^{(i)}(0)=q_i,i=0,1,\dots ,m-1,$$

$$\begin{array}{rl}& I^{\alpha }x(t)=\frac{1}{\Gamma (\alpha )}{\int }_0^t(t-s{)}^{\alpha -1}x(s)ds,\\ & I^{0}x(t)=x(t),\end{array}$$

$$D^{\alpha }x(t)=\frac{1}{\Gamma (m-\alpha )}{\int }_0^t(t-s{)}^{m-\alpha -1}{x}^{(m)}(s)ds,m-1<\alpha \le m.$$

$$D^{\alpha }{I}^{\alpha }x(t)=x(t),$$

$$I^{\alpha }{D}^{\alpha }x(t)=x(t)-\sum _{i=0}^{m-1}{x}^{(i)}(0)\frac{t^i}{i!},t>0,$$

$$I^{\alpha -\beta }{D}^{\alpha }x(t)=D^{\beta }x(t)-\sum _{i=\lceil \beta \rceil }^{m-1}{x}^{(i)}(0)\frac{t^{i-\beta }}{\Gamma (i-\beta +1)},0<\beta <\alpha ,t>0,$$

$$\psi_{0}(t)=\left\{\begin{array}[]{ll}\frac{1}{2h^{2}}(t-h)(t-2h),&0\leq t\leq 2h,\\ &\\ 0,&\text{otherwise};\end{array}\right.$$

$$\psi_{i}(t)=\left\{\begin{array}[]{ll}\frac{-1}{h^{2}}(t-(i-1)h)(t-(i+1)h),&(i-1)h\leq t\leq(i+1)h,\\ &\\ 0,&\text{otherwise};\end{array}\right.$$

$$\psi_{i}(t)=\left\{\begin{array}[]{ll}\frac{1}{2h^{2}}(t-(i-1)h)(t-(i-2)h),&(i-2)h\leq t\leq ih,\\ &\\ \frac{1}{2h^{2}}(t-(i+1)h)(t-(i+2)h),&ih\leq t\leq(i+2)h,\\ &\\ 0,&\text{otherwise},\end{array}\right.$$

$$\psi_{n}(t)=\left\{\begin{array}[]{ll}\frac{1}{2h^{2}}(t-(t_{f}-h))(t-(t_{f}-2h)),&t_{f}-2h\leq t\leq t_{f},\\ &\\ 0,&\text{otherwise}.\end{array}\right.$$

$$\psi_{i}(jh)=\left\{\begin{array}[]{ll}1,&i=j,\\ &\\ 0,&i\neq j,\end{array}\right.$$

$$\sum _{i=0}^n{\psi }_i(t)=1,$$

$$\int_{0}^{t_{f}}\psi_{i}(t)dt=\left\{\begin{array}[]{ll}\frac{h}{3},&\text{if $i=0,n$},\\ &\\ \frac{4h}{3}&\text{if $i$ is odd and $1\leq i\leq n-1$},\\ &\\ \frac{2h}{3},&\text{if $i$ is even and $2\leq i\leq n-2$}.\end{array}\right.$$

$$x(t)\simeq x_n(t)=\sum _{i=0}^n{a}_i{\psi }_i(t)=A^T\Psi (t),$$

$$\Psi (t)=[{\psi }_0(t),{\psi }_1(t),\dots ,{\psi }_n(t){]}^T,$$

$$A=[a_0,a_1,\dots ,a_n{]}^T,$$

$$I^{\alpha }\Psi (t)\simeq P^{(\alpha )}\Psi (t),$$

$$P^{(\alpha)}=\frac{h^{\alpha}}{2\Gamma(\alpha+3)}\left[\begin{array}[]{cccccccc}0&\beta_{1}&\beta_{2}&\beta_{3}&\beta_{4}&\ldots&\beta_{n-1}&\beta_{n}\\ 0&\eta_{0}&\eta_{1}&\eta_{2}&\eta_{3}&\ldots&\eta_{n-2}&\eta_{n-1}\\ 0&\xi_{-1}&\xi_{0}&\xi_{1}&\xi_{2}&\ldots&\xi_{n-3}&\xi_{n-2}\\ 0&0&0&\eta_{0}&\eta_{1}&\ldots&\eta_{n-4}&\eta_{n-3}\\ 0&0&0&\xi_{-1}&\xi_{0}&\ldots&\xi_{n-5}&\xi_{n-4}\\ \vdots&\vdots&\vdots&\vdots&\vdots&&\vdots&\vdots\\ 0&0&0&0&0&\ldots&\eta_{0}&\eta_{1}\\ 0&0&0&0&0&\ldots&\xi_{-1}&\xi_{0}\\ \end{array}\right]$$

$$\begin{array}{rl}& {\beta }_1=\alpha (3+2\alpha ),\\ & {\beta }_i=i^{\alpha +1}(2i-6-3\alpha )+2i^{\alpha }(1+\alpha )(2+\alpha )-(i-2{)}^{\alpha +1}(2i-2+\alpha ),i=2,3,\dots ,n,\\ & {\eta }_0=4(1+\alpha ),\\ & {\eta }_i=4[(i-1{)}^{\alpha +1}(i+1+\alpha )-(i+1{)}^{\alpha +1}(i-1-\alpha )],i=1,2,\dots ,n-1,\\ & {\xi }_{-1}=-\alpha ,\\ & {\xi }_0=2^{\alpha +1}(2-\alpha ),\\ & {\xi }_1=3^{\alpha +1}(4-\alpha )-6(2+\alpha ),\\ & {\xi }_i=(i+2{)}^{\alpha +1}(2i+2-\alpha )-6i^{\alpha +1}(2+\alpha )-(i-2{)}^{\alpha +1}(2i-2+\alpha ),i=2,3,\dots ,n-2.\end{array}$$

$$t={\Theta }^T\Psi (t),$$

$$\Theta =[0,h,2h,\dots ,t_f{]}^T.$$

$$D^{\alpha }x(t)=A^T\Psi (t)$$

$$u(t)\simeq u_n(t)=U^T\Psi (t),$$

$$A=[a_0,a_1,\dots ,a_n{]}^T,U=[u_0,u_1,\dots ,u_n{]}^T,$$

$$\begin{array}{rl}{D}^{{\alpha }_s}x(t)& =I^{\alpha -{\alpha }_s}(D^{\alpha }x(t))+\sum _{i=\lceil {\alpha }_s\rceil }^{m-1}{x}^{(i)}(0)\frac{t^{i-{\alpha }_s}}{\Gamma (j-{\alpha }_s+1)}\\ & \simeq A^T{I}^{\alpha -{\alpha }_s}\Psi (t)+\sum _{i=\lceil {\alpha }_s\rceil }^{m-1}{q}_i\frac{t^{i-{\alpha }_s}}{\Gamma (i-{\alpha }_s+1)}\\ & \simeq (A^T{P}^{(\alpha -{\alpha }_s)}+D_s^T)\Psi (t)\\ & =X_s^T\Psi (t),\end{array}$$

$$\sum _{i=\lceil {\alpha }_s\rceil }^{m-1}{q}_i\frac{t^{i-{\alpha }_s}}{\Gamma (i-{\alpha }_s+1)}\simeq D_s^T\Psi (t),X_s=(A^T{P}^{(\alpha -{\alpha }_s)}+D_s^T{)}^T$$

$$X_s=[x_{s,0},x_{s,1},\dots ,x_{s,n}{]}^T.$$

$$x(t)\simeq x_n(t)=X^T\Psi (t),$$

Full text

1 Introduction

Optimal control theory allows us to choose control functions in a dynamical system to achieve a certain goal. The theory generalizes that of the classical calculus of variations and has found many applications in engineering, physics, economics, and life sciences [1, 2]. In 1996/97, motivated by nonconservative physical processes in mechanics, the subject was extended to the case in which derivatives and integrals are understood as fractional operators of arbitrary order [3, 4].

Fractional optimal control is nowadays an important research area, allowing to apply the power of variational methods to real systems [5, 6, 7]. It is an expanding branch of applied mathematics that attracts more and more attention, having been applied in transportation, electronic, chemical and biological systems [8]. Another applications of fractional optimal control problems are shown in [9]. Besides providing nontrivial generalizations of the classical theory, it opens doors to new and interesting modern scientific problems [10, 11]. Given the complexity of the problems, numerical methods are indispensable [12, 13, 14].

Recently, numerical methods based on modified hat functions have shown to provide an effective way for solving fractional differential equations [15, 16]. Here, for the first time in the literature, we introduce such a method for solving a general class of fractional optimal control problems. More precisely, in this work we begin by considering the following fractional optimal control problem (FOCP):

[TABLE]

subject to the fractional dynamic system

[TABLE]

and the initial conditions

[TABLE]

where $\alpha>\alpha_{k}>\ldots>\alpha_{2}>\alpha_{1}>0$, $D^{\alpha}$ denotes the (left) Caputo derivative of order $\alpha$, $x(\cdot)$ is the state function, $u(\cdot)$ is the control function, $J$ is called the performance index, and $f$ and $g$ are given continuous functions. The aim of problem (1)–(3) consists in finding the optimal state and control functions that minimize the performance index. In this paper, a numerical method based on the modified hat functions is introduced for solving problem (1)–(3) and then the method is developed for FOCPs with inequality constraints. Because of the importance of fractional optimal control in applications, the numerical solution of FOCPs is crucial and our results can be very useful for many researchers.

The paper is organized as follows. In Section 2, some preliminaries of fractional calculus are reviewed and the basis functions and their properties are presented. A numerical method for solving the FOCP (1)–(3) is introduced in Section 3. In Section 4, an error estimate is proved for the proposed method. Section 5 is devoted to the development of the method to solve constrained FOCPs. In Section 6, some FOCPs are considered and solved using the proposed method. Finally, a conclusion is given in Section 7.

2 Preliminaries

In this section, a brief review on necessary definitions and properties of fractional calculus is presented. Furthermore, the modified hat functions and some of their properties are recalled.

2.1 Preliminaries of fractional calculus

We use fractional derivatives in the sense of Caputo. They are defined via Riemann–Liouville fractional integral.

Definition 2.1** (See, e.g., [17]).**

The Riemann–Liouville fractional integral operator with order $\alpha\geq 0$ of a given function $x(\cdot)$ is defined as

[TABLE]

where $\Gamma(\cdot)$ is the Euler gamma function.

Definition 2.2** (See, e.g., [17]).**

The Caputo fractional derivative of order $\alpha>0$ of a function $x(\cdot)$ is defined as

[TABLE]

For $m-1<\alpha\leq m$, $m\in\mathbb{N}$, some useful properties of the Caputo derivative and Riemann–Liouville integral are as follows:

[TABLE]

[TABLE]

[TABLE]

where $\lceil\cdot\rceil$ is the ceiling function.

2.2 Properties of the modified hat functions

By considering an even integer number $n$, the interval $[0,t_{f}]$ is divided into $n$ subintervals $[ih,(i+1)h]$, $i=0,1,2,\ldots,n-1$, with $h=\frac{t_{f}}{n}$. Then the generalized modified hat functions consist of a set of $n+1$ linearly independent functions in $L^{2}[0,t_{f}]$ that are defined as follows [15, 16]:

[TABLE]

if $i$ is odd and $1\leq i\leq n-1$,

[TABLE]

if $i$ is even and $2\leq i\leq n-2$,

[TABLE]

and

[TABLE]

These functions satisfy some interesting properties, which make them very useful. Some of these properties are as follows:

[TABLE]

[TABLE]

[TABLE]

Any function $x(\cdot)\in L^{2}[0,t_{f}]$ can be approximated using the modified hat functions as

[TABLE]

where

[TABLE]

and

[TABLE]

with $a_{i}=x(ih)$.

Let $\Psi(\cdot)$ be the modified hat functions basis vector given by (8) and $\alpha>0$. Then,

[TABLE]

where $P^{(\alpha)}$ is a matrix of dimension $(n+1)\times(n+1)$ called the operational matrix of fractional integration of order $\alpha$. This matrix is given (see [15, 16]) as

[TABLE]

with

[TABLE]

3 Numerical solution of problem (1)–(3)

In this section, we introduce a new numerical scheme for solving the optimal control problem (1)–(3). Our method is based on the properties of the modified hat functions. We begin by expanding the variable $t$ in terms of the modified hat functions. Using (7), this variable can then be rewritten as

[TABLE]

where

[TABLE]

Let us suppose that

[TABLE]

and

[TABLE]

where

[TABLE]

such that the elements of $A$ and $U$ are unknown. By utilizing (5), (9), (13), and the initial conditions given in (3), for $s=1,2,\ldots,k$, we obtain that

[TABLE]

where we have used the notations

[TABLE]

with

[TABLE]

Furthermore, the state function can be approximated using (4), (9), and (13), as

[TABLE]

where we have used

[TABLE]

with

[TABLE]

We rewrite the elements of the vectors $X$ and $X_{s}$, $s=1,2,\ldots,k$, in terms of the elements of the vector $A$. To this aim, suppose that

[TABLE]

Taking into consideration (10), we have

[TABLE]

Therefore, the entries $x_{j}$ and $x_{s,j}$ can be written in details as

[TABLE]

[TABLE]

We employ the approximations given above in order to obtain an approximation of the performance index. For the function $f:\mathbb{R}^{3}\rightarrow\mathbb{R}$, using (7), (11), (14), and (16), we can write

[TABLE]

where

[TABLE]

By substituting (18) into the performance index and employing (6), we obtain that

[TABLE]

where

[TABLE]

In a similar way, taking into account (7), (11), and (14)–(16) for the function $g:\mathbb{R}^{k+3}\rightarrow\mathbb{R}$, we have

[TABLE]

where

[TABLE]

Hence, using (13) and (21), the dynamic system given by (2) is reduced to

[TABLE]

Finally, according to the Lagrange multiplier method for minimizing (19) subject to the conditions given in (22), we define:

[TABLE]

where

[TABLE]

is the unknown Lagrange multiplier. The necessary optimality conditions are as follows:

[TABLE]

System (23) includes $3(n+1)$ nonlinear equations with $3(n+1)$ unknown parameters, which are the elements of $A$, $U$, and $\lambda$. By solving this system, approximations of the optimal control and corresponding state functions are given by (14) and (16), respectively. Furthermore, an approximation of the optimal value of the performance index is given by (19). The speed of our method depends on the speed of solving system (23). Therefore, the form of the system is very important. Here, we detail the closed forms for system (23). For simplicity, suppose that $g:=g(t,x(t),u(t))$. Then, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

As it can be seen, the partial derivatives of $f$ and $g$ with respect to the elements of $X$ and $U$ are needed to provide the final system. Symbolic computation can be used to produce these partial derivatives. Note that our method reduces the fractional optimal control problem (1)–(3) to the solution of a system of nonlinear algebraic equations. In some applications, solving the resulting system can be sensitive to the choice of initial values and may result in a singular Jacobian matrix. In that case, some extra care in choosing the initial values is needed.

4 Error estimate

In this section, we present an error estimate for the method proposed in the previous section. To this end, first let us recall the following theorem from [15].

Theorem 4.1** (See [15]).**

If a function $y(\cdot)\in C^{3}[0,t_{f}]$ is approximated using the set of modified hat functions by $y_{n}(t)=\sum\limits_{i=0}^{n}y(ih)\psi(t)$, then

[TABLE]

According to Theorem 4.1, if the control function $u(\cdot)$ belongs to the space $C^{3}[0,t_{f}]$, then we will have

[TABLE]

where $u_{n}(t)$ has been given in (14). For obtaining an error estimate for the state function $x(\cdot)$, we suppose that $D^{\alpha}x(\cdot)\in C^{3}[0,t_{f}]$ and is approximated using (13). Also, suppose that $D^{T}\Psi(t)$ is the approximation of function $\sum\limits_{i=0}^{m-1}q_{i}\frac{t^{i}}{i!}$. Then, using Theorem 4.1, we have

[TABLE]

Using (4) and the initial conditions in (3), we get

[TABLE]

Now, we consider $x_{n}(t)=X^{T}\Psi(t)$ as an approximation of the state function $x(\cdot)$, where $X$ is given by (17). Then, by neglecting the error of the operational matrix and utilizing (16), (24), and (25), we obtain that

[TABLE]

To give an error estimate for the performance index obtained by the method we propose in Section 3, we make use of the following theorem.

Theorem 4.2** (Composite Simpson’s rule [18]).**

Let $y(\cdot)\in C^{4}[0,t_{f}]$, $n$ be even, $h=t_{f}/n$, and $t_{j}=jh$, $j=0,1,\ldots,n$. Then, the composite Simpson’s rule can be written with its corresponding error as

[TABLE]

where $\mu\in(0,t_{f})$.

To complete the error discussion, we prove the following theorem.

Theorem 4.3**.**

Suppose that $f:\mathbb{R}^{3}\rightarrow\mathbb{R}$ is sufficiently continuously differentiable and satisfies a Lipschitz condition with respect to the second argument. Let $x^{*}(\cdot)\in C^{3}[0,t_{f}]$ and $u^{*}(\cdot)\in C^{3}[0,t_{f}]$ be the optimal state and control functions, respectively, which minimize the performance index $J$ with optimal value $J_{opt}$. If $x_{n}(\cdot)$, $u_{n}(\cdot)$, and $J_{n}$ are the approximation of the optimal state function, control function, and performance index given by our method, then

[TABLE]

Proof.

We put $t_{j}=jh$, $x_{j}:=x_{n}(t_{j})$, and $u_{j}:=u(t_{j})$. By setting $y(t):=f(t,x_{n}(t),u(t))$ and using Theorem 4.2, we have

[TABLE]

Since function $f$ satisfies a Lipschitz condition with respect to the second argument, there exists a real constant $C>0$ such that

[TABLE]

By assumption, we have

[TABLE]

and

[TABLE]

By subtracting (29) from (28) and using (26) and (27), we get

[TABLE]

The proof is complete. ∎

5 Numerical solution of FOCPs with inequality constraints

In this section, we develop the method introduced in Section 3 for solving constrained FOCPs. Consider problem (1)–(3) with an additional inequality constraint

[TABLE]

In order to solve this problem, we substitute the approximations given in (13)–(16) into (31) and get

[TABLE]

We set the following Newton–Cotes nodes as the collocation points:

[TABLE]

Therefore, we get

[TABLE]

Using (19), (22), and (32), the constrained FOCP is reduced to the following nonlinear programming problem:

[TABLE]

It should be noted that the decision variables of the above nonlinear programming problem are: $a_{i}$ and $u_{i}$, $i=0,1,\ldots,n$. By employing a standard optimization solver, problem (33) can be easily solved (see Section 6). After finding the optimal value of the decision variables, the approximations of the state and control functions are given, respectively, by (16) and (14). Moreover, we get an approximation of the optimal performance index by substituting the results in (19).

6 Illustrative examples

In this section, we apply our method to some FOCPs and compare the results with the ones obtained by existing methods in the literature. To this aim, the $l_{2}$ norm of the error and the convergence order of state and control functions are defined, respectively, by

[TABLE]

and

[TABLE]

where $t_{i}=ih$, $x(t_{i})$, and $u(t_{i})$ are the exact state and control functions at $t_{i}$ and $x_{n}(t_{i})$ and $u_{n}(t_{i})$ are the obtained state and control functions by the proposed method at $t_{i}$.

In our implementation, the method was carried out using Mathematica 11.3. For solving the resulting systems of algebraic equations, the function FindRoot was used in two examples, which are non-constrained FOCPs. Moreover, in the case of constrained FOCP, the function Minimize was employed for solving the linear programming problem. In order to have a comparison with other methods, we report the CPU time (seconds) in some of the examples, which have been obtained on a 2.5 GHz Core i7 personal computer with 16 GB of RAM.

Example 6.1**.**

As the first example, we consider the following FOCP borrowed from [11]:

[TABLE]

The optimal value for the performance index is $J=0$, which is obtained by

[TABLE]

where $J_{0}$ is the first kind Bessel function of order zero. This problem has been solved using the method we propose here, with different values of $n$. By applying the method with $n$ equal to some powers of $2$, we report in Table 1 the error of our numerical results in the $l_{2}$ norm and the order of convergence. Table 1 confirms that even when the state and control functions do not belong to $C^{3}[0,20]$ (as considered by Theorem 4.1), our method has high convergence order.

In [11], problem (34) is considered with the boundary condition

[TABLE]

This problem is solved in [11] by approximating the fractional integral by the Grünwald–Letnikov formula (direct GL method), trapezoidal formula (direct TR method), and Simpson formula (direct SI method). A comparison between the results obtained by our method and the aforementioned methods given in [11] is displayed in Table 2. The results show that our method has the same accuracy as the direct SI method. Furthermore, the approximate and exact state and control functions obtained by the present method with $n=50$ and $n=100$, together with the corresponding error functions, are plotted in Figure 1.

Example 6.2**.**

Consider the following constrained FOCP [19, 20, 21]:

[TABLE]

In this problem, we have $t_{f}=1$. By taking $n=2$, the value of $h$ is given as $h=\frac{1}{2}$. We set

[TABLE]

where

[TABLE]

with unknown parameters $a_{i}$, $u_{i}$, $i=0,1,2$, and

[TABLE]

[TABLE]

[TABLE]

Taking into account the initial condition in (35), we obtain that

[TABLE]

where

[TABLE]

By substituting (37) into the performance index $J$ and using (19), we get

[TABLE]

Moreover, using (36) and (37), the dynamic system in (35) is reduced to

[TABLE]

which gives the following system of equations:

[TABLE]

Finally, we take into account the inequality constraints in (35) at $\tau_{i}=\frac{i+1}{6}$, $i=0,1,\ldots,4$, to obtain

[TABLE]

In summary, the problem is reduced to the minimization of (38) subject to the conditions given in (39) and (40).

In the particular case when $\alpha=1$, the exact solution of the problem is given by

[TABLE]

and the optimal value of the performance index with seven significant digits is $J=-0.3068528$. By considering $\alpha=1$, the operational matrix of integration is given by

[TABLE]

Therefore, by substituting the elements of $P^{(1)}$ into (38)–(40), the following linear programming problem is produced:

[TABLE]

By solving this problem, we find

[TABLE]

Substituting $u_{i}=1$, $i=0,1,2$, into (36), the exact control function is obtained. The approximate and exact state functions are plotted in Figure 2. In this case, the approximate value of the performance index is $J_{2}=-0.3063957$.

We have employed our method for solving problem (35) with different values of $n$. For all the considered values of $n$, the method gives the exact control function. Table 3 displays the numerical results for $x(t)$, $\epsilon_{n}(x)$, $J$, and the CPU time consumed for solving the resulting linear programming problem. It is seen that the present method gives the exact solution of $J$, with seven significant digits, with $n=32$ and that the convergence order of the state function is $O(h^{4})$. In Table 3, the method of [21] refers to the hybrid block-pulse and Taylor polynomials method, the method of [20] denotes the method of hybrid block-pulse and Bernoulli polynomials, and the method of [19] refers to the Bernstein polynomials method.

Example 6.3**.**

Let us now consider the following nonlinear FOCP [22, 23]:

[TABLE]

The state and control functions that minimize the performance index $J$ are given by

[TABLE]

with optimal value $J_{opt}=0$. By applying our method with different values of $n$, we obtain the numerical results displayed in Table 4. In this table, we give a comparison between the results for $J$ and CPU times consumed for solving the final algebraic system resulted by our method, the method of [22], which uses Legendre polynomials as the basis functions and the Gauss quadrature rule for integrating the FOCP, and also the method of [23], which employs Bernoulli polynomials and the Gauss quadrature rule for integrating the FOCP. In [22] and [23], $m$ refers to the maximum order of the Legendre and Bernoulli polynomials, respectively, that have been used for approximating the state and control functions. It can be seen from Table 4 that the method we suggest here approximates the state and control functions with convergence of order $O(h^{3})$. In Figure 3 (left), the results for $E_{n}(x)$, $E_{n}(u)$, and $E_{n}(J)=|J_{opt}-J_{n}|$, obtained by employing our method for some selected values of $n$, are plotted in a logarithmic scale. This figure confirms the $O(h^{3})$ accuracy order of the state and control functions. It also shows that the accuracy order of the performance index is higher than $O(h^{3})$. Moreover, the CPU times of the method are plotted in Figure 3 (right).

7 Conclusion

We introduced a new numerical approach for solving fractional optimal control problems (FOCPs). Our scheme uses modified hat functions as basis functions and gives approximations of the control and state variables as linear combinations of such basis functions. The properties of the basis functions, together with those of the Caputo derivative and Riemann–Liouville integral, allow us to reduce the FOCP to a system of nonlinear algebraic equations, which greatly simplifies the problem. An error estimate of the method is proved. Moreover, a generalization is given for solving constrained FOCPs. The method is employed for solving three illustrative examples and the obtained results confirm the efficiency, accuracy, and high performance of our technique when compared with state of the art numerical schemes available in the literature. We claim that our method can be very useful in real applications. We are currently investigating its application to human respiratory syncytial virus infection and to the model introduced in [14]. Another direction of future research consists to generalize our method to variable-order FOCPs, which is an area under development [6, 24].

Acknowledgments

Lima acknowledges support from Fundação para a Ciência e a Tecnologia (FCT, the Portuguese Foundation for Science and Technology) through project UID/MAT/04621/2019; Torres was supported by FCT within project UID/MAT/04106/2019 (CIDMA). The authors are grateful to an anonymous reviewer for suggestions to improve the article.