Optimal control of system governed by nonlinear volterra integral and fractional derivative equations

This work presents a novel formulation for the numerical solution of optimal control problems related to nonlinear Volterra fractional integral equations systems. A spectral approach is implemented based on the new polynomials known as Chelyshkov polynomials. First, the properties of these polynomials are studied to solve the aforementioned problems. The operational matrices and the Galerkin method are used to discretize the continuous optimal control problems. Thereafter, some necessary conditions are defined according to which the optimal solutions of discrete problems converge to the optimal solution of the continuous ones. The applicability of the proposed approach has been illustrated through several examples. In addition, a comparison is made with other methods for showing the accuracy of the proposed one, resulting also in an improved efficiency.


Introduction
Optimal control problems (OCPs) are widely used in several fields, ranging from science and engineering to economics or biomedicine. They are essentially related to the identification of state trajectories for a dynamical system over a time interval that optimizes a specific performance index, by achieving the best possible outcome through endogenous control of a parameter within a mathematical model of the system itself. The associated problem is characterized by a cost or objective function, depending on both the state and control variables, as well as by a group of constraints. There are two important types of OCPs, respectively, subjected to differential equations and to integral equations. Originally, the classical optimal control theory was conceived to solve systems of controlled ordinary differential equations, hence referring to the first type, but the second type of OCPs recently gained significant success for handling a broad class of phenomena and mathematical models, such as, for example, technological, physical, economic, biological, and network control problems, as reported in Fig. 1.
OCPs are typically nonlinear and hence do not admit analytic solutions, especially when they are ruled by Volterra integral or Volterra integral derivative systems (second type). To overcome the difficulties related to obtaining an analytical solution to these problems, several authors have suggested different techniques providing a numerical solution. Belbas described iterative methods with their convergence by assuming some conditions on the kernel of the integral equations involved, to solve optimal control of nonlinear Volterra integral equations (VIEs) Belbas (1999). Also, he discovered a technique to solve OCPs for VIEs which are based on approximating the controlled VIEs by systems of controlled ordinary differential equationsBelbas (2007,2008). The existence and uniqueness of solutions for OCPs governed by VIEs can be found in Angell (1976).
In addition, orthogonal functions have been leveraged for finding the solution OCPs for VIEs. An iterative numerical method for solving optimal control using triangular functions is described in Maleknejad and Almasieh (2011). Maleknejad and Ebrahimzadeh introduced a collocation approach based on rationalized Legendre wavelets to approximate optimal control and state variables in Maleknejad and Ebrahimzadeh (2014). In Tohidi and Samadi (2013), Tohidi and Samadi investigated the use of Lagrange polynomials in solving OCPs for systems governed by VIEs and also analyzed the convergence of their proposed solution, characterized by a considerable efficiency, mainly for problems characterized by smooth solutions. Hybrid functions consisting of block-pulse functions and a Bernoulli polynomial method for OCPs described by integro-differential equations have been investigated by Mashayekhi et al. in Mashayekhi et al. (2013). In Peyghami et al. (2012), the authors proposed some hybrid approaches leveraging steepest descent and two-step Newton methods for achieving optimal control together with the associated optimal state. Some other methods have been described in El-Kady and Moussa (2013); Li (2010); Maleknejad et al. (2012).
In a recent paper Khanduzi et al. [22], proposed a novel revised method based on teachinglearning-based optimization (MTLBO), to gain an approximate solution of OCPs subjected to nonlinear Volterra integro-differential systems.
As said, the OCPs, which is the minimization of a performance index subject to the dynamical system, is one of the most practical subjects in science and engineering. As a generalization of the classical optimal control problems, fractional optimal control problems (FOCPs) involve the minimization of a performance index subject to dynamical systems in which fractional derivatives or integrals are used (See Moradi and Mohammadi (2019) and references therein). Even if fractional calculus is almost so old as the normal integer-order calculus, its application in various fields of science has got increasing attention in the last 3 decades. In related literature, considerable attention has been paid to fractional calculus to have a better description of the behavior of the natural processes Baleanu et al. (2016); Oldham and Spanier (1974); Samko et al. (1993); Srivastava et al. (2017).
Centered on the approach reported in [22], Maleknejad and Ebrahimzadeh (2014) and considering the interest in fractional calculus that has grown over the past few years, the main aim of this analysis is to establish a new computational method for solving the OCPs ruled by nonlinear Volterra integro-fractional differential systems (NVIFs) (1.2) Here, y(t) and u(t) are the state and control functions, D α y(t) denotes the fractional derivative of y(t) in the Caputo sense, a(t), b(t), and c(t) are functions, and ϕ is a linear or nonlinear function. Moreover, L and G are continuously differentiable operators. In this investigation, a new type of orthogonal polynomial, which has been described by Chelyshkov for the first time, is considered. First, the D α (y) is expanded by means of Chelyshkov polynomials vector with unknown coefficients. The fractional integral operational matrix is employed to find the approximate solution of OCPs (1.1) subject to the dynamic system (1.2). By increasing the number of basis functions, the accuracy of numerical results is enhanced. The novelty of this work is that in the dynamic system (1.2), we have considered the order as a fractional, while in the reported works (see Khanduzi et al. (2020),[22], Maleknejad and Ebrahimzadeh (2014) and references therein), the order of the dynamic system is considered α = 1. In fact, we have proposed the new formulation for OCPs subject to nonlinear Volterra integral equations. One of the big advantages of this approach is that by setting α = 1, our scheme can easily be applied to OCPs for NVIFs considered for examples in the work of Khanduzi et al. [22] and Maleknejad et al.Maleknejad and Ebrahimzadeh (2014), and to other similar methods. To verify this notable inference, the new technique is compared with MTLBO, TLBO, Legendre wavelet methods, and also GWO and local methods Khanduzi et al. (2020), [22], Maleknejad and Ebrahimzadeh (2014) when α = 1. Comparing the results of this work with the other relevant ones available in the related literature, as those reported by Khanduzi et al. (2020), [22], Maleknejad and Ebrahimzadeh (2014), revealed that the newly proposed formulation provides better performances with respect to the previous ones.
The rest of the paper is structured as follows: Section 2 deals with the essential definitions and notes of Riemann-Liouville and Caputo fractional derivatives. Section 3 provides the description and the properties of Chelyshkov polynomials. Section 4 evaluates the integration matrix based on these polynomials. Section 5 provides an optimization computational approach for nonlinear Volterra integro-fractional differential systems. The last section gives three numerical examples of the accuracy of the proposed numerical method. Finally, the key final remarks are outlined in Section 6.

