Application of Bernoulli Polynomials for Solving Variable-Order Fractional Optimal Control-Affine Problems

We propose two efficient numerical approaches for solving variable-order fractional optimal control-affine problems. The variable-order fractional derivative is considered in the Caputo sense, which together with the Riemann-Liouville integral operator is used in our new techniques. An accurate operational matrix of variable-order fractional integration for Bernoulli polynomials is introduced. Our methods proceed as follows. First, a specific approximation of the differentiation order of the state function is considered, in terms of Bernoulli polynomials. Such approximation, together with the initial conditions, help us to obtain some approximations for the other existing functions in the dynamical control-affine system. Using these approximations, and the Gauss-Legendre integration formula, the problem is reduced to a system of nonlinear algebraic equations. Some error bounds are then given for the approximate optimal state and control functions, which allow us to obtain an error bound for the approximate value of the performance index. We end by solving some test problems, which demonstrate the high accuracy of our results.


Introduction
The Bernoulli polynomials, named after Jacob Bernoulli (1654-1705), occur in the study of many special functions and, in particular, in relation with fractional calculus, which is a classical area of mathematical analysis whose foundations were laid by Liouville in a paper from 1832 and that is nowadays a very active research area [1]. One can say that Bernoulli polynomials are a powerful mathematical tool in dealing with various problems of dynamical nature [2][3][4][5][6]. Recently, an approximate method, based on orthonormal Bernoulli's polynomials, has been developed for solving fractional order differential equations of Lane-Emden type [7], while in [8] Bernoulli polynomials are used to numerical solve Fredholm fractional integro-differential equations with right-sided Caputo derivatives. Here we are interested in the use of Bernoulli polynomials with respect to fractional optimal control problems.
An optimal control problem refers to the minimization of a functional on a set of control and state variables (the performance index) subject to dynamic constraints on the states and controls. When such dynamic constraints are described by fractional differential equations, then one speaks of fractional optimal control problems (FOCPs) [9]. This is a preprint of a paper whose final and definite form is published Open Access in Axioms, available at https://www.mdpi.com/journal/axioms arXiv:2010.02833v1 [math.OC] 6 Oct 2020 2 of 19 The mathematical theory of fractional optimal control has born in 1996/97 from practical problems of mechanics and began to be developed in the context of the fractional calculus of variations [10][11][12]. Soon after, fractional optimal control theory became a mature research area, supported with many applications in engineering and physics. For the state of the art, see [13][14][15] and references therein. Regarding the use of Bernoulli polynomials to numerically solve FOCPs, we refer to [2], where the operational matrices of fractional Riemann-Liouville integration for Bernoulli polynomials are derived and the properties of Bernoulli polynomials are utilized, together with Newton's iterative method, to find approximate solutions to FOCPs. The usefulness of Bernoulli polynomials for mixed integer-fractional optimal control problems is shown in [16], while the practical relevance of the methods in engineering is illustrated in [17]. Recently, such results have been generalized for two-dimensional fractional optimal control problems, where the control system is not a fractional ordinary differential equation but a fractional partial differential equation [18].
Here we are the first to develop a numerical method, based on Bernoulli polynomials, for FOCPs of variable-order.
The variable-order fractional calculus was introduced in 1993 by Samko and Ross and deals with operators of order α, where α is not necessarily a constant but a function α(t) of time [19]. With this extension, numerous applications have been found in physics, mechanics, control, and signal processing [20][21][22][23][24]. For the state of the art on variable-order fractional optimal control we refer the interested reader to the book [25] and the articles [26,27]. To the best of our knowledge, numerical methods based on Bernoulli polynomials for such kind of FOCPs are not available in the literature. For this reason, in this work we focus on the following variable-order fractional optimal control-affine problem (FOC-AP): subject to the control-affine dynamical system and the initial conditions where φ and ϕ are smooth functions of their arguments, b = 0, n is a positive integer number such that for all is the (left) fractional derivative of variable-order defined in the Caputo sense. We employ two spectral methods based on Bernoulli polynomials in order to obtain numerical solutions to problem (1)-(3). Our main idea consists of reducing the problem to a system of nonlinear algebraic equations. To do this, we introduce an accurate operational matrix of variable-order fractional integration, having Bernoulli polynomials as basis vectors.
The paper is organized as follows. In Section 2, the variable-order fractional calculus is briefly reviewed and some properties of the Bernoulli polynomials are recalled. A new operational matrix of variable-order is introduced for the Bernoulli basis functions in Section 3. Section 4 is devoted to two new numerical approaches based on Bernoulli polynomials for solving problem (1)-(3). In Section 5, some error bounds are proved. Then, in Section 6, some FOC-APs are solved using the proposed methods. Finally, concluding remarks are given in Section 7.

