Optimal Control of Time-Delay Fractional Equations via a Joint Application of Radial Basis Functions and Collocation Method

A novel approach to solve optimal control problems dealing simultaneously with fractional differential equations and time delay is proposed in this work. More precisely, a set of global radial basis functions are firstly used to approximate the states and control variables in the problem. Then, a collocation method is applied to convert the time-delay fractional optimal control problem to a nonlinear programming one. By solving the resulting challenge, the unknown coefficients of the original one will be finally obtained. In this way, the proposed strategy introduces a very tunable framework for direct trajectory optimization, according to the discretization procedure and the range of arbitrary nodes. The algorithm’s performance has been analyzed for several non-trivial examples, and the obtained results have shown that this scheme is more accurate, robust, and efficient than most previous methods.


Introduction
In the last years, the use of fractional calculus has increased significantly due to its attractive applications in physical and engineering systems [1][2][3], materials [4], biology [5], finance [6], and so on. Moreover, fractional differential equations (FDEs) have also recently gained considerable importance in pure and applied mathematics [7], engineering [8], physics [9], and bio-systems [10]. Nonetheless, despite this growing variety of applications, it is often difficult to find numerical methods with low computing cost and enough accuracy for resolving these kinds of equations and analytically framework for solving DFOCPs. The practical importance of the proposed method is that a variety of RBF functions can be applied for interpolation of states and controls, instead of being limited to a specific type of polynomial as in polynomial-based methods. Moreover, a wide range of arbitrary nodes can also be easily employed for discretization of the fractional terms, thus resulting in a flexible RBF framework for solving DFOCPs.
The outline of this paper is as follows. Section 2 demonstrates the problem statement and the basic concepts about fractional derivative. Some preliminaries of RBFs for subsequent developments are presented in Section 3. Moreover, we present a direct RBF collocation scheme to solve DFOCPs in this section. The numerical results obtained by the proposed approach for some non-trivial examples are described and compared with other previous works in Section 4. Finally, the most relevant conclusions are summarized in Section 5, along with some future perspectives.

Statement of the Problem
The aforementioned performance of meshless methods have encouraged some researchers to develop new computing architectures and techniques where the primary focus was on hardware simplicity. In order to lower the implementation cost, we want to explore an applicable numerical scheme to find the approximate solutions of the following DFOCP, subject to dynamic constraints, where and u(t) ∈ U mentions the control variable, in which U ⊂ R m represents the set of continuous functions. Furthermore, it is assumed that In addition, a(t), b(t), e(t), f (t), and g(t) are continuous functions; ϕ 1 (t) and ϕ 2 (t) are known functions; and r(t) and q(t) are two symmetric positive semidefinite and definite matrixes, respectively, which show the time-varying coefficients of the state and control variables in the cost function with continuous functions. Moreover, it is assumed that the dynamic system (2) is at rest from −∞ to t 0 − η. Furthermore, D α is the fractional differentiation operator of order α that is defined as follows.

Definition 1.
For a given function f (t) and α > 0, n − 1 < α ≤ n, n ∈ N, the operators Furthermore, The aforementioned properties of CFDs have led us to use this definition in the following. The main contribution of this paper is thus to suggest a direct method based on RBFs and collocation points to obtain the optimal values of u(t) and x(t), t ∈ [t 0 , t f ], satisfying Equation (2) and minimizing the quadratic performance index in Equation (1). One advantage of this method is that it does not use the maximum principle and calculate pontryagin variations to solve the problem, so there is no need for analytical separation of cost and constraint statements. Moreover, in general terms, the direct methods (such as the proposed one) have a greater convergence radius than indirect methods [33,34]. Moreover, to make the problem significantly simpler, we have tried to reformulated the DFOCP expressed in Equations (1) and (2) as an equivalent NLP by making use of the interpolate approximate of basis functions.

Method of Solution
In this section, a brief description of the proposed method to directly solve the DFOCP modeled by Equations (1) and (2) is introduced.