Notations and definitions
In this section, we will briefly highlight basic definitions and some properties of fractional differentiation and integration operators . While a broad variety of problems can be modeled by fractional order operators, there is no unique definition of fractional derivatives. The definitions of Riemann Liouville and Caputo, which can be described as follows, seem to be the most widely used ones for fractional integral and derivative. For more information on the fractional derivatives and integrals, please refer to Baleanu et al. (2016); Oldham and Spanier (1974); Samko et al. (1993); Srivastava et al. (2017).

Chelyshkov polynomials
In this section, we will report the definition and some properties of Chelyshkov polynomials. These polynomials were introduced in 2006 by Chelyshkov Chelyshkov (2006). They constitute a family of new orthogonal polynomials defined by Moreover, the orthogonality condition for these polynomials is described as follows: where δ pq represents Kronecker delta.

Remark 1
By paying attention to the definition of the Chelyshkov polynomials, we conclude the main difference between the these polynomials and other orthogonal polynomials in the interval [0, 1], where the nth polynomial has a degree n.

Function approximation
Any function f(t) which is integrable on [0, 1) can be approximated by applying the Chelyshkov polynomials as where Υ (t) and C are (N + 1) vectors given by (2.6)

Operational matrices
This section concerns processing operational matrices of the Chelyshkov polynomials vector Υ (t). In the following, some explicit formulations for fractional integration operational matrix in the Riemann-Liouville sense and the product operational matrix for the Chelyshkov polynomials vectors will be given.

Theorem 1
The fractional integration of order α of Chelyshkov polynomials vector can be obtained by , and each element of this matrix can be computed as , i, j = 1, 2..., N + 1.
Proof Let us consider the ith element of the vector Υ (t). The fractional integral of order α for χ i−1 (t), can be obtained as we expand using the Chelyshkov polynomials the expression t α+r +i−1 , and then, we have where θ r , j can be obtained as Now, by substituting (3.3) and (3.4) in (3.2), we have Therefore, the desired outcome is extracted.
Theorem 2 Let Y ∈ R N ×1 be an arbitrary vector

5)
where Υ (t) ∈ R N +1 is the Chelyshkov polynomial vector introduced in (2.5) and the (i, j)th element of the product operational matrixỸ can be obtained as

Proof Consider two Chelyshkov polynomial vectors Υ (t) and Υ T (t).
The product of these two vectors is a matrix described as follows: .
As a consequence, the relation (3.5) can be represented as By multiplying χ j (t) on both sides of the relation (3.6) and integrating results over Finally, the (i, j)th element of product operational matrixỸ provided bỹ

Description of the proposed numerical method
Consider the NVIFs (1.1) with initial conditions (1.2). First of all, all functions involved in NVIFs are approximated as follows:

G(t, s, ε T Υ (s))ds).
Continuing, we can re-write u(t) in the following format using the Gauss-Legendre quadrature formula on [0, 1]: ( 4.7) where s k and w k are Gauss-Legendre quadrature weights and nodes, respectively. Therefore, the performance index (1.1) is approximated as follows: Additionally, the performance indicator J [y 0 , y 1 , ..., y N ] can be approximated by implementing the Gauss-Legendre quadrature formula on [0, 1], as follows: where w k and t k are Gauss-Legendre quadrature nodes and weights, respectively. Ultimately, the conditions required for the optimal performance indicator are We solve the algebraic equation systems for the unknown vector Y to determine the optimal coefficient values y i with i = 0, 1, ..., N . Next, we use the Newton's iterative method to evaluate the coefficients of this modified problem, which is an algebraic equation system for the unknown vector Y . By identifying the vector Y and inserting vector Y in Eqs. (4.6) and (4.7), the state and control functions y(t) and u(t) can be approximated, respectively.

