Access the full text.

Sign up today, get DeepDyve free for 14 days.

Axioms
, Volume 9 (4) – Oct 13, 2020

/lp/multidisciplinary-digital-publishing-institute/application-of-bernoulli-polynomials-for-solving-variable-order-InfhmVzbQ7

- Publisher
- Multidisciplinary Digital Publishing Institute
- Copyright
- © 1996-2020 MDPI (Basel, Switzerland) unless otherwise stated Disclaimer The statements, opinions and data contained in the journals are solely those of the individual authors and contributors and not of the publisher and the editor(s). MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Terms and Conditions Privacy Policy
- ISSN
- 2075-1680
- DOI
- 10.3390/axioms9040114
- Publisher site
- See Article on Publisher Site

axioms Article Application of Bernoulli Polynomials for Solving Variable-Order Fractional Optimal Control-Afﬁne Problems 1 2, Somayeh Nemati and Delﬁm F. M. Torres * Department of Mathematics, Faculty of Mathematical Sciences, University of Mazandaran, Babolsar 47416-13534, Iran; s.nemati@umz.ac.ir Center for Research and Development in Mathematics and Applications (CIDMA), Department of Mathematics, University of Aveiro, 3810-193 Aveiro, Portugal * Correspondence: delﬁm@ua.pt; Tel.: +351-234-370-668 Received: 4 September 2020; Accepted: 6 October 2020; Published: 13 October 2020 Abstract: We propose two efﬁcient numerical approaches for solving variable-order fractional optimal control-afﬁne 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 speciﬁc 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-afﬁne 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. Keywords: variable-order fractional calculus; Bernoulli polynomials; optimal control-afﬁne problems; operational matrix of fractional integration MSC: 34A08; 65M70 (Primary); 11B68 (Secondary) 1. 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–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]. The mathematical theory of fractional optimal control has born Axioms 2020, 9, 114; doi:10.3390/axioms9040114 www.mdpi.com/journal/axioms Axioms 2020, 9, 114 2 of 18 in 1996/97 from practical problems of mechanics and began to be developed in the context of the fractional calculus of variations [10–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–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 ﬁnd 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 ﬁrst 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 a, where a is not necessarily a constant but a function a(t) of time [19]. With this extension, numerous applications have been found in physics, mechanics, control, and signal processing [20–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-afﬁne problem (FOC-AP): min J = f(t, x(t), u(t))dt (1) subject to the control-afﬁne dynamical system a(t) a (t) a (t) 1 s C C C D x(t) = j t, x(t), D x(t), . . . , D x(t) + b(t)u(t) (2) 0 t 0 t 0 t and the initial conditions (i) i x (0) = x , i = 0, 1, . . . , n, (3) where f and j are smooth functions of their arguments, b 6= 0, n is a positive integer number such a(t) that for all t 2 [0, 1], 0 < a (t) < a (t) < . . . < a (t) < a(t) n, and D is the (left) fractional 1 2 0 t derivative of variable-order deﬁned 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 brieﬂy 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. 2. Preliminaries In this section, a brief review on necessary deﬁnitions and properties of the variable-order fractional calculus is presented. Furthermore, Bernoulli polynomials and some of their properties are recalled. Axioms 2020, 9, 114 3 of 18 2.1. The Variable-Order Fractional Calculus The two most commonly used deﬁnitions in fractional calculus are the Riemann–Liouville integral and the Caputo derivative. Here, we deal with generalizations of these two deﬁnitions, which allow the order of the fractional operators to be of variable-order. Deﬁnition 1 (See, e.g., [25]). The left Riemann—Liouville fractional integral of order a(t) is deﬁned by a(t) a(t)1 I y(t) = (t s) y(s)ds, t > 0, G(a(t)) where G is the Euler gamma function, that is, t1 G(t) = t exp(t)dt, t > 0. Deﬁnition 2 (See, e.g., [25]). The left Caputo fractional derivative of order a(t) is deﬁned by a(t) C na(t)1 (n) D y(t) = (t s) y (s)ds, n 1 < a(t) < n, 0 t G(n a(t)) a(t) C (n) D y(t) = y (t), a(t) = n. 0 t For 0 b(t) < a(t) n, n 2 N, g > 0, and n > 1, some useful properties of the Caputo derivative and Riemann–Liouville fractional integral are as follows [25]: G(n + 1) a(t) n n+a(t) I t = t , (4) G(n + 1 + a(t)) dge1 g t C (i) I ( D y(t)) = y(t) y (0) , t > 0, (5) 0 å t 0 t i! i=0 n1 ia(t) a(t) t na(t) (n) C (i) I (y (t)) = D y(t) y (0) , t > 0, (6) t 0 t G(i + 1 a(t)) i=da(t)e da(t)e1 ib(t) a(t) b(t) t a(t)b(t) C C (i) I ( D y(t)) = D y(t) y (0) , t > 0, (7) 0 0 å t t t G(i + 1 b(t)) i=db(t)e where de is the ceiling function. 2.2. Bernoulli Polynomials The set of Bernoulli polynomials, fb (t)g , consists of a family of independent functions that m=0 builds a complete basis for the space L [0, 1] of all square integrable functions on the interval [0, 1]. These polynomials are deﬁned as b (t) = b t , (8) m mi i=0 where b , 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]: ¥ i t t = b . å i e 1 i! i=0 Thus, the ﬁrst few Bernoulli numbers are given by Axioms 2020, 9, 114 4 of 18 1 1 1 1 b = 1, b = , b = , b = 0, b = , b = 0, b = . 0 1 2 3 4 5 6 2 6 30 42 Furthermore, the ﬁrst ﬁve Bernoulli polynomials are b (t) = 1, b (t) = t , b (t) = t t + , 3 1 3 2 b (t) = t t + t, 2 2 4 3 2 b (t) = t 2t + t . For an arbitrary function x 2 L [0, 1], we can write x(t) = a b (t). å m m m=0 Therefore, an approximation of the function x can be given by x(t) ' x (t) = a b (t) = A B(t), (9) M m m m=0 where B(t) = [b (t), b (t), . . . , b (t)] (10) 0 1 M and A = [a , a , . . . , a ] . 0 1 M The vector A in (9) is called the coefﬁcient vector and can be calculated by the formula (see [2]) A = D hx(t), B(t)i, where h,i is the inner product, deﬁned for two arbitrary functions f , g 2 L [0, 1] as h f (t), g(t)i = f (t)g(t)dt, and D = hB(t), B(t)i is calculated using the following property of Bernoulli polynomials [29]: i!j! i1 b (t)b (t)dt = (1) b , i, j 1. i j i+j (i + j)! It should be noted that X = s panfb (t), b (t), . . . , b (t)g 0 1 M is a ﬁnite dimensional subspace of L [0, 1] and x , given by (9), is the best approximation of function x in X. 3. 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: B(t) = QT(t), (11) Axioms 2020, 9, 114 5 of 18 where T is the Taylor basis vector given by h i 2 M T(t) = 1, t, t , . . . , t and Q is the change-of-basis matrix, which is obtained using (8) as 2 3 1 0 0 0 0 . . . 0 6 7 1 0 0 0 . . . 0 6 7 6 7 1 1 0 0 . . . 0 6 7 6 7 Q = 1 3 . 6 0 1 0 . . . 0 7 2 2 6 7 . . . . . . 6 7 . . . . . . 4 5 . . . . . . M M M M b ( )b ( )b ( )b ( )b . . . 1 M M1 M2 M3 M4 1 2 3 4 Since Q is nonsingular, we can write T(t) = Q B(t). (12) By considering (11) and applying the left Riemann–Liouville fractional integral operator of order a(t) to the vector B(t), we get that a(t) a(t) a(t) a(t) I B(t) = I (QT(t)) = Q( I T(t)) = QS T(t), (13) 0 0 0 t t t t a(t) where S is a diagonal matrix, which is obtained using (4) as follows: 2 3 a(t) t 0 0 0 0 G(1+a(t)) 6 7 1 a(t) 6 7 0 t 0 0 0 6 G(2+a(t)) 7 6 7 a(t) a(t) 0 0 t 0 0 6 7 S = G(3+a(t)) . t 6 7 6 . . . . . 7 . . . . . 6 7 . . . . . 4 5 G( M+1) a(t) 0 0 0 0 t G( M+1+a(t)) Finally, by using (12) in (13), we have a(t) a(t) a(t) I B(t) = QS Q B(t) = P B(t), (14) t t t a(t) a(t) where P = QS Q is a matrix of dimension ( M + 1) ( M + 1), which we call the operational t t matrix of variable-order fractional integration a(t) for Bernoulli functions. Since Q and Q are a(t) a(t) lower triangular matrices and S is a diagonal matrix, P is also a lower triangular matrix. In the t t particular case of M = 2, one has 2 3 a(t) t 0 0 G(a(t)+1) 6 7 a(t) 1 1 a(t) 1 a(t) 6 7 t t 0 P = . t 2G(a(t)+2) 2G(a(t)+1) G(a(t)+2) 4 5 1 1 2 a(t) 2 1 a(t) 2 a(t) + t t t 6G(a(t)+1) 2G(a(t)+2) 3G(a(t)+3) G(a(t)+3) G(a(t)+2) G(a(t)+3) 4. Methods of Solution In this section, we propose two approaches for solving problem (1)–(3). To do this, ﬁrst we introduce n = max fda(t)eg . 0<t1 Axioms 2020, 9, 114 6 of 18 Then, we may use the following two approaches to ﬁnd approximations for the state and control functions, which optimize the performance index. 4.1. Approach I In our ﬁrst approach, we consider an approximation of the nth order derivative of the unknown state function x using Bernoulli polynomials. Set (n) T x (t) = A B(t), (15) 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 n1 i n (n) (i) x(t) = I (x (t)) + x (0) 0 å i! i=0 n1 i T n i = A ( I B(t)) + x (16) 0 å t 0 i! i=0 n1 i T n i = A P B(t) + x . t 0 i! i=0 Moreover, using (6), (14), and (15), we get n1 ia(t) a(t) t na(t) C T i D x(t) = A P B(t) + x := F[ A, t] (17) 0 t t 0 G(i + 1 a(t)) i=da(t)e and n1 ia (t) a (t) na (t) t j j C T i D x(t) = A P B(t) + x := F [ A, t], j = 1, . . . , s. (18) 0 å 0 j t t G(i + 1 a (t)) i=da (t)e By substituting (16)–(18) into the control-afﬁne dynamical system given by (2), we obtain an approximation of the control function as follows: " !# n1 i 1 t T n i u(t) = F[ A, t] j t, A P B(t) + x , F [ A, t], . . . , F [ A, t] . (19) å 1 s t 0 b(t) i! i=0 Taking into consideration (16) and (19) in the performance index J, we have " !#! n1 n1 1 i i t 1 t T n i T n i J[ A] = f t, A P B(t) + x , F[ A, t] j t, A P B(t) + x , F [ A, t], . . . , F [ A, t] dt. t å t å 1 0 0 i! b(t) i! i=0 i=0 For the sake of simplicity, we introduce " !#! n1 i n1 i t 1 t T n i T n i G[ A, t] = f t, A P B(t) + x , F[ A, t] j t, A P B(t) + x , F [ A, t], . . . , F [ A, t] . å å 1 s t 0 t 0 i! b(t) i! i=0 i=0 In many applications, it is difﬁcult 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 1 t + 1 J[ A] ' w G A, , (20) 2 2 i=1 where t , i = 1, 2, . . . , N, are the zeros of the Legendre polynomial of degree N, P (t), and w are the i N i corresponding weights [30], which are given by Axioms 2020, 9, 114 7 of 18 w = , i = 1, . . . , N. (21) d 2 P (t ) (1 t ) dt i Finally, the ﬁrst order necessary condition for the optimality of the performance index implies ¶ J[ A] = 0, ¶ A 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). 4.2. Approach II In our second approach, we set a(t) C T D x(t) = A B(t). (22) 0 t Then, using (7) with b(t) 0, we obtain that da(t)e1 a(t) t a(t) C (i) x(t) = I ( D x(t)) + x (0) 0 å t 0 G(i + 1) i=0 da(t)e1 a(t) T i (23) = A ( I B(t)) + x 0 å t 0 i! i=0 da(t)e1 a(t) T i = A P B(t) + x . t 0 i! i=0 Furthermore, we get da(t)e1 ia (t) a(t)a (t) t a (t) j j C T i D x(t) = A P B(t) + x := F [ A, t], j = 1, . . . , s. (24) å j 0 t t 0 G(i + 1 a (t)) i=da (t)e Taking (22)–(24) into consideration, Equation (2) gives " !# da(t)e1 1 t a(t) T T i u(t) = A B(t) j t, A P B(t) + x , F [ A, t], . . . , F [ A, t] . (25) å 1 s t 0 b(t) i! i=0 By substituting the approximations given by (23) and (25) into the performance index, we get da(t)e1 a(t) T i J[ A] = f t, A P B(t) + x , å 0 i! i=0 " !#! da(t)e1 1 t a(t) T T i A B(t) j t, A P B(t) + x , F [ A, t], . . . , F [ A, t] dt. å 1 t 0 b(t) i! i=0 Axioms 2020, 9, 114 8 of 18 By introducing da(t)e1 a(t) T i G[ A, t] = f t, A P B(t) + x , t 0 i! i=0 " !#! da(t)e1 1 t a(t) T T i A B(t) j t, A P B(t) + x , F [ A, t], . . . , F [ A, t] , 1 s å 0 b(t) i! i=0 then this approach continues in the same way of ﬁnding the unknown parameters of the vector A as in Approach I. 5. 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. a(t) Suppose that x is the optimal state function of problem (1)–(3). Let f (t) := D x (t) with 0 t m m T f (t) 2 H (0, 1) (H (0, 1) is a Sobolev space [31]). According to our numerical method, f (t) = A B(t) is the best approximation of function f in terms of the Bernoulli polynomials, that is, 8g 2 X, k f f k k f gk . M 2 2 We recall the following lemma from [31]. Lemma 1 (See [31]). Assume that f 2 H (0, 1) with m 0. Let L ( f ) 2 X be the truncated shifted Legendre series of f . Then, k f L ( f )k C M j fj , m;M M 2 H (0,1) where 0 1 1 (j) 2 @ A j fj = k f k m;M H (0,1) j=minfm,M+1g 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 and L ( f ) are M M both the best approximations of f in X, we have f = L ( f ). Therefore, we get that M M k f f k C M j fj m;M . (26) M 2 H (0,1) 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)–(3) such that a(t) C m f (t) := D x (t) 2 H (0, 1), with m 0, and x be its approximation given by (23). Then, 0 t kx (t) x ˜(t)k C M j fj . (27) m;M H (0,1) a(t) Proof. Let Y = L [0, 1] and I : Y ! Y be the variable-order Riemann–Liouville integral operator. By deﬁnition of the norm for operators, we have a(t) a(t) k I k = sup k I gk . 0 2 0 2 t t kgk =1 2 Axioms 2020, 9, 114 9 of 18 a(t) In order to prove the theorem, ﬁrst we show that the operator I is bounded. Since kgk = 1, 0 2 using Schwarz’s inequality, we get a(t) a(t)1 I g = (t s) g(s)ds 2 G(a(t)) a(t)1 kgk (t s) ds G(a(t)) a(t) G(a(t) + 1) C, a(t) where we have used the assumption a(t) > 0, which gives t < 1 for 0 < t 1, and a particular a(t) property of the Gamma function, which is G(t) > 0.8. Therefore, I is bounded. Now, using (26), and taking into account (7) and (23), we obtain that da(t)e1 da(t)e1 i i t t a(t) a(t) (i) T i kx (t) x ˜(t)k = I f (t) + x (0) I ( A B(t)) + x 0 å 0 å 2 t t 0 G(i + 1) i! i=0 i=0 a(t) = I ( f (t) A B(t)) a(t) I f (t) A B(t) 2 2 C M j fj . m;M H (0,1) The proof is complete. Remark 1. Since we have a(t) a (t) > 0, j = 1, 2, . . . , s, with a similar argument it can be shown that 0 1 da(t)e1 ia (t) a(t)a (t) t a (t) j j C T i m @ A D x (t) A P B(t) + x C M j fj . m;M 0 t t 0 H (0,1) G(i + 1 a (t)) i=da (t)e 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-afﬁne dynamical system given by (2) the function j appears as j := j(t, x) (cf. Remark 2). Theorem 2. Suppose that the assumptions of Theorem 1 are fulﬁlled. Let u and u ˜ be the exact and approximate optimal control functions, respectively. If j : R ! R satisﬁes a Lipschitz condition with respect to the second argument, then ku (t) u ˜(t)k C M j fj . (28) m;M H (0,1) Proof. Using Equation (2), the exact optimal control function is given by u (t) = ( f (t) j (t, x (t))) (29) b(t) and the approximate control function obtained by our Approach II is given by u ˜(t) = A B(t) j (t, x ˜(t)) . (30) b(t) Axioms 2020, 9, 114 10 of 18 By subtracting (30) from (29), we get u (t) u ˜(t) = f (t) j (t, x (t)) A B(t) + j (t, x ˜(t)) . (31) b(t) Since the function j satisﬁes a Lipschitz condition with respect to the second variable, there exists a positive constant K such that jj(t, x ) j(t, x )j < Kjx x j. 1 2 1 2 Therefore, using (26) and (27), and also taking into account b(t) 6= 0, we have T m ku (t) u ˜(t)k f (t) A B(t) + Kkx (t) x ˜(t)k C M j fj , m;M 2 H (0,1) kb(t)k 2 which yields (28). Remark 2. For the general case j := j(t, x, x , . . . , x ), the same result of Theorem 2 can be easily obtained by 1 s assuming that j satisﬁes Lipschitz conditions with respect to the variables x, x , . . . , x . 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 sufﬁciently smooth function. Then, the Gauss–Legendre quadrature rule is given by g(t)dt = w g(t ) + E (g), (32) å i i N i=1 where t , i = 1, . . . , N, are the roots of the Legendre polynomial of degree N, and w are the corresponding i i weights given by (21). In (32), E (g) is the error term, which is given as follows: 2N+1 4 2 (N!) 2N E (g) = g (h), h 2 (1, 1). (2N + 1)[(2N!)] 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) and J be its approximation given by (20). Suppose that the function f : R ! R is a sufﬁciently smooth function with respect to all its variables and satisﬁes Lipschitz conditions with respect to its second and third arguments, that is, jf(t, x , u) f(t, x , u)j K jx x j (33) 1 2 1 1 2 and jf(t, x, u ) f(t, x, u )j K ju u j, (34) 1 2 1 1 2 where K and K are real positive constants. Then, there exist positive constants C and C such that 1 2 1 2 (N!) J J C M j fj + C . (35) m;M 1 2 H (0,1) (2N + 1)[(2N!)] Proof. Using (20) and (32), we have 1 t + 1 t + 1 t + 1 i i i J = w f , x ˜ , u ˜ = f (t, x ˜(t), u ˜(t)) dt x , (36) å i 2 2 2 2 i=1 Axioms 2020, 9, 114 11 of 18 where 2N 2N+1 4 2N 4 2N 1 2 (N!) 1 ¶ f (t, x ˜(t), u ˜(t)) (N!) ¶ f (t, x ˜(t), u ˜(t)) x = = 3 2N 3 2N 2 (2N + 1)[(2N!)] 2 ¶t (2N + 1)[(2N!)] ¶t t=h t=h for h 2 (0, 1). Therefore, taking into consideration (33)–(36), we get Z Z 1 1 J J = f(t, x (t), u (t))dt f (t, x ˜(t), u ˜(t)) dt + x 0 0 Z Z Z Z 1 1 1 1 = f(t, x (t), u (t))dt f(t, x ˜(t), u (t))dt + f(t, x ˜(t), u (t))dt f (t, x ˜(t), u ˜(t)) dt + x 0 0 0 0 Z Z 1 1 K jx (t) x ˜(t)j dt + K ju (t) u ˜(t)j dt + max jx j 1 2 N 0 0 0<t<1 (N!) C M j fj m;M + C , 1 H (0,1) 2 (2N + 1)[(2N!)] 1 2 where we have used the property of equivalence of L and L -norms and 2N ¶ f (t, x ˜(t), u ˜(t)) C = max . 2N 0<t<1 ¶t The proof is complete. (n) Remark 3. A similar error discussion can be considered for Approach I by setting f (t) := x (t) with a (t) m n a(t) f (t) 2 H (0, 1) and taking into account the fact that the operators I , I and I , 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 speciﬁc accuracy e for these functions, we increase the number of basis functions (by increasing M) in our implementation, such that max jF[ A, t ] j (t , x ˜(t ), F [ A, t ], . . . , F [ A, t ]) b(t )u ˜(t )j < e (Approach I), i i i 1 i s i i i 1i M and max A B(t ) j (t , x ˜(t ), F [ A, t ], . . . , F [ A, t ]) b(t )u ˜(t ) < e (Approach II), i i i 1 i s i i i 1i M where t = , i = 1, 2, . . . , M. M + 1 6. Test Problems In this section, some FOC-APs are included and solved by the proposed methods, in order to illustrate the accuracy and efﬁciency 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 ﬁrst example, we consider the following variable-order FOC-AP: " # 1 1 2 2 2a(t) t t t min J = x(t) t + u(t) t e + e dt (37) G(3 a(t)) 2 0 Axioms 2020, 9, 114 12 of 18 subject to a(t) C x(t) t D x(t) = e + 2e u(t), 0 < a(t) 1, 0 t x(0) = 0. The exact optimal state and control functions are given by 1 1 2 2 2a(t) t t t x(t) = t , u(t) = t e e , G(3 a(t)) 2 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 a(t) = 1. For solving this problem with a(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 0 T x (t) = A B(t) = a t + a , where A = [a , a ] and B(t) = 1, t . 0 1 The operational matrix of variable-order fractional integration is given by " # t 0 P = . t t t 4 2 Therefore, we have T 1 x(t) = A P B(t) = a t + a (t 1)t. (38) t 1 Moreover, using the control-afﬁne dynamical system, we get 1 T 1 1 1 1 t T A P B(t) t a t+ a (t1)t 0 1 u(t) = e A B(t) e = e a t + a e . (39) 1 0 2 2 2 By substituting (38) and (39) into (37), using the Gauss–Legendre quadrature for computing J, and, ﬁnally, setting ¶ J ¶ J = 0, = 0, ¶a ¶a 0 1 we obtain a system of two nonlinear algebraic equations in terms of a and a . By solving this system, we ﬁnd 0 1 a = 1, a = 2, 0 1 which gives the exact solution 1 2 2 t t t x(t) = t and u(t) = te e . As it is seen, in the case of a(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 with m ˆ = 6 (six basis functions). Since the optimal state function is a polynomial of degree 2, Approach I gives the exact solution with M = 1 a(t) C 1 for every admissible a(t). On the other hand, if a(t) 6= 1, then D x(t) 2 H (0, 1). Therefore, according 0 t 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 conﬁrmed by the results reported in Table 1 obtained with a(t) = sin(t) and different values of M. Furthermore, by considering a different a(t), and by applying the Axioms 2020, 9, 114 13 of 18 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 t t a (t) = 1, a (t) = sin(t), a (t) = , a (t) = . (40) 2 3 1 4 2 3 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. Table 1. (Example 1.) Numerical results obtained by Approach II for the performance index with different M and a(t) = sin(t). M 1 2 3 2 5 3 3 3 3 3 J 6.80 10 2.33 10 1.76 10 1.57 10 1.56 10 1.0 0.16 0.15 -0.1 0.14 0.8 0.13 0.12 -0.2 0.11 0.6 0.10 0.09 -0.3 0.30 0.32 0.34 0.36 0.38 0.40 0.4 α (t) α1(t) 1 α (t) α (t) 2 2 -0.4 0.2 α (t) α (t) 3 3 α (t) α4(t) 0.0 -0.5 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 t t Figure 1. (Example 1.) Comparison between the approximate state (left) and control (right) functions obtained by Approach I with M = 5 and different a(t) (40). 1.0 0.18 -0.1 0.16 0.8 0.14 -0.2 0.12 0.6 0.10 -0.3 0.30 0.32 0.34 0.36 0.38 0.40 0.4 α (t) α (t) 1 1 α (t) α (t) 2 2 -0.4 0.2 α (t) α3(t) α (t) α4(t) 4 0.0 -0.5 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Figure 2. (Example 1.) Comparison between the approximate state (left) and control (right) functions obtained by Approach II with M = 5 and different a(t) (40). Table 2. (Example 1.) Numerical results for the performance index with M = 5 and different a(t) (40). Method a (t) a (t) a (t) a (t) 1 2 3 4 33 33 33 33 Approach I 3.05 10 3.26 10 6.89 10 2.08 10 33 3 4 5 Approach II 2.74 10 1.56 10 1.71 10 2.50 10 x(t) x(t) u(t) u(t) Axioms 2020, 9, 114 14 of 18 Example 2. Consider now the following FOC-AP borrowed from [32]: " # 5 15 p 2 6 min J = x(t) t + (1 + t ) u(t) + t t dt (41) subject to C 2 2 D x(t) = tx (t) + u(t) (42) 0 t and the initial conditions x(0) = x (0) = 0. For this problem, the state and control functions 5 15 p x(t) = t , u(t) = t + t 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 C T D x(t) = A B(t) = a t + a 1 0 with A = [a , a ] and B(t) = 1, t . 0 1 Therefore, we have 3 5 2 8 2 2 x(t) = A P B(t) = p (2a a )t + p a t , (43) 0 1 1 3 p 15 p where 2 3 p 2 3 t 0 2 3 p 4 5 P = . 3 3 2 8 p 2 p 2 t t 5 p 15 p Using the dynamical control-afﬁne system given by (42), we get 1 2 3 8 5 2 2 u(t) = a t + a t p (2a a )t + p a t 1 0 0 1 1 2 3 p 15 p ! ! (44) 2 2 2 2 64a 32a 64a a 16a a 16a 4a a 0 1 0 1 1 1 6 1 5 0 1 4 = t + t + t + a t + a . 1 0 225p 45p 45p 9p 9p 9p 2 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 p p 15 p 15 p a = , a = . 0 1 16 8 By taking into account these values in (43) and (44), the exact optimal state and control functions are obtained. Lotﬁ 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 ﬁve basis functions has been reported as 7.82 10 while our suggested Approach II gives the exact value only with two basis functions. Axioms 2020, 9, 114 15 of 18 2.5 Exact 1.0 0.12 0.11 M=1 0.10 2.0 M=3 0.8 0.09 0.08 M=5 0.07 1.5 0.6 M=7 0.06 2.20 0.05 2.15 Exact 0.30 0.32 0.34 0.36 0.38 0.40 1.0 0.4 2.10 M=1 2.05 M=3 0.5 0.2 2.00 M=5 1.95 M=7 0.0 0.0 0.60 0.62 0.64 0.66 0.68 0.70 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 t t Figure 3. (Example 2.) Comparison between the exact and approximate state (left) and control (right) functions obtained by Approach I with different values of M. Table 3. (Example 2.) Numerical results for the performance index obtained by Approach I with different M. M 1 3 5 7 4 6 7 8 J 5.24 10 7.59 10 4.65 10 5.86 10 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 inﬁnitely smooth function. Example 3. As our last example, we consider the following FOC-AP [32]: 2 3 0 1 8000 21 6 7 t 4 2 4 @ A min J = e x(t) t + t 1 + (1 + t ) u(t) + 1 t + t t dt 4 5 77G subject to 1.9 D x(t) = x(t) + u(t), 0 t x(0) = 1, x (0) = 1. For this example, the following state and control functions minimize the performance index J with minimum value J = 0: 8000 21 4 4 x(t) = t t + 1, u(t) = t + t + t 1. 77G 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 . By choosing M = 2, according to our numerical method, we have n = 2. Therefore, we set 00 T x (t) = A B(t), where 1 1 A = [a , a , a ], B(t) = 1, t , t t + . 0 1 2 2 6 Hence, using the initial conditions, the state function can be approximated by a a a a a a 2 1 2 0 1 2 T 2 4 3 2 x(t) = A P B(t) t + 1 = t + t + + t t + 1, 12 6 2 4 12 x(t) u(t) Axioms 2020, 9, 114 16 of 18 where 2 3 0 0 2 2 6 7 t t P = . 4 0 5 6 6 2 2 2 t t t 36 12 12 In the continuation of the method, we ﬁnd an approximation of the control function u using the control-afﬁne dynamical system. Then, the method proceeds until solving the resulting system, which yields a = 4, a = 12, a = 12. 0 1 2 These values give us the exact solution of the problem. This problem has been solved in [32] with ﬁve basis functions and the minimum value was obtained as J = 5.42 10 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 conﬁrm 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. Table 4. (Example 3.) Numerical results for the performance index obtained by Approach II with different M. M 2 4 6 8 4 7 8 10 J 3.79 10 5.42 10 1.21 10 7.36 10 1.0 Exact Exact M=1 M=1 0.9 0.8 0.7 0.6 0.5 -2 0.4 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 t t Figure 4. (Example 3.) Comparison between the exact and approximate state (left) and control (right) functions obtained by Approach I with M = 1. 1.0 1.06 0.626 1.04 0.624 0.9 1.02 0.622 1.00 0.98 0.620 0.8 0.96 0.618 4 Exact Exact 0.4 0.405 0.41 0.4 0.405 0.41 0.7 M=2 M=2 M=4 M=4 0.6 M=6 M=6 M=8 M=8 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 t t Figure 5. (Example 3.) Comparison between the exact and approximate state (left) and control (right) functions obtained by Approach II with different values of M. x(t) x(t) u(t) u(t) Axioms 2020, 9, 114 17 of 18 7. Conclusions Two numerical approaches have been proposed for solving variable-order fractional optimal control-afﬁne 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 inﬁnitely smooth function. Moreover, from the results of Example 2, we conclude that Approach II may give a(t) much more accurate results than Approach I in the cases that the smoothness of D x(t) is more 0 t (n) than x (t). Author Contributions: Conceptualization, S.N. and D.F.M.T.; investigation, S.N. and D.F.M.T.; software, S.N.; validation, D.F.M.T.; writing—original draft, S.N. and D.F.M.T.; writing—review and editing, S.N. and D.F.M.T. All authors have read and agreed to the published version of the manuscript. Funding: Torres was funded by Fundação para a Ciência e a Tecnologia (FCT, the Portuguese Foundation for Science and Technology) through grant UIDB/04106/2020 (CIDMA). Acknowledgments: This research was initiated during a visit of Nemati to the Department of Mathematics of the University of Aveiro (DMat-UA) and to the R&D unit CIDMA, Portugal. The hospitality of the host institution is here gratefully acknowledged. The authors are also grateful to three anonymous reviewers for several questions and comments, which helped them to improve the manuscript. Conﬂicts of Interest: The authors declare no conﬂict of interest. References 1. Machado, J.T.; Kiryakova, V.; Mainardi, F. Recent history of fractional calculus. Commun. Nonlinear Sci. Numer. Simul. 2011, 16, 1140–1153. [CrossRef] 2. Keshavarz, E.; Ordokhani, Y.; Razzaghi, M. A numerical solution for fractional optimal control problems via Bernoulli polynomials. J. Vib. Control 2016, 22, 3889–3903. [CrossRef] 3. Bhrawy, A.H.; Tohidi, E.; Soleymani, F. A new Bernoulli matrix method for solving high-order linear and nonlinear Fredholm integro-differential equations with piecewise intervals. Appl. Math. Comput. 2012, 219, 482–497. [CrossRef] 4. Tohidi, E.; Bhrawy, A.H.; Erfani, K. A collocation method based on Bernoulli operational matrix for numerical solution of generalized pantograph equation. Appl. Math. Model. 2013, 37, 4283–4294. [CrossRef] 5. Toutounian, F.; Tohidi, E. A new Bernoulli matrix method for solving second order linear partial differential equations with the convergence analysis. Appl. Math. Comput. 2013, 223, 298–310. [CrossRef] 6. Bazm, S. Bernoulli polynomials for the numerical solution of some classes of linear and nonlinear integral equations. J. Comput. Appl. Math. 2015, 275, 44–60. [CrossRef] 7. Sahu, P.K.; Mallick, B. Approximate solution of fractional order Lane-Emden type differential equation by orthonormal Bernoulli’s polynomials. Int. J. Appl. Comput. Math. 2019, 5, 89. [CrossRef] 8. Loh, J.R.; Phang, C. Numerical solution of Fredholm fractional integro-differential equation with right-sided Caputo’s derivative using Bernoulli polynomials operational matrix of fractional derivative. Mediterr. J. Math. 2019, 16, 28. [CrossRef] 9. Rosa, S.; Torres, D.F.M. Optimal control of a fractional order epidemic model with application to human respiratory syncytial virus infection. Chaos Solitons Fractals 2018, 117, 142–149. [CrossRef] 10. Malinowska, A.B.; Torres, D.F.M. Introduction to the Fractional Calculus of Variations; Imperial College Press: London, UK, 2012. [CrossRef] 11. Malinowska, A.B.; Odzijewicz, T.; Torres, D.F.M. Advanced methods in the fractional calculus of variations. In Springer Briefs in Applied Sciences and Technology; Springer: Cham, Switzerland, 2015. 12. Almeida, R.; Pooseh, S.; Torres, D.F.M. Computational Methods in the Fractional Calculus of Variations; Imperial College Press: London, UK, 2015. [CrossRef] Axioms 2020, 9, 114 18 of 18 13. Ali, M.S.; Shamsi, M.; Khosravian-Arab, H.; Torres, D.F.M.; Bozorgnia, F. A space-time pseudospectral discretization method for solving diffusion optimal control problems with two-sided fractional derivatives. J. Vib. Control 2019, 25, 1080–1095. [CrossRef] 14. Nemati, S.; Lima, P.M.; Torres, D.F.M. A numerical approach for solving fractional optimal control problems using modiﬁed hat functions. Commun. Nonlinear Sci. Numer. Simul. 2019, 78, 104849. [CrossRef] 15. Salati, A.B.; Shamsi, M.; Torres, D.F.M. Direct transcription methods based on fractional integral approximation formulas for solving nonlinear fractional optimal control problems. Commun. Nonlinear Sci. Numer. Simul. 2019, 67, 334–350. [CrossRef] 16. Rabiei, K.; Ordokhani, Y.; Babolian, E. Numerical solution of 1D and 2D fractional optimal control of system via Bernoulli polynomials. Int. J. Appl. Comput. Math. 2018, 4, 7. [CrossRef] 17. Behroozifar, M.; Habibi, N. A numerical approach for solving a class of fractional optimal control problems via operational matrix Bernoulli polynomials. J. Vib. Control 2018, 24, 2494–2511. [CrossRef] 18. Rahimkhani, P.; Ordokhani, Y. Generalized fractional-order Bernoulli-Legendre functions: An effective tool for solving two-dimensional fractional optimal control problems. IMA J. Math. Control Inf. 2019, 36, 185–212. [CrossRef] 19. Samko, S.G.; Ross, B. Integration and differentiation to a variable fractional order. Integral Transform. Spec. Funct. 1993, 1, 277–300. [CrossRef] 20. Lorenzo, C.F.; Hartley, T.T. Variable order and distributed order fractional operators. Nonlinear Dyn. 2002, 29, 57–98. [CrossRef] 21. Abdeljawad, T.; Mert, R.; Torres, D.F.M. Variable order Mittag-Lefﬂer fractional operators on isolated time scales and application to the calculus of variations. In Fractional Derivatives with Mittag-Lefﬂer Kernel; Springer: Cham, Switzerland, 2019; Volume 194, pp. 35–47. 22. Hassani, H.; Naraghirad, E. A new computational method based on optimization scheme for solving variable-order time fractional Burgers’ equation. Math. Comput. Simul. 2019, 162, 1–17. [CrossRef] 23. Odzijewicz, T.; Malinowska, A.B.; Torres, D.F.M. Fractional variational calculus of variable order. In Advances in Harmonic Analysis and Operator Theory; Birkhäuser/Springer Basel AG: Basel, Switzerland, 2013; Volume 229, pp. 291–301. [CrossRef] 24. Yan, R.; Han, M.; Ma, Q.; Ding, X. A spectral collocation method for nonlinear fractional initial value problems with a variable-order fractional derivative. Comput. Appl. Math. 2019, 38, 66. [CrossRef] 25. Almeida, R.; Tavares, D.; Torres, D.F.M. The variable-order fractional calculus of variations. In Springer Briefs in Applied Sciences and Technology; Springer: Cham, Switzerland, 2019. [CrossRef] 26. Heydari, M.H.; Avazzadeh, Z. A new wavelet method for variable-order fractional optimal control problems. Asian J. Control 2018, 20, 1804–1817. [CrossRef] 27. Mohammadi, F.; Hassani, H. Numerical solution of two-dimensional variable-order fractional optimal control problem by generalized polynomial basis. J. Optim. Theory Appl. 2019, 180, 536–555. [CrossRef] 28. Costabile, F.; Dell’Accio, F.; Gualtieri, M.I. A new approach to Bernoulli polynomials. Rend. Mat. Appl. 2006, 26, 1–12. 29. Arfken, G. Mathematical Methods for Physicists; Academic Press: New York, NY, USA; London, UK, 1966. 30. Shen, J.; Tang, T.; Wang, L.L. Spectral Methods; Springer Series in Computational Mathematics; Springer: Heidelberg, Germany, 2011; Volume 41. [CrossRef] 31. Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. Spectral Methods. In Scientiﬁc Computation; Springer: Berlin, Germany, 2006. 32. Lotﬁ, A.; Youseﬁ, S.A.; Dehghan, M. 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. J. Comput. Appl. Math. 2013, 250, 143–160. [CrossRef] c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Axioms – Multidisciplinary Digital Publishing Institute

**Published: ** Oct 13, 2020

Loading...

You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!

Read and print from thousands of top scholarly journals.

System error. Please try again!

Already have an account? Log in

Bookmark this article. You can see your Bookmarks on your DeepDyve Library.

To save an article, **log in** first, or **sign up** for a DeepDyve account if you don’t already have one.

Copy and paste the desired citation format or use the link below to download a file formatted for EndNote

Access the full text.

Sign up today, get DeepDyve free for 14 days.

All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.