RBF Definition and Collocation Method
Any function Φ that satisfies Φ(x) = φ(|x|), with φ ∈ C[0, ∞), is a radial function. This function is positive definite or m-order conditionally positive definite on R n , when in which all nonzero a ∈ R n satisfying ∑ N i=1 a i p(x i ) = 0, for all p ∈ Π m , and Π m is the set of polynomials of degree m − 1 or less. The primal space related to the nodal points X N is constructed as follows, Commonly used types of RBFs include the following forms, in which r = x − x i and the shape parameter ε controls their flatness [35].
• φ(r) = e −(εr) 2 , Gaussian RBF. Now, we briefly introduce the RBFs collocation method. Let Ω ⊆ R d and consider a boundary value problem as follows, where L is a linear differential operator and d is the dimension of the problem. We distinguish in our notation center X = {x 1 , ..., x N } and the collocation points Ξ = {α 1 , ..., α N }. Then, we have the approximate solution of Equations (6) and (7) in the form where λ i , i = 1, 2, · · · , N, are unknown coefficients that determined by collocation, φ is a RBF, . is the Euclidean norm and x i is the centers of the RBFs. Now, let Ξ divided into two subsets. One subset contains N I centers, Ξ 1 , where Equation (6) is enforced and the other subset contains N B centers, Ξ 2 , where boundary conditions are enforced. The collocation matrix is obtained by applying the collocation points to differential equation, and its boundary condition is as follows, By solving the linear system Aλ = F, we can obtain the unknown coefficients λ i , in which F is a vector included f (α i ), α i ∈ Ξ 1 , and g(α i ), α i ∈ Ξ 2 .

Application of RBF Collocation Method
For solving a DFOCP by the RBF collocation method, without loss of generality, it has to be assumed that η ≤ δ. Then, we can rewrite the problem expressed in Equations (1) and (2) as follows, subject to For simplicity and clarity, the method is only derived for Cubic RBFs and equally spaced nodes. Therefore, we choose the same number of RBF functions and collocation points (N) for the following approximation, Also, for the delay terms we have: Now, fractional derivation from the sides of Equation (11) with respect to t yields Obtaining a closed form analytic expression for the fractional derivative of a radial function may lead to a challenge. Accordingly, Mohammadi and Schaback [36] provided the required formulas for the fractional derivatives of RBFs, which allow us to use high order approximation methods for solving fractional problems. Now, we can approximate the continuous cost function described in Equation (9) with a trapezoidal quadrature rule as follows, where w i and t i are weight and nodes of integral quadrature rule, respectively. Now, by substituting Equations (11)-(15) into the problem modeled in Equations (9) and (10) and evaluating the dynamic constraints expressed in Equation (10) at the collocation nodes, we have the following NLP problem, subject to The purpose is to find Λ = (λ 1 , λ 2 , · · · , λ N ) and Γ = (γ 1 , γ 2 , · · · , γ N ) from Equation (18) such that minimize the cost function expressed in Equation (17). The described solution is called the RBF collocation method, developed as a set of MATLAB functions to transcribe the FOCP modeled in Equations (1) and (2) into an NLP optimization problem, and then use SNOPT [37] (i.e., a sparse NLP solver) to find the optimal trajectory. SNOPT uses a gradient-based optimization algorithm to solve the NLP, meaning that derivatives of cost and constraints must be provided. The proposed method has been developed in such a way that it automatically computes those gradients using the Symbolic Math Toolbox in MATLAB.

Numerical Implementation
Here, we apply the Cubic RBFs which is discussed in Section 3 for solving several DFOCPs. We test the performance of the proposed scheme on some test problems, and also present the results for different values of fractional order α and number of Cubic RBFs N. All numerical computations have been coded in Matlab R2015b on a 2.30 MHz Alpha Machine with 2GB RAM. Note that, in a minimization problem, the minimum value of the objective function is the best comparison to decide which the most efficient method is. This comparison between the proposed method and other previous algorithms can be found in the conclusion section. Moreover, comparison of these methods in terms of computational time (i.e., CPU time in seconds) is also provided along this section. Example 1. Let us consider the first DFOCP as follows, subjected to the dynamical system where 0 ≤ t ≤ 2 and x(t) = 0 at t < −1.
This problem was introduced by Moradi and Mohammadi [38], who proposed a solution based on discrete Chebyshev polynomials. More precisely, the authors solved this problem for different choices of α [26,28]. Moreover, for α = 1, Tohidi et al. [39] solved the problem using Müntz-Legendre spectral collocation method, and Ghomanjani et al. [40] used the Bezier curves for approximating the trajectory and control functions. However, the proposed RBF collocation method was more efficient than these and other previous algorithms, as Table 1 shows. From the perspective of cost values for various basis functions, our suggested approach is more effective by increasing N. Figures 1 and 2 show the graphs of x(t) and u(t), respectively, for N = 20. Moreover, these figures show that as α approaches 1, the solution for the integer order system is recovered.
In direct methods, initial guesses must be offered only for some quantities, like the states and possibly controls which are physically intuitive. As can be seen in Figure 1, the initial condition x(0) = 1 is achieved with the proposed method. By contrast, that condition was not reached in previous works [26,28,38,39], thus increasing their error.

