On fractional optimal control problems with an application in fractional chaotic systems using a Legendre collocation-optimization technique

Abstract

In this paper, an intelligence method based on single layer legendre neural network is proposed to solve fractional optimal control problems where the dynamic control system depends on Caputo fractional derivatives. First, with the help of an approximation, the Caputo derivative is replaced to integer order derivative. According to the Pontryagin minimum principle for optimal control problems and by constructing an error function, an unconstrained minimization problem is then defined. In the optimization problem, trial solutions are used for state, costate and control functions, where these trial solutions are constructed by using Legendre polynomial based functional link artificial neural network. In the following, error back propagation algorithm is used for updating the network parameters (weights). At the end, some illustrative examples are included to demonstrate the validity and capability of the proposed method. Three applicable examples about chaos control of Malkus waterwheel, finance fractional chaotic models and fractional-order geomagnetic field models are also considered.

Keywords

Pontryagin minimum principle fractional optimal control problem Caputo fractional derivative Legendre neural network fractional chaotic system

Introduction

Fractional calculus is a classical mathematical notion and a generalization of the integer-order differentiation and integration to arbitrary non-integer order. The formation of fractional calculus is firstly projected by Leibniz in 1695 (Yong et al., 2016). Fractional calculus is developing faster during the last few years and many phenomena were described accurately with fractional models (see Podlubny, 1998; Kilbas et al., 2006; Machado et al., 2011; Magin, 2006; Oustaloup et al., 2000; Samko et al., 1993; Torvik and Bagley, 1984). For example, it has been described that materials with memory and hereditary effects and dynamical processes including gas diffusion and heat conduction in fractal porous media can be more adequately modeled by fractional order models (Zamani et al., 2007). Recently, studies implemented on new properties of real world problems including fractional differential equations have attracted attention of scientists from all over the world. Therefore, some new published papers are cited in Prakasha et al. (2019); Veeresha et al. (2019); Esen et al. (2018); Bulut et al. (2018a); Ravichandran et al. (2018); Gencoglu et al. (2017); Baskonus and Bulut (2015, 2016); Kumar et al. (2017); Baskonus et al. (2017); Bulut et al. (2018b); Dokuyucu et al. (2018); Sulaiman et al. (2018); and Yavuz et al. (2018).

If fractional differential equations contain a control variable and a performance index is given, a fractional optimal control problem (FOCP) is obtained that is interesting and very active in fractional systems. A FOCP can be defined with respect to different definitions of fractional differentials, but the most important types of derivatives are the Caputo and Riemann-Liouville. In most cases, the analytical solution of fractional problems may not be obtained easily and numerical methods are often applied, which have proven to be efficient and reliable (Alizadeh and Effati, 2018; Bhrawy et al., 2017; Dehghan et al., 2016; Ejlali and Hosseini, 2017; Ghasemi et al., 2017; Heydari et al., 2016; Hosseinpour and Nazemi, 2015, 2016; Jafari and Tajadodi, 2014; Lotfi et al., 2013; Rabiei et al., 2017; Rakhshan et al., 2018; Singha and Nahak, 2017; Tang et al., 2017; Tohidi and Nik, 2015). Among them, the authors in Singha and Nahak (2017) discuss a class of FOCP where the dynamical system comprises classical and fractional derivative. In this work, for numerical solution of the problem, authors use from Laguerre polynomials. In Rakhshan et al. (2018), a general formulation and solution scheme based on finding the numerical solution of the Hamilton-Jacobi-Bellman equation by the Legendre- Gauss collocation method is presented. In Tang et al. (2017), FOCP is solved by a new fractional pseudo spectral method. Although these methods provide good approximations to the solution, they require a discretization of domain via meshing, which may be challenging in two or higher dimension problems. The approximate solution derivatives are discontinuous and can seriously impact on the stability of the solution. In order to achieve satisfactory solution accuracy, it may be necessary to deal with finite meshes that significantly increase the computational cost. Also, some numerical methods use the operational matrices for solving FOCPs (see Heydari et al., 2016; Hosseinpour and Nazemi, 2015; Jafari and Tajadodi, 2014; Lotfi et al., 2013). However, finding the operational matrices is usually difficult and these matrices with high dimensions cause the complexity in computations. One promising approach for overcoming these limitations is to employ artificial neural networks based on a network topology, a connection pattern, neural activation properties, train strategy and ability to process data.

There are many references in theory and applications, modeling, design, structure and mathematics of neural networks (see Daniel, 2013; Haykin, 1994; Müller et al., 2012; Tang et al., 2007). In particular, the numerical solution of ordinary and partial differential equations (Beidokhti and Malek, 2009; Kumar and Yadav, 2011; Shirvany et al., 2009), optimal control problems (Cheng et al., 2007; Effati and Pakdaman, 2013; Vrabie and Lewis, 2009), mathematical programming problems (Nazemi, 2011, 2012, 2013, 2014; Nazemi and Effati, 2013; Nazemi and Nazemi, 2014; Nazemi and Omidi, 2012, 2013; Nazemi and Sharifi, 2013) and numerical solution of FOCP and fractional differential equations (Ghasemi and Nazemi, 2018; Ghasemi et al., 2017; Jafarian et al., 2017; Kheyrinataj and Nazemi, 2020a, 2020b; Yavari and Nazemi, 2019, 2020; Zuniga-Aguilar et al., 2018). Among the papers mentioned and in order to solve FOCPs, by using perceptron neural network’s ability in approximating a nonlinear function and the Grunwald–Letnikov approximation of Riemann–Liouville fractional derivatives, Ghasemi et al. (2017) contribute to presenting an indirect method for solving FOCPs. Kheyrinataj and Nazemi (2020b) have presented a new technique for some types of delay FOCPs, by means of a combination of the fractional power series method and ANNs. Yavari and Nazemi (2019) have also introduced a novel optimal control problem in fractional order as fractional infinite-horizon optimal control problems, where the dynamic control system depends on Caputo fractional derivatives, and solve it by a neural network scheme. Ghasemi and Nazemi (2018) use a computational intelligence method for solution of FOCPs with equality and inequality constraints based on a combination of the Mittag-Leffler function and an artificial neural network. Recently, Kheyrinataj and Nazemi (2020a) propose a higher order neural network, namely the functional link neural network, for the model of linear and nonlinear delay FOCPs with mixed control state constraints. Also Yavari and Nazemi (2020) introduce a novel fractional infinite horizon optimal control problem with a combination of conformable and Caputo–Fabrizio fractional derivatives and solve it using a network network scheme. In this paper, to reduce the complexity of some existing methods for solving FOCPs and to remove some disadvantages, a computational intelligence method based on single layer Legendre neural network is presented. This network model consists of single input node, one output layer and a functional expansion block based on Legendre polynomials. This is the basic motivation of this paper.

It should be noted that in fractional problems, for numerical purposes, approximations are used to deal with the fractional operators. For example, the Riemann-Liouville derivatives are expandable in a power series involving integer order derivatives only. For an analytic function, in Kilbas et al. (2006) an approximation for Riemann- Liouville derivative is defined. However, the obvious disadvantage of using this approximation is that in order to have a small error, one has to sum a large number of terms and thus the function has to possess higher order derivatives, which is not suitable for computation. In this paper for approximation of fractional derivatives, the approximation carried out in Almeida and Torres (2015) and Pooseh et al. (2013) is used, where a good approximation is obtained without the requirement of such higher-order smoothness on the admissible functions. Using this approximation the FOCP problem is converted into an OCP problem in integer order. The solution of the obtained OCP can be calculated either by using Pontryagin minimum principle (PMP), which provides a necessary condition for optimality, or by solving the Hamilton-Jacobi-Bellman (HJB) partial differential equation (PDE), which is a sufficient condition (see Kirk, 2012; Lewis and Syrmos, 1987). After imposing PMP to the considered OCP, a two-point boundary value problem (TPBVP) is solved by a Legendre neural network scheme.

The rest of this paper is recognized as follows: In Section 2, some basic definitions and results on fractional calculus are reviewed. In Section 3, mathematical modeling of the problem is introduced. The structure of feed-forward neural networks and the related necessary optimality condition are presented in Section 4. In Section 5, a learning algorithm for updating the neural network weights is suggested. In Section 6, the optimal control of the chaotic fractional systems is provided. Numerical simulations are also presented in Section 7 and Section 8 consists of brief summary and conclusion.

Fractional derivatives and integrals

In this section, some definitions in fractional calculus (Kilbas et al., 2006; Pooseh et al., 2013) are provided, which are required throughout the paper.

Definition 1: Let $f : [a, b] \to R$ , $α > 0$ be the order of the integral or the derivative and $n = [α] + 1$ . For $t \in [a, b]$ ,

• The left and right Riemann-Liouville fractional integral, respectively, by

{{}_{a}I}_{t}^{α} f (t) = \frac{1}{Γ (α)} \int_{a}^{t} (t - τ)^{α - 1} f (τ) d τ, (t > a),

(1)

{{}_{t}I}_{b}^{α} f (t) = \frac{1}{Γ (α)} \int_{t}^{b} (τ - t)^{α - 1} f (τ) d τ, (t < b) .

(2)

The left and right Riemann-Liouville fractional derivatives, respectively, by

{{}_{a}^{R}D}_{t}^{α} f (t) = \frac{1}{Γ (n - α)} {(\frac{d}{dt})}^{n} \int_{a}^{t} (t - τ)^{n - α - 1} f (τ) d τ, (t > a),

(3)

{{}_{t}^{R}D}_{b}^{α} f (t) = \frac{1}{Γ (n - α)} {(- \frac{d}{dt})}^{n} \int_{t}^{b} (τ - t)^{n - α - 1} f (τ) d τ, (t < b) .

(4)

The left and right Caputo fractional derivatives of $f (t)$ of order $α$ , respectively, by

{{}_{a}^{c}D}_{t}^{α} f (t) = \frac{1}{Γ (n - α)} \int_{a}^{t} (t - τ)^{n - α - 1} f^{(n)} (τ) d τ, (t > a),

(5)

{{}_{t}^{c}D}_{b}^{α} f (t) = \frac{{(- 1)}^{n}}{Γ (n - α)} \int_{t}^{b} (τ - t)^{n - α - 1} f^{(n)} (τ) d τ, (t < b),

(6)

where $f^{(n)} (τ)$ is the usual derivative of $f (τ)$ of order $n$ .

The following relations show the connection between Caputo and Riemann-Liouville fractional derivatives

{{}_{a}^{c}D}_{t}^{α} f (t) =_{a}^{R} D_{t}^{α} f (t) - \sum_{k = 0}^{n - 1} \frac{f^{(n)} (a)}{Γ (k + 1 - α)} (t - a)^{k - α},

(7)

{{}_{t}^{c}D}_{b}^{α} f (t) =_{t}^{R} D_{b}^{α} f (t) - \sum_{k = 0}^{n - 1} \frac{f^{(n)} (b)}{Γ (k + 1 - α)} (b - t)^{k - α} .

(8)

When $0 < α < 1$ ,

{{}_{a}^{c}D}_{t}^{α} f (t) =_{a}^{R} D_{t}^{α} f (t) - \frac{f (a)}{Γ (1 - α)} (t - a)^{- α},