Preliminaries
In this section, a brief review on necessary definitions and properties of the variable-order fractional calculus is presented. Furthermore, Bernoulli polynomials and some of their properties are recalled.

The variable-order fractional calculus
The two most commonly used definitions in fractional calculus are the Riemann-Liouville integral and the Caputo derivative. Here, we deal with generalizations of these two definitions, which allow the order of the fractional operators to be of variable-order.

Bernoulli polynomials
The set of Bernoulli polynomials, {β m (t)} ∞ m=0 , consists of a family of independent functions that builds a complete basis for the space L 2 [0, 1] of all square integrable functions on the interval [0, 1]. These polynomials are defined as where b k , k = 0, 1, . . . , m, are the Bernoulli numbers [28]. These numbers are seen in the series expansion of trigonometric functions and can be given by the following identity [29]: Thus, the first few Bernoulli numbers are given by Furthermore, the first five Bernoulli polynomials are β 0 (t) = 1, For an arbitrary function x ∈ L 2 [0, 1], we can write Therefore, an approximation of the function x can be given by where and A = [a 0 , a 1 , . . . , a M ] T .
The vector A in (9) is called the coefficient vector and can be calculated by the formula (see [2]) where ·, · is the inner product, defined for two arbitrary functions f , g ∈ L 2 [0, 1] as and D = B(t), B(t) is calculated using the following property of Bernoulli polynomials [29]: It should be noted that is a finite dimensional subspace of L 2 [0, 1] and x M , given by (9), is the best approximation of function x in X.

Operational matrix of variable-order fractional integration
In this section, we introduce an accurate operational matrix of variable-order fractional integration for Bernoulli functions. To this aim, we rewrite the Bernoulli basis vector B given by (10) in terms of the Taylor basis functions as follows: where T is the Taylor basis vector given by and Q is the change-of-basis matrix, which is obtained using (8) as Since Q is nonsingular, we can write By considering (11) and applying the left Riemann-Liouville fractional integral operator of order α(t) to the vector B(t), we get that where S is a diagonal matrix, which is obtained using (4) as follows: Finally, by using (12) in (13), we have where P , which we call the operational matrix of variable-order fractional integration α(t) for Bernoulli functions. Since Q and Q −1 are lower triangular matrices and S is also a lower triangular matrix. In the particular case of M = 2, one has

Methods of solution
In this section, we propose two approaches for solving problem (1)-(3). To do this, first we introduce Then, we may use the following two approaches to find approximations for the state and control functions, which optimize the performance index.

Approach I
In our first approach, we consider an approximation of the nth order derivative of the unknown state function x using Bernoulli polynomials. Set where A is a 1 × (M + 1) vector with unknown elements and B is the Bernoulli basis vector given by (10). Then, using the initial conditions given in (3), and equations (5), (14), and (15), we get Moreover, using (6), (14), and (15), we get and By substituting (16)-(18) into the control-affine dynamical system given by (2), we obtain an approximation of the control function as follows: Taking into consideration (16) and (19) in the performance index J, we have For the sake of simplicity, we introduce In many applications, it is difficult to compute the integral of function G[A,t]. Therefore, it is recommended to use a suitable numerical integration formula. Here, we use the Gauss-Legendre quadrature formula to obtain where t i , i = 1, 2, . . . , N, are the zeros of the Legendre polynomial of degree N, P N (t), and ω i are the corresponding weights [30], which are given by Finally, the first order necessary condition for the optimality of the performance index implies which gives a system of M + 1 nonlinear algebraic equations in terms of the M + 1 unknown elements of the vector A. By solving this system, approximations of the optimal state and control functions are, respectively, given by (16) and (19).

Approach II
In our second approach, we set Then, using (7) with β (t) ≡ 0, we obtain that Furthermore, we get Taking (22)-(24) into consideration, equation (2) gives By substituting the approximations given by (23) and (25) into the performance index, we get By introducing then this approach continues in the same way of finding the unknown parameters of the vector A as in Approach I.