Example 2.
Here, we consider the following FOCP with delay in control, subjected to the dynamical system where t ∈ [0, 1 4 ]. The values of J obtained by the proposed algorithm and other previous works [27,39,40] are presented in Table 2. As can be seen, the best performance was obtained by our approach, which also achieved good approximation results with small values of N. Figures 3 and 4 displays the graphs of x(t) and u(t), respectively, for N = 20. These figures corroborate the validity and efficacy of our method for this problem. Again, it can be seen that the initial condition x(0) = 1 is achieved with the proposed method, while that condition was not obtained in other previous reports.

Example 3. Consider a DFOCPs with two different delays in the form
where 0 ≤ t ≤ 1. Table 3 shows the obtained values of J for α = 1 with our scheme, Chelyshkov wavelets [27], Bernoulli polynomials [41], fractional-order Lagrange polynomials [42], Bernoulli wavelets basis [28], Müntz-Legendre polynomials [39], the least square method [40], and fractional-order Boubaker functions [29]. Again, the proposed algorithm also reported a very efficient performance. In addition, Table 4 illustrates the effect of the parameters α and N on the performance of the proposed method for this problem. In this case, we can see that good approximation results are also achieved by the proposed method with small values of N. The graphs of x(t) and u(t) with different values of α are shown in Figures 5 and 6. It should be noted that, as α approaches 1, the numerical results converge to that of an integer-order differential equation. Moreover, the initial conditions x(0) = 1 and u(0) = 0 are achieved with the proposed method, while they were not reached in [28,29].
subject to: where 0 ≤ t ≤ 2. This example have been solved by Rahimkhani et al. [28], Haddadi et al. [41], Ordukhani et al. [42], Moradi et al. [27,38], and Rabiei et al. [29], but any of them reached the initial condition x(0) = 1. A comparison of the values of J obtained by these methods and that reported by the proposed scheme is presented in Table 5.
Moreover, the effect of the parameters α and N on the proposed algorithm performance is displayed in Table 6. Both comparisons reveal that the accuracy of our method was higher than all previously proposed ones. Figures 7 and 8 show the approximation graphs of x(t) and u(t) for N = 20, respectively. We can see that, as α approaches 1, the numerical results converge to those obtained for an integer-order differential equation.  Table 6. Values of J obtained by the proposed algorithm for the problem modeled in Equation (22) when different parameters α and N are analyzed. Example 5. Consider the following DFOCP, subject to: The exact solution of this problem is unavailable. Table 7 displays the numerical results achieved by the proposed method for various values of N and α = 1, as well as for other previous algorithms dealing with the same problem. As can be seen, the obtained results corroborate the validity and efficacy of our method for this problem.

Conclusions
This paper has introduced a new technique based on the collocation method to solve DFOCPs. The proposed design first uses collocation approximations of RBFs for control and state variables in the problem. In the next step, both the context of these basis functions and a joint application of the direct method allow us to turn a DFOCP into an NLP for finally choosing the coefficients and optimal control parameters. The numerical results obtained from several non-trivial examples, with a small number of N and some values of α, confirm the efficiency, accuracy, and high performance of the proposed approximation, which would remove ill-conditioning in most systems of discrete equations. Moreover, our results have also shown that using RBFs via a collocation method bears some advantages, such as simple evaluation of fractional derivatives and delay terms of given differential equations, and less expensive of computational costs. Moreover, as the necessary conditions need not be derived, the proposed direct method does not contain the difficulties of indirect approaches for DFOCPs. Consequently, other significant merits of the proposed approach are swift calculations, ease of implementation, and robustness. Indeed, it has provided satisfactory results when a small number of RBFs has been used. To this respect, comparison of cost values for different number of nodes discloses that the accuracy of the proposed RBF collocation method is higher than most previous methods, additionally requiring less CPU time (Please see Table 8). Table 8. Summary of the values of J obtained by several algorithms for the tested problems with α = 1.