(9)

{{}_{t}^{c}D}_{b}^{α} f (t) =_{t}^{R} D_{b}^{α} f (t) - \frac{f (b)}{Γ (1 - α)} (b - t)^{- α} .

(10)

the Left Riemann-Liouville fractional derivatives can be explained in the following way

\begin{matrix} {{}_{a}^{R}D}_{t}^{α} x (t) = A (α) (t - a)^{- α} x (t) + B (α) (t - a)^{1 - α} \overset{\cdot}{x} (t) \\ - \sum_{p = 2}^{\infty} C (α, p) (t - a)^{1 - p - α} y_{p} (t), \end{matrix}

(11)

where $α \in (0, 1)$ , $x \in C^{2} [a, b]$ and $y_{p} (t)$ for $p = 2, 3, \dots$ is the solution of the system

{\begin{matrix} {\overset{\cdot}{y}}_{p} (t) = (1 - p) (t - a)^{p - 2} y (t), \\ y_{p} (t) = 0 . \end{matrix}

(12)

$A (α)$ , $B (α)$ and $C (α)$ in (11) are also given by

{\begin{matrix} A (α) = \frac{1}{Γ (1 - α)} [1 + \sum_{p = 2}^{\infty} \frac{Γ (p - 1 + α)}{Γ (α) (p - 1)!}], \\ B (α) = \frac{1}{Γ (2 - α)} [1 + \sum_{p = 1}^{\infty} \frac{Γ (p - 1 + α)}{Γ (α - 1) p!}], \\ C (α, p) = \frac{Γ (p - 1 + α)}{Γ (2 - α) Γ (α - 1) (p - 1)!}, p = 2, 3, \dots . \end{matrix}

(13)

A similar formula can be deduced for the Caputo fractional derivative by using relationship between Riemann-Liouville and Caputo fractional derivatives. When finite sums are only considered, the following approximation is used

\begin{matrix} {{}_{a}^{R}D}_{t}^{α} x (t) \approx A (α, k) (t - a)^{- α} x (t) + B (α, k) (t - a)^{1 - α} \overset{\cdot}{x} (t) \\ - \sum_{p = 2}^{k} C (α, p) (t - a)^{1 - p - α} y_{p} (t), \end{matrix}

(14)

where $k \geq 2$ and

\begin{matrix} A (α, k) = \frac{1}{Γ (1 - α)} [1 + \sum_{p = 2}^{k} \frac{Γ (p - 1 + α)}{Γ (α) (p - 1)!}], \\ B (α, k) = \frac{1}{Γ (2 - α)} [1 + \sum_{p = 1}^{k} \frac{Γ (p - 1 + α)}{Γ (α - 1) p!}], \end{matrix}

(15)

The error is bounded by

| E_{tr} (t) | \leq max_{τ \in [a, t]} | \overset{\cdot\cdot}{x} (τ) | \frac{\exp ((1 - α)^{2} + 1 - α)}{Γ (2 - α) (1 - α) k^{1 - α}} (t - a)^{2 - α} .

(16)

For proofs and other details, see Pooseh et al. (2013).

Mathematical modeling of the problem

In this paper, a class of FOCPs with the Caputo fractional derivative is investigated as follows

\begin{matrix} minimize J (x, u) = \int_{a}^{b} L (t, x, u) dt, \\ subject to M \overset{\cdot}{x} (t) + N_{a}^{c} D_{t}^{α} x (t) = f (t, x, u), \\ x (a) = x_{a}, x (b) = x_{b}, \end{matrix}

(17)

where $x (t) \in R^{p}$ , $(M, N) \neq (0, 0)$ and $0 < α < 1$ . First, the operator ${{}_{a}^{c}D}_{t}^{α} x (t)$ is replaced with the help of relation (9) and approximation (14). For simplicity of our discussion, assume $A = A (α, k)$ , $B = B (α, k)$ and $C_{p} = C (α, p), p = 2, 3, \dots$ . Thus

\overset{\cdot}{x} (t) = \frac{f (t, x (t), u (t)) - NA {(t - a)}^{- α} x (t) + \sum_{p = 2}^{k} N C_{p} {(t - a)}^{1 - p - α} y_{p} (t) + \frac{N x_{a} {(t - a)}^{- α}}{Γ (1 - α)}}{M + NB {(t - a)}^{1 - α}} .

(18)

Define $Y (t) = (y_{2} (t), y_{3} (t), \dots, y_{k} (t))$ and

\begin{matrix} F (t, x (t), Y (t), u (t)) = \frac{f (t, x (t), u (t)) - NA {(t - a)}^{- α} x (t) + \sum_{p = 2}^{k} N C_{p} {(t - a)}^{1 - p - α} y_{p} (t) + \frac{N x_{a} {(t - a)}^{- α}}{Γ (1 - α)}}{M + NB {(t - a)}^{1 - α}} . \end{matrix}

The problem (17) is then converted into the following optimization problem

\begin{matrix} minimize J (x, u) = \int_{a}^{b} L (t, x, u) dt, \\ subject to \overset{\cdot}{x} (t) = F (t, x (t), Y (t), u (t)), \\ {\overset{\cdot}{y}}_{p} (t) = (1 - p) (t - a)^{p - 2} x (t), p = 2, \dots, k, \\ x (a) = x_{a}, \\ x (b) = x_{b}, \\ y_{p} (a) = 0, p = 2, \dots, k . \end{matrix}

(19)

According to the discussions in Nazemi and Karami (2017), if $(x, u)$ be a minimum solution of the problem (19), then there exist $λ (t)$ and $ϒ = (γ_{2}, \dots, γ_{k})$ such that $(x, Y, u, λ, ϒ)$ satisfies

{\begin{matrix} \begin{matrix} \frac{\partial H}{\partial x} (t, x, Y, u, λ, ϒ) = - \overset{\cdot}{λ} (t), \\ \frac{\partial H}{\partial λ} (t, x, Y, u, λ, ϒ) = \overset{\cdot}{x} (t), \\ \frac{\partial H}{\partial Y} (t, x, Y, u, λ, ϒ) = - \overset{\cdot}{ϒ} (t), \\ \frac{\partial H}{\partial ϒ} (t, x, Y, u, λ, ϒ) = \overset{\cdot}{Y} (t), \\ \frac{\partial H}{\partial u} (t, x, Y, u, λ, ϒ) = 0, \\ x (a) = x_{a}, \\ x (b) = x_{b}, \\ y_{p} (a) = 0, p = 2, \dots, k, \end{matrix} \end{matrix}

(20)

where $λ$ and $γ_{p}$ , $p = 2, \dots, k$ are the Lagrange multipliers and $H$ denotes the Hamiltonian function in the form

\begin{matrix} H (t, x, Y, u, λ, ϒ) = L (t, x, u) + λ F (t, x, Y, u) \\ + \sum_{p = 2}^{k} γ_{p} (1 - p) (t - a)^{p - 2} x (t) . \end{matrix}

(21)

In the next section, a neural network scheme for solving (20) is constructed.

Neural network construction and optimization

Artificial neural network methods have broad applicability to be a flexible and powerful scheme for the mathematical modeling and solving natural phenomena. Figure 1 shows the structure of single layer Legendre neural network, which consists of single input node, one output layer and a functional expansion block based on Legendre polynomials (Mall and Chakraverty, 2016). The hidden layer is eliminated by transforming the input pattern to a higher dimensional space using Legendre polynomials. In the presented architecture, the output value with input data $t$ and the weight vector $\bar{P}$ are defined as

{\begin{matrix} N (t, \bar{P}) = ν σ (z) \\ z = \sum_{i = 0}^{I} w_{i} L_{i} (t), \end{matrix}

(22)

where $L_{0} (t), L_{1} (t), \dots, L_{I} (t)$ are processing elements, their respective weights are $w_{i}, i = 0, 1, 2, \dots, I$ and $σ$ is an activation function as

σ (z) = \tanh (z) = \frac{e^{z} - e^{- z}}{e^{z} + e^{- z}} .

(23)

Figure 1.

A simplified block diagram of the Legendre neural network.

In order to solve numerically optimality system (20), consider neural networks for each functions: state, Lagrangian multipliers and the control as

{\begin{matrix} N_{x} (t, {\bar{P}}_{x}) = ν_{x} σ (z^{x}), z^{x} = \sum_{i = 0}^{I} w_{i}^{x} L_{i} (t), \\ N_{Y} (t, {\bar{P}}_{Y}) = ν_{Y} σ (z^{Y}), z^{Y} = \sum_{i = 0}^{I} w_{i}^{Y} L_{i} (t), \\ N_{λ} (t, {\bar{P}}_{λ}) = ν_{λ} σ (z^{λ}), z^{λ} = \sum_{i = 1}^{I} w_{i}^{λ} L_{i} (t), \\ N_{ϒ} (t, {\bar{P}}_{ϒ}) = ν_{ϒ} σ (z^{ϒ}), z^{ϒ} = \sum_{i = 0}^{I} w_{i}^{ϒ} L_{i} (t), \\ N_{u} (t, {\bar{P}}_{u}) = ν_{u} σ (z^{u}), z^{u} = \sum_{i = 0}^{I} w_{i}^{u} L_{i} (t), \end{matrix}

(24)

where $σ$ is an activation function as (23). It can be established a relationship between solution of the optimality conditions (20) and solution of an equivalent unconstrained minimization problem via the trial solutions (Dua, 2011; Kumar and Yadav, 2011). These trial functions are written as a sum of two parts and is constructed based on two facts: The first term contains no adjustable parameters and satisfies the initial/boundary conditions. The second term employs a single output feed forward neural network with adjustable parameters and input units fed with the input vector. It should be noted that in the second term, weights are adjusted to deal with the minimization problem and it is constructed so as not to contribute with initial/boundary conditions. Based on these facts and using the proposed neural networks in (24), the trial solutions (for state, Lagrangian multipliers and control) can be defined in the following structures

{\begin{matrix} x_{T} = x_{0} + (t - a) N_{x}, \\ Y_{T} = t N_{Y}, \\ λ_{T} = N_{λ}, \\ ϒ_{T} = (t - b) N_{ϒ}, \\ u_{T} = N_{u} . \end{matrix}

(25)

Note that it may be $λ (.) = 0$ for free end points. Thus, it can be defined $λ_{T}$ in (25) as $λ_{T} = (t - b) N_{λ}$ . For other initial (or boundary) conditions, appropriate trial functions can be constructed.

The trial solutions in (25) are the universal approximation and must satisfy condition (20). Consequently

{\begin{matrix} \frac{\partial H}{\partial x_{T}} (t, x_{T}, Y_{T}, u_{T}, λ_{T}, ϒ_{T}) = - {\overset{\cdot}{λ}}_{T} (t), \\ \frac{\partial H}{\partial λ_{T}} (t, x_{T}, Y_{T}, u_{T}, λ_{T}, ϒ_{T}) = {\overset{\cdot}{x}}_{T} (t), \\ \frac{\partial H}{\partial Y_{T}} (t, x_{T}, Y_{T}, u_{T}, λ_{T}, ϒ_{T}) = - {\overset{\cdot}{ϒ}}_{T} (t), \\ \frac{\partial H}{\partial ϒ_{T}} (t, x_{T}, Y_{T}, u_{T}, λ_{T}, ϒ_{T}) = {\overset{\cdot}{Y}}_{T} (t), \\ \frac{\partial H}{\partial u_{T}} (t, x_{T}, Y_{T}, u_{T}, λ_{T}, ϒ_{T}) = 0 . \end{matrix}

(26)

In order to reformulate (26) as an unconstrained minimization problem, the optimality system (26) can be collocated on the $m$ points $t_{l}, l = 1, . . ., m$ of the interval $[a, b]$ . General form of the corresponding error function for the ordinary differential equation (26) may be formulated as

\begin{matrix} {minimize}_{P} E (P) = \frac{1}{2} \sum_{l = 1}^{m} \\ {E_{1} (t_{l}, P) + E_{2} (t_{l}, P) + E_{3} (t_{l}, P) + E_{4} (t_{l}, P) + E_{5} (t_{l}, P)}, \end{matrix}

(27)

where $P = ({\bar{P}}_{x}, {\bar{P}}_{Y}, {\bar{P}}_{λ}, {\bar{P}}_{ϒ}, {\bar{P}}_{u})$ and

{\begin{matrix} E_{1} (t_{l}, P) = {[\frac{\partial H}{\partial x_{T}} (t, x_{T}, Y_{T}, u_{T}, λ_{T}, ϒ_{T}) + {\overset{\cdot}{λ}}_{T} (t)]}^{2}, \\ E_{2} (t_{l}, P) = {[\frac{\partial H}{\partial λ_{T}} (t, x_{T}, Y_{T}, u_{T}, λ_{T}, ϒ_{T}) - {\overset{\cdot}{x}}_{T} (t)]}^{2}, \\ E_{3} (t_{l}, P) = {[\frac{\partial H}{\partial Y_{T}} (t, x_{T}, Y_{T}, u_{T}, λ_{T}, ϒ_{T}) + {\overset{\cdot}{ϒ}}_{T} (t)]}^{2}, \\ E_{4} (t_{l}, P) = {[\frac{\partial H}{\partial ϒ_{T}} (t, x_{T}, Y_{T}, u_{T}, λ_{T}, ϒ_{T}) - {\overset{\cdot}{Y}}_{T} (t)]}^{2}, \\ E_{5} (t_{l}, P) = {[\frac{\partial H}{\partial u_{T}} (t, x_{T}, Y_{T}, u_{T}, λ_{T}, ϒ_{T})]}^{2} . \end{matrix}

(28)

Lemma 1: If $P = ({\bar{P}}_{x}, {\bar{P}}_{Y}, {\bar{P}}_{λ}, {\bar{P}}_{ϒ}, {\bar{P}}_{u})$ satisfies the following equation

\begin{matrix} η (P) = [E_{1} (t_{1}, P), \dots, E_{1} (t_{m}, P), E_{2} (t_{1}, P), \dots, E_{2} (t_{m}, P), \\ E_{3} (t_{1}, P), \dots, E_{3} (t_{m}, P), E_{4} (t_{1}, P), \dots, \\ E_{4} (t_{m}, P), E_{5} (t_{1}, P), \dots, E_{5} (t_{m}, P)] = 0, \end{matrix}

(29)

then $P^{*}$ is an optimal solution of (27).

Proof: Let $η (P^{*}) = 0 .$ Then $E_{i} (t_{l}, P^{*}) = 0,$ for $i = 1, \dots, 5$ and $l = 1, . . ., m .$ Since in (27), $E (P) \geq 0,$ thus $P^{*}$ is an optimal solution of (27).

Now by Lemma 1, it can be easily verified that the minimization problem (27) is equivalent to the following problem

{minimize}_{P} E (P) = \frac{1}{2} ‖ η (P) ‖^{2} .

(30)

In order to solve (30) as an unconstrained optimization problem, one can use some optimization and heuristic algorithms (see Bazaraa et al., 2013; Lee and El-Sharkawi, 2008; Nocedal and Wright, 2006; Zhang, 2013).

A learning algorithm

Our aim is to determine trial solutions in the form of (25) such that

E (P) = \frac{1}{2} ‖ η (P) ‖^{2},

(31)

is minimized. Here, an unsupervised version of back propagation algorithm has been used for learning and for updating the weights of trial solutions. Using the gradient descent method, weights are initialized randomly and those are updated as follows

P_{i} (j + 1) = P_{i} (j) - κ \frac{\partial E (P)}{\partial P_{i}},

where $κ$ is the learning parameter lies between 0 and 1, $j$ is the iteration step, $P_{i}, i = 1, 2, \dots, 5 (I + 2)$ are all adjustable parameters (weights) in trial solutions (25). For instance, unsupervised error back propagation algorithm of the $N_{x}$ , for updating the weights is as follows

w_{i}^{x} (j + 1) = w_{i}^{x} (j) - κ \frac{\partial E (P)}{\partial w_{i}^{x}}, i = 0, 1, \dots, I,

(32)

ν_{x} (j + 1) = ν_{x} (j) - κ \frac{\partial E (P)}{\partial ν_{x}} .

(33)

Similarly, other weights used to develop the $N_{λ}, N_{Y}, N_{ϒ}$ and $N_{u}$ are optimized using (32). The flowchart of the neural network learning algorithm is shown in Figure 2. In the following steps, the learning algorithm is described:

Step 1: Randomly select the initial values of trial solutions parameters $P_{i}, i = 1, 2, \dots, 5 (I + 2)$ and select an error tolerance parameter $ε > 0$ .

Step 2: Initialize the input vector $t = (t_{1}, t_{2}, \dots, t_{m})$ .

Step 3: Compute the output values of $N_{x}, N_{λ}, N_{Y}, N_{ϒ}$ and $N_{u}$ .

Step 4: Calculate the error function $E (P)$ using (31).

Step 5: Update the parameters using unsupervised back propagation algorithm

P_{i} (j + 1) = P_{i} (j) - κ \frac{\partial E (P)}{\partial P_{i}} .

Step 6: If the error function $E (P) \leq ε$ or iterations >1000, then go to step 7, otherwise go to step 2.

Step 7: After completing the learning algorithm, the final parameters are stored and then the converged neural networks can be used for testing.

Figure 2.

Flowchart of the neural network learning algorithm.

Theorem 1: Suppose that ${P (j), j = 1, 2, . . .}$ is the constructed sequence of the proposed learning algorithm in which the initial point is $P (0)$ and the level set $L (P_{0}) = {P (j) | E (P (j)) \leq E (P (0))}$ is bounded. Then

$γ^{+} (P_{0}) = {P (j), j = 1, 2, . . .}$ is bounded.

There exists $\hat{P}$ such that $lim_{j \to \infty} P (j) = \hat{P} .$

Proof:

Since in the learning scheme, the weight parameters and biases are optimized by the gradient descent algorithm, thus $E (P (j))$ in along ${P (j), j = 1, 2, . . .}$ is monotone nonincreasing. Therefore $γ^{+} (P (0)) \subseteq L (P (0)$ , that is to say $γ^{+} (P (0)) = {P (j) | j = 1, 2, . . .}$ is bounded.

By (a), $γ^{+} (P (0)) = {P (j) | j = 1, 2, . . .}$ is a bounded set of points. Thus there exists limiting point $\hat{P}$ , and there exists an index set $J$ such that ${P (j)}_{J} \to \hat{P}, j \in J$ which indicates that $\hat{P}$ is $ω$ -limit point of $γ^{+} (P (0))$ . Using the LaSalle invariant set theorem for discrete time dynamical systems (see Mei and Bullo, 2017), one has that ${P (j)}_{J} \to \hat{P} \in M$ as, where $M$ is the largest invariant set in $K = {P (j) | | E (P (j + 1)) - E (P (j)) | \leq ϵ}$ .

Optimal control of fractional chaotic systems

Chaos, as a very interesting phenomena, has complex dynamical behaviors, such as depending sensitively on tiny variations of initial conditions, having bounded trajectories in the phase space, etc. On the other side, the control of chaos, or control of chaotic systems, is the boundary field between control theory and dynamical systems theory studying when and how it is possible to control systems exhibiting irregular, chaotic behavior (Yu et al., 1998). Chaos control of chaotic system are initially put forward by Ott et al. (1990). As a comparison with the presented method in Baskonus et al. (2015), a Lyapunov function for the controlled fractional order autonomous system fractional order chaotic system was first constructed. According to the Lyapunov theory, to any equilibrium point, an active controller is found. Taking into account the described controllers, the equilibrium points are stabilized, and then, the chaos is controlled in original system of fractional order autonomous system. However, construction a Lyapunov function and finding the desired control is a routine method in chaos control. Moreover, constructing a suitable Lyapunov function in some cases for some control systems is difficult. In recent years, chaos control has been attracted increasingly attentions due to their potential applications in the fields of secure communications, optical, chemical, physical, biological systems, neural networks, and so forth (Das et al., 2000; Rong and Xiaoning, 1998). Many techniques have used about achieving chaos control of fractional order systems. For some recent works, Mahmoud et al. (2017) utilize the time delayed feedback control to control the integer and fractional version of Burke-Shaw chaotic system. Bigdeli and Ziazi (2017) precedes chaos control of fractional-order chaotic systems in presence of uncertainty and external disturbances based on a hybrid adaptive intelligent backstepping-sliding mode controller. Yadav et al. (2017) perform the stability analysis, chaos control and function projective synchronization between a fractional order chaotic chemical reactor system and a chaotic chemical reactor system with uncertain parameters with the help of a new lemma for the Caputo derivative and the Lyapunov stability theory. Avalos-Ruiz et al. (2018) present the simulation and control implementation on a field programmable gate array (FPGA) for a class of variable-order fractional chaotic systems by using sliding mode control strategy. Al-khedhairi et al. (2019) introduce a novel linear control technique based on Lyapunov stability theory to stabilize the chaotic states of the fractional-order geomagnetic field model to its five equilibrium points. Eshaghi et al. (2020) use the linear feedback control technique to control and stabilize the system to equilibrium point in order to eliminate the chaotic vibration.

In order to handle neural networks to chaos control of fractional system, Dehghani and Khanlo (2019) propose an adaptive chaos control based on the inverse dynamic method and radial basis function network, for a typical vibratory piezomagnetoelastic energy harvesting system to return the chaotic behavior to a periodic one. Lu and Wang (2019) address an observer-based command filtered adaptive neural network tracking control problem for a fractional-order chaotic permanent magnet synchronous motor (PMSM) with the immeasurable state, parameter uncertainties, and external load disturbance. Li and Sun (2019) consider the control of uncertain fractional-order Chua-Hartley chaotic systems by means of adaptive neural network backstepping control. However, until now, an infinite time optimal control approach has been used to to control of fractional chaotic systems. It is also denoted that Legendre neural networks have not used for the control of fractional-order chaotic systems yet. In this section, the FOCP of the chaotic systems about its equilibrium states using an infinite horizon optimal control strategy is studied.

Consider a dynamic system modeled by the state equations

{\begin{matrix} {{}_{0}^{C}D}_{t}^{α} x_{i} (t) = f_{i} (x_{1} (t), x_{2} (t), \dots, x_{n} (t)), i = 1, 2, \dots, s, \\ x_{i} (0) = x_{i, 0}, \end{matrix}

(34)

where $x = (x_{1}, x_{2}, \dots, x_{s}) : [0, \infty) \to A \subseteq R^{s}$ is the state variable and $f_{i}$ is a given nonlinear continuous functions. From Petráš (2011), it is clear that, in order to control of chaos in the chaotic system (34), it must be determined the control functions $u_{i}, i = 1, \dots, s$ so that

{\begin{matrix} {{}_{0}^{C}D}_{t}^{α} x_{i} (t) = f_{i} (x_{1} (t), x_{2} (t), \dots, x_{n} (t)) + u_{i} (t), i = 1, 2, \dots, s, \\ lim_{t \to \infty} x_{i} (t) = {\bar{x}}_{i}, \end{matrix}

(35)

where ${\bar{x}}_{i}, i = 1, 2, \dots, s$ is equilibrium state and $u = (u_{1}, u_{2}, \dots, u_{s}) : [0, \infty) \to B \subseteq R^{s}$ is control input which satisfies the conditions of the dynamical system (34) about its equilibrium points. It can be suggested a fractional infinite horizon optimal control (FIHOC) problem as

minimize J = \int_{0}^{\infty} \sum_{i = 0}^{n} (x_{i} (t) - {\bar{x}}_{i})^{2} dt,

(36)

\begin{matrix} subject to \\ {\begin{matrix} {{}_{0}^{C}D}_{t}^{α} x_{i} (t) = f_{i} (x_{1} (t), x_{2} (t), \dots, x_{n} (t)) + u_{i} (t), i = 1, 2, \dots, s, \\ x_{i} (0) = x_{i, 0} . \end{matrix} \end{matrix}

(37)

It is now shown that the optimal control problem (36) and (37) can solve the controlled chaotic system (35).

Definition 1: Input function $u (.) = \hat{u} (t), t \geq 0$ is control of the chaotic system (34) to unstable equilibrium point $\bar{x}$ , if the corresponding output $x (.)$ of system (34) satisfies $lim_{t \to \infty} x (t) = \bar{x}$ (i.e. $lim_{t \to \infty} ‖ x (t) - \bar{x} (t) ‖_{2} = 0$ ).

Theorem 2: Assume that there exists a solution $u (.) = \hat{u} (t)$ for control of the chaotic system (34) such that its corresponding output $x (.) = \hat{x} (t)$ is in $L_{2} ([0, \infty))$ . Consider FIHOC problem (36) and (37). If $(x^{*} (.), u^{*} (.))$ be an optimal solution of FIHOC problem (36) and (37), then optimal control $u^{*} (.)$ reaches the chaotic system (34) to its unique unstable equilibrium point.

Proof: Since $u (.) = \hat{u} (t)$ is a solution for control of the chaotic system (34), its corresponding output $x (.) = \hat{x} (t)$ satisfies $lim_{t \to \infty} \hat{x} (t) = \bar{x} (t)$ . So

lim_{t \to \infty} ‖ \hat{x} (t) - \bar{x} (t) ‖_{2}^{2} = ‖ lim_{t \to \infty} (\hat{x} (t) - \bar{x} (t)) ‖_{2}^{2} = 0,

(38)

is finite. On the other hand, $(x^{*} (.), u^{*} (.))$ is an optimal solution of FIHOC problem (36) and (37). Thus

\begin{matrix} 0 \leq J (x^{*} (.)) = ‖ x^{*} (.) - \bar{x} (.) ‖_{L_{2}}^{2} \\ = \int_{0}^{\infty} ‖ x^{*} (.) - \bar{x} (.) ‖_{2}^{2} dt \leq J (\hat{x} (.)) . \end{matrix}

(39)

Hence $J (\hat{x} (.))$ is finite. Now assume that $lim_{t \to \infty} x^{*} (t) = K, \bar{x} \neq K \in R^{n}$ . Then $lim_{t \to \infty} ‖ x^{*} (t) - \bar{x} (t) ‖_{2} = ‖ K - \bar{x} (t) ‖_{2}$ and by (39), $‖ K - \bar{x} (t) ‖_{2}$ must be a finite number. Hence, there are $n_{1} \in$ and $t_{1} > 0$ such that for all $t \geq t_{1},$ $‖ x^{*} (t) - \bar{x} (t) ‖_{2} \geq ‖ K - \bar{x} (t) ‖_{2} - \frac{1}{n_{1}}$ . Thus

\begin{matrix} J (x^{*} (.)) = \int_{0}^{\infty} ‖ x^{*} (t) - \bar{x} (t) ‖_{2}^{2} dt \\ = \int_{0}^{t_{1}} ‖ x^{*} (t) - \bar{x} (t) ‖_{2}^{2} dt + \int_{t_{1}}^{\infty} ‖ x^{*} (t) - \bar{x} (t) ‖_{2}^{2} dt \\ \geq \int_{0}^{t_{1}} ‖ x^{*} (t) - \bar{x} (t) ‖_{2}^{2} dt + (‖ K - \bar{x} (t) ‖_{2} - \frac{1}{n_{1}}) \int_{t_{1}}^{\infty} dt . \end{matrix}

(40)

By the last integral in the above relation, $J (x^{*} (.))$ cannot be finite which is a contradiction. Thus, $lim_{t \to \infty} x^{*} (t) = \bar{x} (t)$ .

Numerical examples

This section is devoted to illustrate the accuracy of the proposed numerical approach by several examples and application in fractional chaotic systems.

Example 1: Consider an optimization problem as follows

minimize J (x, u) = \int_{0}^{1} (u^{2} (t) - 4 x (t))^{2} dt,

subject to the dynamic constraint

\overset{\cdot}{x} (t) +_{0}^{c} D_{t}^{0.5} x (t) = u (t) + \frac{2}{Γ (2.5)} t^{1.5}, t \in [0, 1],

and the boundary conditions

x (0) = 0 and x (1) = 1 .

The exact solution for $α = 0.5$ in Almeida and Torres (2015) is the pair

(x (t), u (t)) = (t^{2}, 2 t) .

Using approximation (18), the problem is transformed into as

minimize J (x, u) = \int_{0}^{1} (u^{2} (t) - 4 x (t))^{2} dt,

subject to the dynamic constraint

{\begin{matrix} \begin{matrix} \overset{\cdot}{x} (t) = \frac{u (t) + \frac{2}{Γ (2.5)} t^{1.5} - A t^{- 0.5} x (t) + \sum_{p = 2}^{k} c_{p} t^{0.5 - p} y_{p} (t)}{1 + B t^{0.5}}, \\ {\overset{\cdot}{y}}_{p} (t) = (1 - p) t^{p - 2} x (t), p = 2, \dots, k, \end{matrix} \end{matrix}

and the boundary conditions

{\begin{matrix} \begin{matrix} x (0) = 0, x (1) = 1, \\ y_{p} (0) = 0, p = 2, \dots, k . \end{matrix} \end{matrix}

(41)

According to the initial conditions (41) with $k = 3,$ trial solutions may be selected as

{\begin{matrix} x_{T} = t + t (1 - t) N_{x}, \\ λ_{T} = N_{λ}, \\ y_{2_{T}} = t N_{y_{2}}, \\ y_{3_{T}} = t N_{y_{3}}, \\ γ_{2_{T}} = (1 - t) N_{γ_{2}}, \\ γ_{3_{T}} = (1 - t) N_{γ_{3}}, \\ u_{T} = N_{u} . \end{matrix}

Using the presented scheme in this article, the optimal trajectory, control and the absolute error for $α = 0.5$ are plotted in Figures 3 and 4, respectively. In Figure 5, the approximate solutions for several values of $α$ and the exact solution are depicted. Also, the values of the error function $E (P)$ are shown in the Table 1 for several values of $0 < α < 1$ .

Figure 3.

Left (the trajectory $x (t)$ for $α = 0.5$ ), right (the absolute error of $x (t)$ for $α = 0.5$ ).

Figure 4.

Left (the control $u (t)$ for $α = 0.5$ ), right (the absolute error of $u (t)$ for $α = 0.5$ ).

Figure 5.

The optimal trajectory $x (.)$ and the optimal control $u (.)$ with different values of $α$ .

Table 1.

The values of the error function $E (P)$ for several values $0 < α < 1$ .

$α$	0.3	0.5	0.7	0.9
$E$	$7.33 \times 10^{- 3}$	$8.75 \times 10^{- 4}$	$1.94 \times 10^{- 3}$	$1.56 \times 10^{- 3}$

Example 2: Consider the following FOCP

m inimize J (x, u) = \int_{0}^{1} (tu (t) - (α + 2) x (t))^{2} dt,

subject to dynamical system and boundary conditions

{\begin{matrix} \overset{\cdot}{x} (t) +_{0}^{c} D_{t}^{α} x (t) = u (t) + t^{2}, \\ x (0) = 0, x (1) = \frac{2}{Γ (3 + α)} . \end{matrix}

The exact solution for $α = 0.5$ is given in Sweilam and Al-Ajami (2015) as

(x (t), u (t)) = (\frac{2 t^{(α + 2)}}{Γ (α + 3)}, \frac{2 t^{α + 1}}{Γ (α + 2)}) .

Using approximation (18), the problem is transformed into as

minimize J (x, u) = \int_{0}^{1} (tu (t) - (α + 2) x (t))^{2} dt,

subject to dynamical system

{\begin{matrix} \begin{matrix} \overset{\cdot}{x} (t) = \frac{u (t) + t^{2} - A t^{- α} x (t) + \sum_{p = 2}^{k} c_{p} t^{1 - p - α} y_{p} (t)}{1 + B t^{1 - α}}, \\ {\overset{\cdot}{y}}_{p} (t) = (1 - p) t^{p - 2} x (t), p = 2, \dots, k, \end{matrix} \end{matrix}

and boundary conditions

{\begin{matrix} \begin{matrix} x (0) = 0, & x (1) = \frac{2}{Γ (3 + α)}, \\ y_{p} (0) = 0, & p = 2, \dots, k . \end{matrix} \end{matrix}

(42)

According to the initial conditions in (42) with $k = 3,$ trial solutions may be selected as

{\begin{matrix} x_{T} = \frac{2}{Γ (3 + α)} t + t (1 - t) N_{x}, \\ λ_{T} = N_{λ}, \\ y_{2_{T}} = t N_{y_{2}}, \\ y_{3_{T}} = t N_{y_{3}}, \\ γ_{2_{T}} = (1 - t) N_{γ_{2}}, \\ γ_{3_{T}} = (1 - t) N_{γ_{3}}, \\ u_{T} = N_{u} . \end{matrix}

Using the suggested approach, the optimal trajectory, control and the absolute error for $α = 0.2, 0.3, 0.5$ and $0.7$ are depicted in Figures 6 and 7, respectively.

Figure 6.

Left (the trajectory $x (t)$ for $α = 0.2, 0.3, 0.5, 0.7$ ), right (the absolute error of $x (t)$ for $α = 0.2, 0.3, 0.5, 0.7$ ).

Figure 7.

Left (the control $u (t)$ for $α = 0.2, 0.3, 0.5, 0.7$ ), right (absolute error of $u (t)$ for $α = 0.2, 0.3, 0.5, 0.7$ ).

Example 3. (Petráš, 2011) Consider the fractional chaotic system as

{\begin{matrix} {{}_{0}^{C}D}_{t}^{α} X_{1} (t) = - X_{1} (t) - X_{2}^{2} (t), \\ {{}_{0}^{C}D}_{t}^{α} X_{1} (t) = 2.5 X_{2} (t) - 4 X_{1} (t) X_{3} (t), \\ {{}_{0}^{C}D}_{t}^{α} X_{1} (t) = - 5 X_{3} (t) + 4 X_{1} (t) X_{2} (t), \\ X_{1} (0) = 0.2, X_{2} (0) = 0, X_{3} (0) = 0.5 . \end{matrix}

(43)

In order to control of chaos in the chaotic system (43), it must be determined the controls $U_{1} (.), U_{2} (.)$ and $U_{3} (.)$ so that

{\begin{matrix} {{}_{0}^{C}D}_{t}^{α} X_{1} (t) = - X_{1} (t) - X_{2}^{2} (t) + U_{1} (t), \\ {{}_{0}^{C}D}_{t}^{α} X_{2} (t) = 2.5 X_{2} (t) - 4 X_{1} (t) X_{3} (t) + U_{2} (t), \\ {{}_{0}^{C}D}_{t}^{α} X_{3} (t) = - 5 X_{3} (t) + 4 X_{1} (t) X_{2} (t) + U_{3} (t), \\ \lim_{t \to \infty} X_{i} (t) = {\bar{X}}_{i}, i = 1, 2, 3, \end{matrix}

(44)

where ${\bar{X}}_{i}$ are equilibrium states. This is equivalent to solving the following infinite horizon optimization problem

\begin{array}{l} minimize J = \int_{0}^{\infty} [{(X_{1} (t) - {\bar{X}}_{1})}^{2} + {(X_{2} (t) - {\bar{X}}_{2})}^{2} + {(X_{3} (t) - {\bar{X}}_{3})}^{2}] dt, \\ subject to {\begin{matrix} {{}_{0}^{C}D}_{t}^{α} X_{1} (t) = - X_{1} (t) - X_{2}^{2} (t) + U_{1} (t), \\ {{}_{0}^{C}D}_{t}^{α} X_{2} (t) = 2.5 X_{2} (t) - 4 X_{1} (t) X_{3} (t) + U_{2} (t), \\ {{}_{0}^{C}D}_{t}^{α} X_{3} (t) = - 5 X_{3} (t) + 4 X_{1} (t) X_{2} (t) + U_{3} (t), \\ X_{1} (0) = 0.2, X_{2} (0) = 0, X_{3} (0) = 0.5 . \end{matrix} \end{array}

(45)

It is clear that interval $[0, \infty)$ can be transformed into $[0, 1)$ , by a suitable change of variable

t = \frac{τ}{(1 - τ)} .

(46)

Now by (46), approximation (18) and assuming

{\begin{matrix} x_{i} (τ) = X_{i} (\frac{τ}{1 - τ}), & i = 1, 2, 3, \\ u_{i} (τ) = U_{i} (\frac{τ}{1 - τ}), & i = 1, 2, 3, \\ y_{i} (τ) = Y_{i} (\frac{τ}{1 - τ}), & i = 1, 2, 3, \end{matrix}

the problem (45) for $p = 2, 3$ is reduced to the following problem

\begin{matrix} minimize J = \int_{0}^{1} \frac{1}{{(1 - τ)}^{2}} \\ [(x_{1} (τ) - {\bar{X}}_{1})^{2} + (x_{2} (τ) - {\bar{X}}_{2})^{2} + (x_{3} (τ) - {\bar{X}}_{3})^{2}] d τ, \\ subject to {\begin{matrix} {\overset{\cdot}{x}}_{1} = \frac{1}{{(1 - τ)}^{2} B {(\frac{τ}{1 - τ})}^{1 - α}} (C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) \\ - (A {(\frac{τ}{1 - τ})}^{- α} + 1) x_{1} (τ) + x_{2}^{2} (τ) + u_{1} (τ)), \\ {\overset{\cdot}{x}}_{2} = \frac{1}{{(1 - τ)}^{2} B {(\frac{τ}{1 - τ})}^{1 - α}} (C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) \\ + (A (\frac{τ}{1 - τ} - 2.5)^{- α}) x_{2} (τ) - 4 x_{1} (τ) - x_{1} (τ) x_{3} (τ) + u_{2} (τ)), \\ {\overset{\cdot}{x}}_{3} = \frac{1}{{(1 - τ)}^{2} B {(\frac{τ}{1 - τ})}^{1 - α}} (C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) \\ - (A {(\frac{τ}{1 - τ})}^{- α} + 5) x_{3} + 4 x_{1} (τ) x_{2} (τ) + u_{3} (τ)), \\ {\overset{\cdot}{y}}_{21} (t) = - \frac{1}{{(1 - τ)}^{2}} x_{1} (τ), \\ {\overset{\cdot}{y}}_{22} (t) = - \frac{1}{{(1 - τ)}^{2}} x_{2} (τ), \\ {\overset{\cdot}{y}}_{23} (t) = - \frac{1}{{(1 - τ)}^{2}} x_{3} (τ), \\ {\overset{\cdot}{y}}_{31} (t) = - \frac{2 τ}{{(1 - τ)}^{3}} x_{1} (τ), \\ {\overset{\cdot}{y}}_{32} (t) = - \frac{2 τ}{{(1 - τ)}^{3}} x_{1} (τ), \\ {\overset{\cdot}{y}}_{33} (t) = - \frac{2 τ}{{(1 - τ)}^{3}} x_{1} (τ), \\ x_{1} (0) = 0.2, x_{2} (0) = 0, x_{3} (0) = 0.5, \\ y_{21} (0) = 0, y_{22} (0) = 0, y_{23} (0) = 0, \\ y_{31} (0) = 0, y_{32} (0) = 0, y_{33} (0) = 0 . \end{matrix} \end{matrix}

With $h = \frac{1}{{(1 - τ)}^{2}}$ and $K = \frac{1}{{(1 - τ)}^{2} B {(\frac{τ}{1 - τ})}^{1 - α}}$ , Hamiltonian function is written as

\begin{matrix} H = h [(x_{1} (τ) - {\bar{X}}_{1})^{2} + (x_{2} (τ) - {\bar{X}}_{2})^{2} + (x_{3} (τ) - {\bar{X}}_{3})^{2}] \\ + λ_{1} K [C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) \\ - (A {(\frac{τ}{1 - τ})}^{- α} + 1) x_{1} (τ) + x_{2}^{2} (τ) + u_{1} (τ)] \\ + λ_{2} K [C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) \\ + (A {(\frac{τ}{1 - τ})}^{- α} - 2.5) x_{2} (τ) - 4 x_{1} (τ) x_{3} (τ) + u_{2} (τ)] \\ + λ_{3} K [C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) \\ - (A {(\frac{τ}{1 - τ})}^{- α} + 5) x_{3} + 4 x_{1} (τ) x_{2} (τ) + u_{3} (τ)] \\ + γ_{21} [- \frac{1}{{(1 - τ)}^{2}} x_{1} (τ)] + γ_{22} [- \frac{1}{{(1 - τ)}^{2}} x_{2} (τ)] \\ + γ_{23} [- \frac{1}{{(1 - τ)}^{2}} x_{3} (τ)] \\ + γ_{31} [- \frac{2 τ}{{(1 - τ)}^{3}} x_{1} (τ)] + γ_{32} [- \frac{2 τ}{{(1 - τ)}^{3}} x_{1} (τ)] \\ + γ_{33} [- \frac{2 τ}{{(1 - τ)}^{3}} x_{1} (τ)] . \end{matrix}

Therefore, the optimality condition is obtained as follows

{\begin{matrix} 2 h (x_{1} (τ) - {\bar{x}}_{1})^{2} - λ_{1} K (A {(\frac{τ}{1 - τ})}^{- α} + 1) - λ_{2} K x_{3} (τ) \\ + 4 λ_{3} K x_{2} (τ) + λ_{3} K x_{2} (τ) - h ϒ_{1} - \frac{2 τ}{1 - τ} ϒ_{4} = - {\overset{\cdot}{λ}}_{1}, \\ 2 h (x_{2} (τ) - {\bar{x}}_{2})^{2} + 2 λ_{1} K x_{2} (τ) - λ_{2} K (A (\frac{τ}{1 - τ} - 2.5)^{- α}) \\ + 4 λ_{3} K x_{1} (τ) - h ϒ_{2} - \frac{2 τ}{1 - τ} ϒ_{5} = - {\overset{\cdot}{λ}}_{2}, \\ 2 h (x_{3} (τ) - {\bar{x}}_{3})^{2} - 4 λ_{2} K x_{1} (τ) - λ_{3} K (A {(\frac{τ}{1 - τ})}^{- α} + 5) - h ϒ_{3} \\ - \frac{2 τ}{1 - τ} ϒ_{6} = - {\overset{\cdot}{λ}}_{3}, \\ h [C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) \\ - (A {(\frac{τ}{1 - τ})}^{- α} + 1) x_{1} (τ) - x_{2}^{2} (τ) + u_{1} (τ)] = {\overset{\cdot}{x}}_{1}, \\ h [C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) + (A {(\frac{τ}{1 - τ})}^{- α} \\ - 2.5) x_{2} (τ) - 4 x_{1} (τ) x_{3} (τ) + u_{2} (τ)] = {\overset{\cdot}{x}}_{2}, \\ h [C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) - (A {(\frac{τ}{1 - τ})}^{- α} + 5) x_{3} \\ + 4 x_{1} (τ) x_{2} (τ) + u_{3} (τ)] = {\overset{\cdot}{x}}_{3}, \\ λ_{1} K C_{2} ({\frac{τ}{1 - τ}}^{- 1 - α}) = - {\overset{\cdot}{ϒ}}_{1}, λ_{2} K C_{2} ({\frac{τ}{1 - τ}}^{- 1 - α}) \\ = - {\overset{\cdot}{ϒ}}_{2}, λ_{3} K C_{2} ({\frac{τ}{1 - τ}}^{- 1 - α}) = - {\overset{\cdot}{ϒ}}_{3}, \\ λ_{1} K C_{3} ({\frac{τ}{1 - τ}}^{- 1 - α}) = - {\overset{\cdot}{ϒ}}_{4}, λ_{2} K C_{3} ({\frac{τ}{1 - τ}}^{- 1 - α}) \\ = - {\overset{\cdot}{ϒ}}_{5}, λ_{3} K C_{3} ({\frac{τ}{1 - τ}}^{- 1 - α}) = - {\overset{\cdot}{ϒ}}_{6}, \\ - h x_{1} = {\overset{\cdot}{y}}_{21}, - h x_{1} = {\overset{\cdot}{y}}_{22}, - h x_{1} = {\overset{\cdot}{y}}_{23}, \\ {\frac{- 2 τ}{1 - τ}}^{3} = {\overset{\cdot}{y}}_{31}, {\frac{- 2 τ}{1 - τ}}^{3} = {\overset{\cdot}{y}}_{32}, {\frac{- 2 τ}{1 - τ}}^{3} = {\overset{\cdot}{y}}_{33}, \\ \frac{λ_{1}}{{(1 - τ)}^{2}} = 0, \frac{λ_{2}}{{(1 - τ)}^{2}} = 0, \frac{λ_{3}}{{(1 - τ)}^{2}} = 0 . \end{matrix}

According to the initial conditions, trial solutions may be selected as

{\begin{matrix} \begin{matrix} x_{1_{T}} = 0.2 + τ N_{x_{1}} & x_{2_{T}} = τ N_{x_{2}} & x_{3_{T}} = 0.5 + τ N_{x_{3}} \\ y_{21_{T}} = τ N_{y_{21}}, & y_{22_{T}} = τ N_{y_{22}}, & y_{23_{T}} = τ N_{y_{23}}, \\ y_{31_{T}} = τ N_{y_{31}}, & y_{32_{T}} = τ N_{y_{32}}, & y_{33_{T}} = τ N_{y_{33}}, \\ λ_{1_{T}} = (1 - τ) N_{λ_{1}}, & λ_{2_{T}} = (1 - τ) N_{λ_{2}}, & λ_{3_{T}} = (1 - τ) N_{λ_{3}}, \\ ϒ_{1_{T}} = (1 - τ) N_{ϒ_{1}}, & ϒ_{2_{T}} = (1 - τ) N_{ϒ_{2}}, & ϒ_{3_{T}} = (1 - τ) N_{ϒ_{3}}, \\ u_{1_{T}} = N_{u_{1}}, & u_{2_{T}} = N_{u_{2}}, & u_{3_{T}} = N_{u_{3}} . \end{matrix} \end{matrix}

It can be seen the approximate solutions for the state and control functions for $α = 0.9, 0.8$ and $α = 0.5$ in Figures 8 –10, respectively.

Figure 8.

The stabilized behaviors of state and control for $α = 0.9$ .

Figure 9.

The stabilized behaviors of state and control for $α = 0.8$ .

Figure 10.

The stabilized behaviors of state and control for $α = 0.5$ .

Example 4: Consider the following system

{\begin{matrix} {{}_{0}^{C}D}_{t}^{α} X_{1} (t) = 35 (X_{2} (t) - X_{1} (t)), \\ {{}_{0}^{C}D}_{t}^{α} X_{1} (t) = - 7 X_{1} (t) - X_{1} (t) X_{3} (t) + 28 X_{2} (t), \\ {{}_{0}^{C}D}_{t}^{α} X_{1} (t) = X_{1} (t) X_{2} (t) - 3 X_{3} (t), \\ X_{1} (0) = - 9, X_{2} (0) = - 5, X_{3} (0) = 14 . \end{matrix}

(47)

By choosing controls $U_{1} (t), U_{2} (t), U_{3} (t)$ , an equivalent optimization problem is defined as

minimize J = \int_{0}^{\infty} [(X_{1} (t) - {\bar{X}}_{1})^{2} + (X_{2} (t) - {\bar{X}}_{2})^{2} + (X_{3} (t) - {\bar{X}}_{3})^{2}] dt,

(48)

subject to {\begin{matrix} {{}_{0}^{C}D}_{t}^{α} X_{1} (t) = 35 (X_{2} (t) - X_{1} (t)) + U_{1} (t), \\ {{}_{0}^{C}D}_{t}^{α} X_{2} (t) = - 7 X_{1} (t) - X_{1} (t) X_{3} (t) + 28 X_{2} (t) + U_{2} (t), \\ {{}_{0}^{C}D}_{t}^{α} X_{3} (t) = X_{1} (t) X_{2} (t) - 3 X_{3} (t) + U_{3} (t), \\ X_{1} (0) = - 9, X_{2} (0) = - 5, X_{3} (0) = 14 . \end{matrix}

(49)

Now by the change of variable (46), approximation (18) and assuming

{\begin{matrix} x_{i} (τ) = X_{i} (\frac{τ}{1 - τ}), & i = 1, 2, 3, \\ u_{i} (τ) = U_{i} (\frac{τ}{1 - τ}), & i = 1, 2, 3, \\ y_{i} (τ) = Y_{i} (\frac{τ}{1 - τ}), & i = 1, 2, 3 . \end{matrix}

The problem (48) and (49) for $p = 2, 3$ is reduced to the following problem

\begin{matrix} minimize J = \int_{0}^{1} \frac{1}{{(1 - τ)}^{2}} \\ [(x_{1} (τ) - {\bar{X}}_{1})^{2} + (x_{2} (τ) - {\bar{X}}_{2})^{2} + (x_{3} (τ) - {\bar{X}}_{3})^{2}] d τ \end{matrix}

subject to {\begin{matrix} {\overset{\cdot}{x}}_{1} = \frac{1}{{(1 - τ)}^{2} B {(\frac{τ}{1 - τ})}^{1 - α}} (C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) \\ - (A {(\frac{τ}{1 - τ})}^{- α} + 35) x_{1} (τ) + 35 x_{2} (τ) + u_{1} (τ)), \\ {\overset{\cdot}{x}}_{2} = \frac{1}{{(1 - τ)}^{2} B {(\frac{τ}{1 - τ})}^{1 - α}} (C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) \\ + (28 - A {(\frac{τ}{1 - τ})}^{- α}) x_{2} (τ) - 7 x_{1} (τ) - x_{1} (τ) x_{3} (τ) + u_{2} (τ)), \\ {\overset{\cdot}{x}}_{3} = \frac{1}{{(1 - τ)}^{2} B {(\frac{τ}{1 - τ})}^{1 - α}} (C_{2} {(\frac{τ}{1 - τ})}^{- 1 - α} y_{21 (τ)} + C_{3} {(\frac{τ}{1 - τ})}^{- 2 - α} y_{31} (τ) \\ - (A {(\frac{τ}{1 - τ})}^{- α} + 3) x_{3} + x_{1} (τ) x_{2} (τ) + u_{3} (τ)), \\ {\overset{\cdot}{y}}_{21} (t) = - \frac{1}{{(1 - τ)}^{2}} x_{1} (τ), \\ {\overset{\cdot}{y}}_{22} (t) = - \frac{1}{{(1 - τ)}^{2}} x_{2} (τ), \\ {\overset{\cdot}{y}}_{23} (t) = - \frac{1}{{(1 - τ)}^{2}} x_{3} (τ), \\ {\overset{\cdot}{y}}_{31} (t) = - \frac{2 τ}{{(1 - τ)}^{3}} x_{1} (τ), \\ {\overset{\cdot}{y}}_{32} (t) = - \frac{2 τ}{{(1 - τ)}^{3}} x_{1} (τ), \\ {\overset{\cdot}{y}}_{33} (t) = - \frac{2 τ}{{(1 - τ)}^{3}} x_{1} (τ), \\ x_{1} (0) = - 9, x_{2} (0) = - 5, x_{3} (0) = 14, \\ y_{21} (0) = 0, y_{22} (0) = 0, y_{23} (0) = 0, \\ y_{31} (0) = 0, y_{32} (0) = 0, y_{33} (0) = 0 . \end{matrix}

According to the initial conditions, the trial solutions may be selected as

{\begin{matrix} x_{1_{T}} = - 9 + τ N_{x_{1}}, x_{2_{T}} = - 5 + τ N_{x_{2}}, x_{3_{T}} = 14 + τ N_{x_{3}}, \\ y_{21_{T}} = τ N_{y_{21}}, y_{22_{T}} = τ N_{y_{22}}, y_{23_{T}} = τ N_{y_{23}}, \\ y_{31_{T}} = τ N_{y_{31}}, y_{32_{T}} = τ N_{y_{32}}, y_{33_{T}} = τ N_{y_{33}}, \\ λ_{1_{T}} = (1 - τ) N_{λ_{1}}, λ_{2_{T}} = (1 - τ) N_{λ_{2}}, λ_{3_{T}} = (1 - τ) N_{λ_{3}}, \\ ϒ_{1_{T}} = (1 - τ) N_{ϒ_{1}}, ϒ_{2_{T}} = (1 - τ) N_{ϒ_{2}}, ϒ_{3_{T}} = (1 - τ) N_{ϒ_{3}}, \\ u_{1_{T}} = N_{u_{1}}, u_{2_{T}} = N_{u_{2}}, u_{3_{T}} = N_{u_{3}} . \end{matrix}

(50)

The approximate solutions for state and control functions for $α = 0.9, 0.8$ and $α = 0.5$ can be seen in Figures 11 –13, respectively.

Figure 11.

The stabilized behaviors of state and control for $α = 0.9$ .

Figure 12.

The stabilized behaviors of state and control for $α = 0.8$ .

Figure 13.

The stabilized behaviors of state and control for $α = 0.5$ .

Example 5: Akinlar et al. (2020) Consider the fractional-order Malkus waterwheel model in the form of

{\begin{matrix} {{}_{0}^{c}D}_{t}^{α} (X_{1} (t)) = X_{2} (t) - X_{1} (t), \\ {{}_{0}^{c}D}_{t}^{α} (X_{2} (t)) = X_{1} (t) X_{3} (t) - a X_{2} (t), \\ {{}_{0}^{c}D}_{t}^{α} (X_{3} (t)) = - X_{1} (t) X_{2} (t) - b X_{3} (t) + b, \end{matrix}

(51)

where $X_{1} (t), X_{2} (t), X_{3} (t)$ are spatial coordinates and $a > 0$ and $b > 0$ are constant parameters. The initial conditions of the system are given by ${[X_{1} (0), X_{2} (0), X_{3} (0)]}^{T} = {[0.51, 0.60407, 0.07]}^{T}$ . $p^{*} = (0, 0, 1)$ is considered as the equilibrium point of the fractional order system (51). In this case, the input controls $U_{1} (t), U_{2} (t)$ and $U_{3} (t)$ are added to the all state equation in order to control chaos in the fractional-order Malkus waterwheel system. Figure 14 presents the approximate solutions of the states and controls for $α = 0.98$ with $I = 5, m = 10$ and $k = 2$ . The values of the error function $E (P)$ is $1.300131 \times 10^{- 5}$ , which shows that the approximate solutions satisfy the necessary condition of optimality.

Figure 14.

The stabilized behaviors of state and control for $α = 0.98$ .

Example 6: Xin and Zhang (2015) Consider the fractional-order financial model as follows

{\begin{matrix} {{}_{0}^{c}D}_{t}^{α} (X_{1} (t)) = X_{3} (t) + (X_{2} (t) - a) X_{1} (t) + m_{1} X_{4} (t), \\ {{}_{0}^{c}D}_{t}^{α} (X_{2} (t)) = 1 - b X_{2} (t) - X_{1}^{2} (t) + m_{2} X_{4} (t), \\ {{}_{0}^{c}D}_{t}^{α} (X_{3} (t)) = - X_{1} (t) - c X_{3} (t) + m_{3} X_{4} (t), \\ {{}_{0}^{c}D}_{t}^{α} (X_{3} (t)) = - X_{1} (t) X_{2} (t) X_{3} (t), \end{matrix}

(52)

where $X_{1}$ denotes the interest rate, $X_{2}$ denotes the investment demand, $X_{3}$ denotes the price index, $X_{4}$ denotes the market confidence, $a \geq 0$ is the saving amount, $b \geq 0$ is the cost per investment, $c \geq 0$ is the demand elasticity of commercial markets and $m_{1}, m_{2}, m_{3}$ are the impact factors. We fix $α = 0.95, a = 2.1, b = 0.01, c = 2.6, m_{1} = 8.4, m_{2} = 6.4, m_{3} = 2.2,$ initial value $[X_{1} (0), X_{2} (0), X_{3} (0), X_{4} (0)]^{T} = [1, 5, 4, 3]^{T}$ and equilibrium point $p^{*} = (0, 0, 0, 0)$ . In this case, the input controls $U_{1} (t), U_{2} (t), U_{3} (t)$ and $U_{4} (t)$ are added to the all state equation in order to control chaos in the fractional-order financial system. In Figure 15, the approximate solution of the states and controls with $I = 5, m = 10$ and $k = 2$ are shown. The values of the error function $E (P)$ is $2.74561 \times 10^{- 2}$ , which shows that the approximate solutions satisfy the necessary condition of optimality.

Figure 15.

The stabilized behaviors of state and control for $α = 0.95$ .

Example 7: Al-Khedhairi et al. (2019) Consider the fractional-order geomagnetic field model in the form of

{\begin{matrix} {{}_{0}^{c}D}_{t}^{α} (X_{1} (t)) = a X_{1} (t) - X_{2} (t) X_{3} (t), \\ {{}_{0}^{c}D}_{t}^{α} (X_{2} (t)) = - b X_{2} (t) + X_{1} (t) X_{3} (t), \\ {{}_{0}^{c}D}_{t}^{α} (X_{3} (t)) = c - X_{3} (t) + X_{1} (t) X_{2} (t), \end{matrix}

(53)

where $a, b, c$ are all real valued parameters. The initial conditions of the system are given by $[X_{1} (0), X_{2} (0), X_{3} (0)]^{T} = [1, 1, 1]^{T}$ . $p^{*} = (0, 0, 1)$ is considered as the equilibrium point of the fractional order system (53). In this case, the input controls $U_{1} (t), U_{2} (t)$ and $U_{3} (t)$ are added to the all state equation in order to control chaos in the fractional-order geomagnetic field system. Figure 16 presents the approximate solutions of the states and controls for $α = 0.95$ with $I = 5, m = 15$ and $k = 2$ . The values of the error function $E (P)$ is $5.014544 \times 10^{- 4}$ , which shows that the approximate solutions satisfy the necessary condition of optimality.

Figure 16.

The stabilized behaviors of state and control for $a = 3, b = 5, c = 1, α = 0.95$ .

To end this section, some advantages of the proposed method and the theoretical results are summarized as below:

For approximation of fractional derivatives, the approximation that was carried out in Almeida and Torres (2015) and Pooseh et al. (2013) is used, where a good approximation is obtained without the requirement of such higher-order smoothness on the admissible functions.

Some numerical methods use the operational matrices for solving FOCPs (see Heydari et al., 2016; Hosseinpour and Nazemi, 2015; Jafari and Tajadodi, 2014). However, finding the operational matrices is usually difficult. Moreover, the operational matrices with high dimensions cause the complexity in computations. But, in the proposed neural network scheme, the role of the operational matrices has been removed. Thus, the computational complexity has significantly reduced.

By providing an error table for different values of $α$ , it is confirmed that the solutions for $0 < α < 1$ are accurate.

The solution of the problem is continuous over all the domain of the problem. In contrast, the numerical methods provide solutions only over discrete points; and the solution between these points must be interpolated.

The rounding-off error propagation of standard numerical methods does not affect the neural network solution.

A simple and reliable learning algorithm, we use a gradient descent method with good convergence properties compared with some other optimization algorithms. The convergence analysis of the proposed algorithm is also provided.

Conclusion

In this article that has been studied, a good approximation is applied. It is obtained without the requirement of such higher-order smoothness on the admissible functions. By transforming a fractional problem into a classical one, one can benefit from the vast number of available techniques in the field. An approach based on feed-forward neural networks, optimization techniques and collocation methods is then proposed to determine the approximate solution of the OCPs. Neural networks, due to lower computational complexity have more benefits than other methods for solving the fractional differential equations and related optimal control problems. The most interesting characteristic of the neural networks is their capability of the formulation of problems by means of training. After sufficient training, the neural networks can solve problems of the same class, since training algorithms converge to the optimal solutions.

As shown in numerical results, unlike some other methods that use operational matrices with high computational complexity, in the proposed method with choosing the trial solutions and putting them in optimality condition, an unconstrained optimization problem is easily solved. Moreover, without the requirement of higher-order smoothness on the functions, with the good approximation in (14), the FOCP problem is replaced with the OCP problem without fractional derivative order. This, in fact, helps to convert the fractional problem to a non-fractional one, which is very effective in reducing computational complexity of the neural network scheme. These can be noticed as some innovations of the proposed scheme. Another novelty as a result of the figures is that the stabilized behaviors of state and control for different values of $α$ are plotted. It also provided the values of error function $E (P)$ for several values $0 < α < 1$ , which proves the solutions for $0 < α < 1$ are accurate. For future works, the proposed method can be utilized for synchronization and stabilization of various fractional chaotic systems.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Alireza Nazemi

References

Akinlar

Tchier

Inc

(2020) Chaos control and solutions of fractional-order malkus waterwheel model. Chaos, Solitons & Fractals 135: 109746.

Al-Khedhairi

Matouk

Khan

(2019) Chaotic dynamics and chaos control for the fractional-order geomagnetic field model. Chaos, Solitons and Fractals 128: 390–401.

Alizadeh

Effati

(2018) An iterative approach for solving fractional optimal control problems. Journal of Vibration and Control 24(1): 18–36.

Almeida

Torres

(2015) A discrete method to solve fractional optimal control problems. Nonlinear Dynamics 80(4): 1811–1816.

Ávalos-Ruiz

Zúñiga-Aguilar

Gómez-Aguilar

, et al. (2018) FPGA implementation and control of chaotic systems involving the variable-order fractional operator with Mittag–Leffler law. Chaos, Solitons & Fractals 115: 177–189.

Baskonus

Mekkaoui

Hammouch

Bulut

(2015) Active control of a chaotic fractional order economic system. Entropy 17(8): 5771–5783.

Baskonus

Bulut

(2015) On the numerical solutions of some fractional ordinary differential equations by fractional adams-bashforth-moulton method. Open Mathematics 13(1): 547–556.

Baskonus

Bulut

(2016) Regarding on the prototype solutions for the nonlinear fractional-order biological population model. In: AIP Conference Proceedings, volume 1738, pp. 2900041-2900046. College Park, Maryland: AIP Publishing

Baskonus

Yel

Bulut

(2017) Novel wave surfaces to the fractional zakharov-kuznetsov-benjamin-bona-mahony equation. In: AIP Conference Proceedings, 21 July 2017, volume 1863, pp. 5600841-5600845. College Park, Maryland: AIP Publishing.

10.

Bazaraa

Sherali

Shetty

(2013) Nonlinear Programming: Theory and Algorithms. New Jersey, USA: John Wiley & Sons.

11.

Beidokhti

Malek

(2009) Solving initial-boundary value problems for systems of partial differential equations using neural networks and optimization techniques. Journal of the Franklin Institute 346(9): 898–913.

12.

Bhrawy

Ezz-Eldien

Doha

, et al. (2017) Solving fractional optimal control problems within a chebyshev–legendre operational technique. International Journal of Control 90(6): 1230–1244.

13.

Bigdeli

Ziazi

(2017) Finite-time fractional-order adaptive intelligent backstepping sliding mode control of uncertain fractional-order chaotic systems. Journal of the Franklin Institute 354(1): 160–183.

14.

Bulut

Kumar

Singh

, et al. (2018a) Analytic study for a fractional model of hiv infection of cd4+ tcd4+ t lymphocyte cells. Mathematics in Natural Science 2(1): 33–43.

15.

Bulut

Sulaiman

Baskonus

(2018b) Dark, bright optical and other solitons with conformable space-time fractional second-order spatiotemporal dispersion. Optik 163: 1–7.

16.

Cheng

Lewis

Abu-Khalaf

(2007) Fixed-final-time-constrained optimal control of nonlinear systems using neural network hjb approach. IEEE Transactions on Neural Networks 18(6): 1725–1737.

17.

Daniel

(2013) Principles of Artificial Neural Networks, volume 7. Singapore: World Scientific.

18.

Das

Roy

Das

(2000) Chaos in a three dimensional neural network. Applied Mathematical Modelling 24(7): 511–522.

19.

Dehghani

Khanlo

(2019) Radial basis function neural network chaos control of a piezomagnetoelastic energy harvesting system. Journal of Vibration and Control 25(16): 2191–2203.

20.

Dehghan

Hamedi

Khosravian-Arab

(2016) A numerical scheme for the solution of a class of fractional variational and optimal control problems using the modified Jacobi polynomials. Journal of Vibration and Control 22(6): 1547–1559.

21.

Dokuyucu

Celik

Bulut

Baskonus

(2018) Cancer treatment model with the caputo-fabrizio fractional derivative. The European Physical Journal Plus 133(3): 92: 1–6.

22.

Dua

(2011) An artificial neural network approximation based decomposition approach for parameter estimation of system of ordinary differential equations. Computers & Chemical Engineering 35(3): 545–553.

23.

Effati

Pakdaman

(2013) Optimal control problem via neural networks. Neural Computing and Applications 23(7–8): 2093–2100.

24.

Ejlali

Hosseini

(2017) A pseudospectral method for fractional optimal control problems. Journal of Optimization Theory and Applications 174(1): 83–107.

25.

Esen

Sulaiman

Bulut

Baskonus

(2018) Optical solitons to the space-time fractional (1+ 1)-dimensional coupled nonlinear schrödinger equation. Optik 167: 150–156.

26.

Eshaghi

Ghaziani

Ansari

(2020) Hopf bifurcation, chaos control and synchronization of a chaotic fractional-order system with chaos entanglement function. Mathematics and Computers in Simulation 172: 321–340.

27.

Gencoglu

Baskonus

Bulut

(2017) Numerical simulations to the nonlinear model of interpersonal relationships with time fractional derivative. In: AIP Conference Proceedings, 27 January 2017, volume 1798, p. 0201031-0201039. College Park, Maryland: AIP Publishing.

28.

Ghasemi

Nazemi

(2018) A neural network method based on mittag-leffler function for solving a class of fractional optimal control problems. AUT Journal of Modeling and Simulation 50(2): 211–218.

29.

Ghasemi

Nazemi

Hosseinpour

(2017) Nonlinear fractional optimal control problems with neural network and dynamic optimization schemes. Nonlinear Dynamics 89(4): 2669–2682.

30.

Haykin

(1994) Neural Networks: A Comprehensive Foundation. New Jersey, United States: Prentice Hall PTR.

31.

Heydari

Hooshmandasl

Ghaini

Cattani

(2016) Wavelets method for solving fractional optimal control problems. Applied Mathematics and Computation 286: 139–154.

32.

Hosseinpour

Nazemi

(2015) Solving fractional optimal control problems with fixed or free final states by haar wavelet collocation method. IMA Journal of Mathematical Control and Information 33(2): 543–561.

33.

Hosseinpour

Nazemi

(2016) A collocation method via block-pulse functions for solving delay fractional optimal control problems. IMA Journal of Mathematical Control and Information 34(4): 1215–1237.

34.

Jafari

Tajadodi

(2014) Fractional order optimal control problems via the operational matrices of bernstein polynomials. UPB Scientific Bulletin, Series A: Applied Mathematics and Physics 76(3): 115–128.

35.

Jafarian

Mokhtarpour

Baleanu

(2017) Artificial neural network approach for a class of fractional ordinary differential equation. Neural Computing and Applications 28(4): 765–773.

36.

Kheyrinataj

Nazemi

(2020a) Fractional Chebyshev functional link neural network-optimization method for solving delay fractional optimal control problems with atangana-baleanu derivative. Optimal Control Applications and Methods 41(3): 808–832.

37.

Kheyrinataj

Nazemi

(2020b) Fractional power series neural network for solving delay fractional optimal control problems. Connection Science 32(1): 53–80.

38.

Kilbas

AAA

Srivastava

Trujillo

(2006) Theory and Applications Of Fractional Differential Equations, volume 204. Netherlands: Elsevier Science Limited.

39.

Kirk

(2012) Optimal Control Theory: An Introduction. Mineola, New York: Courier Corporation.

40.

Kumar

Yadav

(2011) Multilayer perceptrons and radial basis function neural network methods for the solution of differential equations: a survey. Computers & Mathematics with Applications 62(10): 3796–3811.

41.

Kumar

Singh

Baskonus

Bulut

(2017) An effective computational approach for solving local fractional telegraph equations. Nonlinear Science Letters A: Mathematics, Physics and Mechanics 8(2): 200–206.

42.

Lee

El-Sharkawi

(2008) Modern Heuristic Optimization Techniques: Theory and Applications to Power Systems, volume 39. New Jersey, USA: John Wiley & Sons.

43.

Lewis

Syrmos

(1987) Optimal Control, volume 31. Atlanta: School of Electrical Engineering, Georgia Institute of Technology.

44.

Sun

(2019) Adaptive neural network backstepping control of fractional-order chua–hartley chaotic system. Advances in Difference Equations 2019(1): 148: 1–14.

45.

Lotfi

Yousefi

Dehghan

(2013) Numerical solution of a class of fractional optimal control problems via the legendre orthonormal basis combined with the operational matrix and the gauss quadrature rule. Journal of Computational and Applied Mathematics 250: 143–160.

46.

Wang

(2019) Observer-based command filtered adaptive neural network tracking control for fractional-order chaotic pmsm. IEEE Access 7: 88777–88788.

47.

Machado

Kiryakova

Mainardi

(2011) Recent history of fractional calculus. Communications in Nonlinear Science and Numerical Simulation 16(3): 1140–1153.

48.

Magin

(2006) Fractional Calculus in Bioengineering. USA: Begell House Redding.

49.

Mahmoud

Arafa

Abed-Elhameed

Mahmoud

(2017) Chaos control of integer and fractional orders of chaotic burke–shaw system using time delayed feedback control. Chaos, Solitons & Fractals 104: 680–692.

50.

Mall

Chakraverty

(2016) Application of legendre neural network for solving ordinary differential equations. Applied Soft Computing 43: 347–356.

51.

Mei

Bullo

(2017) Lasalle invariance principle for discrete-time dynamical systems: a concise and self-contained tutorial. Available at: arXiv:1710.03710v5

52.

Müller

Reinhardt

Strickland

(2012) Neural Networks: An Introduction. Heidelberg, Germany: Springer Science & Business Media.

53.

Nazemi

(2011) A dynamical model for solving degenerate quadratic minimax problems with constraints. Journal of Computational and Applied Mathematics 236(6): 1282–1295.

54.

Nazemi

(2012) A dynamic system model for solving convex nonlinear optimization problems. Communications in Nonlinear Science and Numerical Simulation 17(4): 1696–1705.

55.

Nazemi

(2013) Solving general convex nonlinear optimization problems by an efficient neurodynamic model. Engineering Applications of Artificial Intelligence 26(2): 685–696.

56.

Nazemi

(2014) A neural network model for solving convex quadratic programming problems with some applications. Engineering Applications of Artificial Intelligence 32: 54–62.

57.

Nazemi

Effati

(2013) An application of a merit function for solving convex programming problems. Computers & Industrial Engineering 66(2): 212–221.

58.

Nazemi

Karami

(2017) A neural network approach for solving optimal control problems with inequality constraints and some applications. Neural Processing Letters 45(3): 995–1023.

59.

Nazemi

(2014) A gradient-based neural network method for solving strictly convex quadratic programming problems. Cognitive Computation 6(3): 484–495.

60.

Nazemi

Omidi

(2012) A capable neural network model for solving the maximum flow problem. Journal of Computational and Applied Mathematics 236(14): 3498–3513.

61.

Nazemi

Omidi

(2013) An efficient dynamic model for solving the shortest path problem. Transportation Research Part C: Emerging Technologies 26: 1–19.

62.

Nazemi

Sharifi

(2013) Solving a class of geometric programming problems by an efficient dynamic model. Communications in Nonlinear Science and Numerical Simulation 18(3): 692–709.

63.

Nocedal

Wright

(2006) Numerical Optimization. Heidelberg, Germany: Springer Science & Business Media.

64.

Ott

Grebogi

Yorke

(1990) Controlling chaos. Physical Review Letters 64(11): 1196–1199.

65.

Oustaloup

Levron

Mathieu

Nanot

(2000) Frequency-band complex noninteger differentiator: characterization and synthesis. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications 47(1): 25–39.

66.

Petráš

(2011) Fractional-order Nonlinear Systems: Modeling, Analysis and Simulation. Heidelberg, Germany: Springer Science & Business Media.

67.

Podlubny

(1998). Fractional differential equations. San Diego, USA: Academic Press.

68.

Pooseh

Almeida

Torres

(2013) Numerical approximations of fractional derivatives with applications. Asian Journal of Control 15(3): 698–712.

69.

Prakasha

Veeresha

Baskonus

(2019) Residual power series method for fractional swift–hohenberg equation. Fractal and Fractional 3(1): 1–16.

70.

Rabiei

Ordokhani

Babolian

(2017) The boubaker polynomials and their application to solve fractional optimal control problems. Nonlinear Dynamics 88(2): 1013–1026.

71.

Rakhshan

Effati

Vahidian Kamyad

(2018) Solving a class of fractional optimal control problems by the hamilton–jacobi–bellman equation. Journal of Vibration and Control 24(9): 1741–1756.

72.

Ravichandran

Jothimani

Baskonus

Valliammal

(2018) New results on nondensely characterized integrodifferential equations with fractional order. The European Physical Journal Plus 133(3): 109.

73.

Rong

Xiaoning

(1998) From Chaos to Order: Methodologies, Perspectives and Applications, volume 24. Singapore: World Scientific.

74.

Samko

Kilbas

Marichev

, et al. (1993) Fractional Integrals and Derivatives, volume 1993. Switzerland: Gordon and Breach Science Publishers, Yverdon Yverdon-les-Bains.

75.

Shirvany

Hayati

Moradian

(2009) Multilayer perceptron neural networks with novel unsupervised training method for numerical solution of the partial differential equations. Applied Soft Computing 9(1): 20–29.

76.

Singha

Nahak

(2017) An efficient approximation technique for solving a class of fractional optimal control problems. Journal of Optimization Theory and Applications 174(3): 785–802.

77.

Sulaiman

Baskonus

Bulut

(2018) Optical solitons and other solutions to the conformable space–time fractional complex ginzburg–landau equation under Kerr law nonlinearity. Pramana 91(4): 1–8.

78.

Sweilam

Al-Ajami

(2015) Legendre spectral-collocation method for solving some types of fractional optimal control problems. Journal of Advanced Research 6(3): 393–403.

79.

Tang

Tan

(2007) Neural Networks: Computational Models and Applications, volume 53. Heidelberg, Germany: Springer Science & Business Media.

80.

Tang

Shi

Wang

(2017) A new framework for solving fractional optimal control problems using fractional pseudospectral methods. Automatica 78: 333–340.

81.

Tohidi

Nik

(2015) A Bessel collocation method for solving fractional optimal control problems. Applied Mathematical Modelling 39(2): 455–465.

82.

Torvik

Bagley

(1984) On the Appearance of the Fractional Derivative in the Behavior of Real Materials. Journal of Applied Mechanics 51(2): 294–298.

83.

Veeresha

Prakasha

Baskonus

(2019) New numerical surfaces to the mathematical model of cancer chemotherapy effect in caputo fractional derivatives. Chaos: An Interdisciplinary Journal of Nonlinear Science 29(1): 0131191-01311914.

84.

Vrabie

Lewis

(2009) Neural network approach to continuous-time direct adaptive optimal control for partially unknown nonlinear systems. Neural Networks 22(3): 237–246.

85.

Xin

Zhang

(2015) Finite-time stabilizing a fractional-order chaotic financial system with market confidence. Nonlinear Dynamics 79(2): 1399–1409.

86.

Yadav

Das

Bhadauria

, et al. (2017) Stability analysis, chaos control of a fractional order chaotic chemical reactor system and its function projective synchronization with parametric uncertainties. Chinese Journal of Physics 55(3): 594–605.

87.

Yavari

Nazemi

(2019) An efficient numerical scheme for solving fractional infinite-horizon optimal control problems. ISA Transactions 94: 108–118.

88.

Yavari

Nazemi

(2020) On fractional infinite-horizon optimal control problems with a combination of conformable and caputo–fabrizio fractional derivatives. ISA Transactions 101: 78–90.

89.

Yavuz

Ozdemir

Baskonus

(2018) Solutions of partial differential equations using the fractional operator involving Mittag-Leffler kernel. The European Physical Journal Plus 133(6): 1–11.

90.

Yong

Jinrong

(2016) Basic Theory of Fractional Differential Equations. Singapore: World Scientific.

91.

, et al. (1998) Introduction to Control Of Oscillations and Chaos, volume 35. Singapore: World Scientific.

92.

Zamani

Karimi-Ghartemani

Sadati

(2007) Fopid controller design for robust performance using particle swarm optimization. Fractional Calculus and Applied Analysis 10(2): 169–187.

93.

Zhang

(2013) Neural Networks in Optimization, volume 46. Heidelberg, Germany: Springer Science & Business Media.

94.

Zuniga-Aguilar

Coronel-Escamilla

Gómez-Aguilar

, et al. (2018) New numerical approximation for solving fractional delay differential equations of variable order using artificial neural networks. The European Physical Journal Plus 133(2): 1–16.