Error bounds
The aim of this section is to give some error bounds for the numerical solution obtained by the proposed methods of Section 4. We present the error discussion for Approach II, which can then be easily extended to Approach I.
Suppose that x * is the optimal state function of problem (1) [31]). According to our numerical method, f M (t) = A T B(t) is the best approximation of function f in terms of the Bernoulli polynomials, that is, We recall the following lemma from [31].
and C is a positive constant independent of function f and integer M.
Since the best approximation of function f in the subspace X is unique and f M and L M ( f ) are both the best approximations of f in X, we have f M = L M ( f ). Therefore, we get that Hereafter, C denotes a positive constant independent of M and n.
Theorem 1. Suppose x * to be the exact optimal state function of problem (1)- x * (t) ∈ H µ (0, 1), with µ ≥ 0, andx be its approximation given by (23). Then, : Y → Y be the variable-order Riemann-Liouville integral operator. By definition of the norm for operators, we have In order to prove the theorem, first we show that the operator 0 I α(t) t is bounded. Since g 2 = 1, using Schwarz's inequality, we get where we have used the assumption α(t) > 0, which gives t α(t) < 1 for 0 < t ≤ 1, and a particular property of the Gamma function, which is Γ(t) > 0.8. Therefore, 0 I α(t) t is bounded. Now, using (26), and taking into account (7) and (23), we obtain that The proof is complete.
Remark 1. Since we have α(t) − α j (t) > 0, j = 1, 2, . . . , s, with a similar argument it can be shown that With the help of Theorem 1, we obtain the following result for the error of the optimal control function. For simplicity, suppose that in the control-affine dynamical system given by (2) the function ϕ appears as ϕ := ϕ(t, x) (cf. Remark 2).
Theorem 2. Suppose that the assumptions of Theorem 1 are fulfilled. Let u * andũ be the exact and approximate optimal control functions, respectively. If ϕ : R 2 −→ R satisfies a Lipschitz condition with respect to the second argument, then Proof. Using equation (2), the exact optimal control function is given by and the approximate control function obtained by our Approach II is given bỹ By subtracting (30) from (29), we get Since the function ϕ satisfies a Lipschitz condition with respect to the second variable, there exists a positive constant K such that |ϕ(t, Therefore, using (26) and (27), and also taking into account b(t) = 0, we have which yields (28).

Remark 2.
For the general case ϕ := ϕ(t, x, x 1 , . . . , x s ), the same result of Theorem 2 can be easily obtained by assuming that ϕ satisfies Lipschitz conditions with respect to the variables x, x 1 , . . . , x s .
As a result of Theorems 1 and 2, we obtain an error bound for the approximate value of the optimal performance index J given by (20). First, we recall the following lemma in order to obtain the error of the Gauss-Legendre quadrature rule.
Lemma 2 (See [30]). Let g be a given sufficiently smooth function. Then, the Gauss-Legendre quadrature rule is given by where t i , i = 1, . . . , N, are the roots of the Legendre polynomial of degree N, and ω i are the corresponding weights given by (21). In (32), E N (g) is the error term, which is given as follows: Now, by considering the assumptions of Theorems 1 and 2, we prove the following result.
Theorem 3. Let J * be the exact value of the optimal performance index J in problem (1)-(3) andJ be its approximation given by (20). Suppose that the function φ : R 3 −→ R is a sufficiently smooth function with respect to all its variables and satisfies Lipschitz conditions with respect to its second and third arguments, that is, and where K 1 and K 2 are real positive constants. Then, there exist positive constants C 1 and C 2 such that Proof. Using (20) and (32), we havẽ where t=η for η ∈ (0, 1). Therefore, taking into consideration (33)-(36), we get where we have used the property of equivalence of L 1 and L 2 -norms and The proof is complete.

Remark 3.
A similar error discussion can be considered for Approach I by setting f (t) := x * (n) (t) with f (t) ∈ H µ (0, 1) and taking into account the fact that the operators I n , I α(t) and I α j (t) , for j = 1, 2, . . . , s, are bounded.

Remark 4.
In practice, since the exact control and state functions that minimize the performance index are unknown, in order to reach a given specific accuracy ε for these functions, we increase the number of basis functions (by increasing M) in our implementation, such that

Test problems
In this section, some FOC-APs are included and solved by the proposed methods, in order to illustrate the accuracy and efficiency of the new techniques. In our implementation, the method was carried out using Mathematica 12. Furthermore, we have used N = 14 in employing the Gauss-Legendre quadrature formula. Example 1. As first example, we consider the following variable-order FOC-AP: The exact optimal state and control functions are given by which minimize the performance index J with the minimum value J = 0. In [26], a numerical method based on the Legendre wavelet has been used to solve this problem with α(t) = 1. For solving this problem with α(t) = 1, according to our methods, we have n = 1. In this case, both approaches introduced in Section 4 give the same result. By setting M = 1, we suppose that The operational matrix of variable-order fractional integration is given by Therefore, we have Moreover, using the control-affine dynamical system, we get By substituting (38) and (39) into (37), using the Gauss-Legendre quadrature for computing J, and, finally, setting we obtain a system of two nonlinear algebraic equations in terms of a 0 and a 1 . By solving this system, we find which gives the exact solution x(t) = t 2 and u(t) = te −t − 1 2 e t 2 −t .
As it is seen, in the case of α(t) = 1, our approaches give the exact solution with M = 1 (only two basis functions) compared to the method introduced in [26] based on the use of Legendre wavelets withm = 6 (six basis functions).
Since the optimal state function is a polynomial of degree 2, Approach I gives the exact solution with M = 1 for every admissible α(t). On the other hand, if α(t) = 1, then C 0 D α(t) t x(t) ∈ H 1 (0, 1). Therefore, according to the theoretical discussion and the error bound given by (35), the numerical solution given by Approach II converges to the exact solution, very slowly, that can be confirmed by the results reported in Table 1 obtained with α(t) = sin(t) and different values of M. Furthermore, by considering a different α(t), and by applying the two proposed approaches with M = 5, the numerical results for the functions x and u are displayed in Figures 1 and 2. Figure 1 displays the numerical results obtained by Approach I, while Figure 2 shows the numerical results given by Approach II. For these results, we have used Moreover, the numerical results for the performance index, obtained by our two approaches, are shown in Table 2. It can be easily seen that, in this case, Approach I gives higher accuracy results than Approach II. This is caused by the smoothness of the exact optimal state function x.   Example 2. Consider now the following FOC-AP borrowed from [32]: and the initial conditions x(0) = x (0) = 0. For this problem, the state and control functions minimize the performance index with the optimal value J = 0. We have solved this problem by both approaches.
The numerical results of applying Approach I to this problem, with different values of M, are presented in Figure 3 and Table 3. Figure 3 displays the approximate state (left) and control (right) functions obtained by M = 1, 3, 5, 7, together with the exact ones, while Table 3 reports the approximate values of the performance index. Here, we show that Approach II gives the exact solution by considering M = 1. To do this, we suppose that Therefore, we have where Using the dynamical control-affine system given by (42), we get By substituting (43) and (44) into (41), the value of the integral can be easily computed. Then, by taking into account the optimality condition, a system of nonlinear algebraic equations is obtained. Finally, by solving this system, we obtain By taking into account these values in (43) and (44), the exact optimal state and control functions are obtained. Lotfi et al. have solved this problem using an operational matrix technique based on the Legendre orthonormal functions combined with the Gauss quadrature rule. In their method, the approximate value of the minimum performance index with five basis functions has been reported as 7.82 × 10 −9 while our suggested Approach II gives the exact value only with two basis functions.   As we see, in this example, Approach II yields the exact solution with a small computational cost, while the precision of the results of Approach I increases by enlarging M. Note that here the optimal state function is not an infinitely smooth function.
Example 3. As our last example, we consider the following FOC-AP [32]: For this example, the following state and control functions minimize the performance index J with minimum value J = 0: This problem has been solved using the proposed methods with different values of M. By considering M = 1, the numerical results of Approach I are shown in Figure 4. In this case, an approximation of the performance index is obtained as J = 7.21 × 10 −1 . By choosing M = 2, according to our numerical method, we have n = 2. Therefore, we set In the continuation of the method, we find an approximation of the control function u using the control-affine dynamical system. Then, the method proceeds until solving the resulting system, which yields a 0 = 4, a 1 = 12, a 2 = 12.
These values give us the exact solution of the problem. This problem has been solved in [32] with five basis functions and the minimum value was obtained as J = 5.42 × 10 −7 while our suggested Approach I gives the exact value with only three basis functions.
In the implementation of Approach II, we consider different values of M and report the results in Table 4 and Figure 5. These results confirm that the numerical solutions converge to the exact one by increasing the value of M. Nevertheless, we see that since the exact state function x is a smooth function, it takes much less computational effort to solve this problem by using Approach I.

Conclusions
Two numerical approaches have been proposed for solving variable-order fractional optimal control-affine problems. They use an accurate operational matrix of variable-order fractional integration for Bernoulli polynomials, to give approximations of the optimal state and control functions. These approximations, along with the Gauss-Legendre quadrature formula, are used to reduce the original problem to a system of algebraic equations. An approximation of the optimal performance index and an error bound were given. Some examples have been solved to illustrate the accuracy and applicability of the new techniques. From the numerical results of Examples 1 and 3, it can be seen that our Approach I leads to very high accuracy results with a small number of basis functions for optimal control problems in which the state function that minimizes the performance index is an infinitely smooth function. Moreover, from the results of Example 2, we conclude that Approach II may give much more accurate results than Approach I in the cases that the smoothness of C 0 D α(t) t x(t) is more than x (n) (t).