Selected numerical examples and comparisons
In this section, to investigate the effectiveness of the proposed method, the numerical results based on three examples are exhibited. In these examples, the exact solutions are compared with the numerical solutions. Moreover, the obtained results are compared with the results of  Maleknejad and Ebrahimzadeh (2014). All the algorithms have been implemented using Maple 17 with 16 digits and M (number of Gauss-Legendre quadrature nodes and weights).
Example 1 Consider the following NVIFs: minimize the performance indicator subjected to the initial dynamical system For α = 1,ỹ(t) = e t 2 ,ũ(t) = 1 + 2t are the exact solutions. Hence, the solution of these NVIFs using the presented Chelyshkov polynomials-based approach for various values of N and α has been approximated. As it can be seen from Fig. 2, the approximate solutions (for N = 8, M = N + 2 and α = 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 1) have been determined. The absolute errors of the numerical solution for y(t) and u(t) for α = 1, N = 10 are also shown in Fig. 3. Table 1 summarizes the results obtained by using the presented method and the ones reported in other papers [22], Maleknejad and Ebrahimzadeh (2014) with various values of N , where α = 1 and M = N + 2. In addition, as α approaches 1, the numerical solutions converge to the exact one and agree well with it. That is, as the fractional order α approaches 1, the optimal performance indicator J gets close to the optimal value (J = 0) of the integer-order α = 1. Based on the results, it can be concluded that the approach has been very successful in solving the above problem and outperformed the other analyzed techniques.
Example 2 Consider the following NVIFs: minimize the performance indicator:  considering both state and control functions together is shown in Fig. 4, whereas the absolute errors for α = 1 and N = 10 are plotted in Fig. 5. The solution of these NVIFs using the presented Chelyshkov polynomials approach for various values of N and α has been approximated and a comparison between the obtained optimal performance indicator J results obtained with the presented method and the other ones referred in [22] Maleknejad and Ebrahimzadeh (2014) with different values of N , where α = 1 and M = N + 2 are reported in Table 2.
Based on the numerical findings presented in these tables, the utility of the method for solving NVIFs is obvious, and in contrast to other approaches, the implementation of Chelyshkov polynomials is effective and accurate. In addition, as α approaches 1, the numerical solutions converge to the exact one and agree well with it. That is, as the fractional order α approaches 1, the optimal performance indicator J get close to the optimal value (J = 0) of the integerorder α = 1. From the outcome of our investigation, it is possible to conclude that also this experiment has given good results.   Example 3 Now, consider the following NVIFs: minimize the performance indicator: subjected to the initial dynamical system Here,ỹ(t) = e t ,ũ(t) = e 3t are the exact solutions. The approximate solutions (related to N = 6, M = N + 2, and α = 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 1) are shown in Fig. 6 for both state and control functions. The absolute errors of the numerical solution for y(t) and u(t) for α = 1, N = 10 are also shown in Fig. 7. Hence, the solution of these NVIFs using the suggested Chelyshkov polynomials approaches for various values of N and α has been approximated. A comparison between the obtained optimal performance indicator J obtained with the suggested method and the other ones reported in [22], Maleknejad and Ebrahimzadeh (2014) with different values of N , where α = 1 and M = N + 2 is reported in Table 3. Based on the presented results, the utility of the method for solving NVIFs is obvious, and in contrast to other approaches, the implementation of Chelyshkov polynomials is efficient and accurate. In addition, as α approaches 1, the numerical solutions converge to the exact one and agree well with it. That is, as the fractional order α approaches 1, the optimal performance indicator J gets close to the optimal value (J = 0) of the integer-order α = 1. The findings of our research are quite convincing, and thus, it is possible to assert that the method is accurate and successful.
At the end, a comparison between the obtained optimal performance indicator J with the suggested method and the other ones reported in Khanduzi et al. (2020) for N = 7, where α = 1 and M = N + 2 is reported in Table 4 (for Examples 1, 2 and 3). As can be seen, the superiority of the method for solving NVIFs is clear that the implementation of Chelyshkov polynomials is efficient and accurate.

Conclusion
In this paper, an effective approach was introduced to approximate solutions of systems of Volterra fractional integral equations. The key characteristic of the proposed method is based on new polynomials as named Chelyshkov polynomials and their fractional operational  The absolute errors of numerical results for y(t) and u(t) for α = 1, N = 10 Example 3 matrix and it helps to reduce system of Volterra fractional integral equations into systems of algebraic equations to obtain approximate solutions. Three examples illustrating the usefulness and precision of the suggested method have been presented. In addition, a summary of our numerical findings and the numerical solutions obtained with some other methods already show that the Chelyshkov method of polynomials is more precise than other approaches. The obtained results by the proposed Chelyshkov polynomials emphasized that -The main contribution is that a new type of polynomials is applied to obtain numerical solutions. -Chelyshkov polynomials are efficient and successful for solving NVIFs. -Application of Chelyshkov polynomials is accurate and the results, as α approach to 1, are better in comparison with other reported results. -The current strategy has ended well and with good results.