A general numerical scheme for the optimal control of fractional Birkhoffian systems

This paper gives a general numerical scheme for the optimal control problem of fractional Birkhoffian systems. The fractional forced Birkhoff equations within Riemann–Liouville fractional derivatives are derived from the fractional Pfaff–Birkhoff–d’Alembert principle which includes the control as an external force term. Following the strategy of variational integrators, the fractional Pfaff–Birkhoff–d’Alembert principle is directly discretized to develop the equivalent discrete fractional forced Birkhoff equations that served as the equality constraints of the optimization problem. Together with the initial and final state constraints on the configuration space, the original optimal control problem is converted into a nonlinear optimization problem subjected to a system of algebraic constraints, which can be solved by existing algorithms. An illustrative example is given to show the efficiency and simplicity of the proposed method.


Introduction
In recent years, fractional calculus (FC, derivatives and integrals of arbitrary order) has been a hot topic in different aspects of science and engineering, such as viscoelasticity, control theory, heat conduction, electricity, biomechanics, economics, polymers, liquid crystals, chaos and fractals, etc. [1][2][3][4][5][6]. Among the diversified fields of FC, fractional optimal control problems (FOCPs, optimal control problems which contain fractional derivatives in the performance index or system dynamic constraints or both) are of great importance due to the fact that it usually has a more accurate description of dynamical systems when using fractional differential equations (FDEs) compared to classical differential equations.
Since "the Birkhoffian mechanics is the most general possible mechanics that can be constructed from the Hamiltonian mechanics via the transformation theory" [41] (see Ref. [42] and references therein for a detailed introduction of the Birkhoffian mechanics), the fractional Birkhoffian mechanics should have a wider application than the fractional Lagrangian or Hamiltonian mechanics when applied to fractional dynamics. There has been limited work in the area of fractional Birkhoffian mechanics. Ref. [43] constructed a general fractional dynamical model of fractional Birkhoffian mechanics. Ref. [44] gave a numerical scheme to solve the fractional Birkhoff equations using variational integrators. Other existing literature mostly deals with the Noether symmetries and conserved quantities for fractional Birkhoffian systems [45][46][47]. Ref. [48] studied the optimal control of a Birkhoffian system based on variational discretization. For a wider range of theoretical and practical needs, this work will focus on the FOCPs in the Birkhoffian sense whose dynamic constraints are described by the fractional forced Birkhoff equations. The formulation is given in terms of Riemann-Liouville fractional derivative and can be easily generalized to Caputo and Riesz fractional derivatives or their combinations. In order to incorporate control into the fractional Birkhoffian systems, we formulate the fractional forced Birkhoff equations starting from the fractional Pfaff-Birkhoffd'Alembert principle, which is a generalization of the fractional Pfaff-Birkhoff principle governing the fractional Birkhoff equations. Subsequently, the cost functional depending on the state variables and their fractional derivatives along with the external control forces is minimized under the dynamic constraints governed by the fractional forced Birkhoff equations.
Numerical schemes such as finite difference methods [49], orthonormal polynomials basis [50][51][52][53], collocation methods [54,55], Oustaloup's approximation [56] and so on have been applied to solve FOCPs. These indirect methods (take variations first and then discretize) deal with the resulting fractional differential equations as necessary optimality conditions derived from the Pontryagin maximum principle [57]. Another efficient way to solve the FOCPs formulated by FCVs is implemented by directly discretizing the fractional variational principles, or more precisely, the fractional Pfaff-Birkhoff-d'Alembert principle in this paper. In this way, the discrete equations derived from the fractional Pfaff-Birkhoff-d'Alembert principle are naturally an approximation of the fractional forced Birkhoff equations. This kind of numerical treatment (discretize first and then take variations), also known as the variational integrators (VIs), has been widely used to derive structure-preserving algorithms for Lagrangian and Hamiltonian mechanical systems [58,59]. Since the numerical schemes proposed in the framework of VIs usually exhibit excellent numerical behaviors, for instance, high accuracy, long-time stability, conserved quantity preservation and precise energy prediction, the main ideas of VIs have been extended to optimal control and give rise to the method of discrete mechanics and optimal control (DMOC) [60][61][62], such that the discrete solution directly inherits characteristic structural properties from the continuous one like symmetries and integrals of the motion. Benefiting from the advantages of VIs, we will propose a general framework for the discrete optimal control of the fractional Birkhoffian system by variational discretization of the Pfaff-Birkhoff-d'Alembert principle. To the best of our knowledge, there is no research on the optimal control of fractional Birkhoffian systems, so it is of great significance in regard to the formulation and numerical solution of optimal control problem for this system. The main contributions and highlights of this paper are as follows: (i) A general formulation of optimal control of fractional Birkhoffian systems is given. (ii) The strategy of variational integrators is extended to deal with the optimal control of fractional Birkhoffian systems, which gives rise to the discrete optimal control of fractional Birkhoffian systems under a direct discretization of the fractional Pfaff-Birkhoff-d'Alembert principle. (iii) The final discrete optimal control problem we formulate and need to solve is a quadratic optimization problem with constraints in terms of a system of algebraic equations, which can be easily solved by existing algorithms. (iv) An illustrated example is given to show the effectiveness of the proposed method; a few numerical experiments are made to investigate the dynamical behaviors of the fractional linear damped oscillator under control.
An outline of the paper is as follows. Section 2 introduces the basic definitions and properties of fractional calculus needed in the following sections. In Sect. 3, we formulate the forced fractional Birkhoff equations based on the fractional Pfaff-Birkhoff-d'Alembert principle. The original optimal control problem is transformed into a constrained nonlinear optimization problem. The main contribution is presented in Sect. 4. The discrete fractional forced Birkhoff equations are derived and served as constraints of the discrete cost functional, which result in a nonlinear constrained optimization problem. An example is given in Sect. 5 to illustrate our results, and some conclusions are drawn in Sect. 6.

Preliminaries
To begin with, we recall some definitions and basic properties of fractional calculus.
Definition 1 [1] Let f (t) ∈ L 1 (0, T ), the left and right Riemann-Liouville fractional integrals of function f (t) are defined by and respectively, where α > 0 is called the fractional order and (z) is the Euler gamma function defined by There are various kinds of definitions for fractional derivatives. Here we mainly focus on the Riemann-Liouville and Caputo ones for later use. and respectively, where n = [α] + 1 is the smallest integer larger than α, i.e., n − 1 α < n, and [α] represents the largest integer not greater than α. respectively.

Remark 1
The fractional integral and fractional derivative are linear operators, i.e., we have where c 1 , c 2 ∈ R, I α denotes one of the Riemann-Liouville fractional integral operators 0 I α t , t I α T , and D α denotes one of the fractional derivative operators The relationship between Riemann-Liouville and Caputo fractional derivatives is as follows: Obviously . This statement is very important for the following derivation, as we will see in the formulation of optimal control problems with fixed endpoints.
In the following, we will always assume α ∈ (0, 1] for simplicity unless emphasized. For the case of α > 1, one can still obtain similar results.
Proof The proof can be done by direct computation. Please see Ref. [1] for details.
To develop the fractional forced Birkhoff equations in the next section, we introduce the following fractional integration by parts formula.

Lemma 2 (Fractional integration by parts)[1] Sup-
Proof We start from the following calculation: where we apply Lemma 1 in the third equality.
Similarly, in the case of which is equivalent to Eq. (9). This completes the proof.

Fractional forced Birkhoff equations
The aim of optimal control is to steer a dynamical system from the initial state to the desired final state such that a quantity given in advance, for example, the control effort in our case, is minimal. Moreover, the given objective functional is optimized under the dynamical constraints governed by the dynamic equations. Mathematically, we consider a dynamical system with a 2n-dimensional configuration space M equipped with the local coordinates a = (a 1 , a 2 , . . . , a 2n ). Given a time interval [0, T ], the system is driven from an initial state a(0) = a 0 to the final state (12) is minimal.
In this work, we study the optimal control of fractional Birkhoffian systems, in which the dynamical constraints are described by following fractional forced Birkhoff equations (14). During the optimization process under control, the system is also subject to the following fractional Pfaff-Birkhoff-d'Alembert principle a(t)) are Birkhoff functions and B(t, a(t)) is the Birkhoffian of the system [42]. The conditions δa i (0) = δa i (T ) = 0 are equivalent to the fixed endpoint assumptions a(0) = a 0 , a(T ) = a T for the underlying control problem. From fractional Pfaff-Birkhoff-d'Alembert principle (13), we can derive the equations of motion satisfying the endpoints conditions, namely (13) is equivalent to the following fractional forced Birkhoff equations: Proof We prove the proposition in a straightforward way. Firstly, we have Secondly, following the same calculation in the proof of fractional integration by parts formula (11), Substituting Eq. (16) into Eq. (15) and rearranging similar items, it follows from the fractional Pfaff-Birkhoff-d'Alembert principle that vanishes for any variations δa j , j = 1, 2, ..., 2n.
Since we consider the variational problem with fixed endpoints, the formula δa j (0) = 0 in the last term of the above equation is naturally satisfied. Thus, Eq. (17) can be simplified to Finally, we obtain following the arbitrariness of variations δa j , j = 1, 2, ..., 2n and the fundamental theorem in calculus of variations, which completes the proof.
Now we can formulate the optimal control problem as follows: t a(t), F(t))dt w.r. t. a(t), F(t), As we can see from the above derivation, the external force term in the fractional Pfaff-Birkhoff-d'Alembert principle allows us to add control to the system, which means that the fractional Pfaff-Birkhoff-d'Alembert principle can provide a nice framework to deal with the optimal control problem of fractional Birkhoffian systems. Problem (20) can be seen as a nonlinear constrained optimization problem. The necessary conditions for optimality can be derived from the Pontryagin maximum principle [57], which results in a system of FDEs and can be solved by numerical methods mentioned in the introduction. In the next section, we will not adopt this kind of treatment but start from a direct discretization of the fractional variational principle.

Discrete optimal control for fractional Birkhoffian systems
In this section, we will develop a general framework based on variational integrators to give a numerical solution for the optimal control problem of fractional Birkhoffian systems. The main strategy can be summarized as follows: (i) Choose a numerical scheme to approximate the Riemann-Liouville fractional derivative in fractional Pfaff-Birkhoff-d'Alembert principle (13); (ii) Take variations of the discrete fractional Pfaff-Birkhoff-d'Alembert principle with respect to discrete variations δa k i , i = 1, 2, . . . , 2n, k = 0, 1, . . . , N , and obtain the discrete fractional forced Birkhoff equations; (iii) Compare the coefficients of variations in the derivations of both continuous and discrete fractional forced Birkhoff equations, the discrete control forces are incorporated into the formulation as an additional equality constraint; (iv) Put all the resulting equations and the initial/final state constraints together, and discretize the cost functional to form the discrete optimal control problem.
Given a time step h, we define an increasing sequence of discrete times T = t k = kh|k = 0, 1, . . . , N } by splitting the time interval [0, T ] into N subintervals of equal length. The state variable a(t) : [0, T ] → M is replaced by its discrete counterpart a d : T → M, and a k = a d (kh) is viewed as an approximation of a(kh). Similarly, the discrete force F k = F d (kh) is regarded as the approximation of the control force F(kh).
We also denote by R k i = (t k , a k ) and B k = (t k , a k ), for i = 1, 2, . . . , 2n, k = 0, 1, . . . , N as discrete Birkhoff functions and the discrete Birkhoffian, respectively. Meanwhile, the left Riemann-Liouville fractional derivative 0 D α t a i (t) at t = t k is approximated by the shifted Grünwald-Letnikov fractional derivative, i.e., we have where w 0 α = 1 and w k Here the shifted Grünwald-Letnikov fractional derivative is used to approximate the left Riemann-Liouville fractional derivative; we must point out that it also works by using a more complex numerical integration formula to deal with the Riemann-Liouville fractional derivative. However, this will introduce additional complexity, which is not conducive to the understanding of how variational integrators are used in the numerical solution for optimal control problem of fractional Birkhoffian systems and are beyond the scope of this paper. This means that higher-order approximations of the fractional derivatives, which leads to higher-order discretizations of the optimal control problem, can be done in future work. Now we have gotten all the discrete settings to move on to the next step. As mentioned above, we start from the direct discretization of the fractional Pfaff-Birkhoff-d'Alembert principle. In the discrete case, the discrete fractional Pfaff-Birkhoff-d'Alembert principle is given by ) for all k = 0, 1, . . . , N − 1 are called the left and right forces, respectively.
Parallel to the continuous case, we have the following proposition.

Remark 2
Based on a similar discussion in Ref. [44], it is easy to see that Eq. (23) is an α-order approximation of Eq. (14). Here the Grünwald-Letnikov fractional derivative is shifted by one, which is consistent with the fact that the numerical approximation is equivalent to the forward difference method when α equals one.
In particular, the cost functional in Eq. (12) on the time slice [kh, (k + 1)h] is approximated by

(t), F(t))dt,
which yields the discrete objective function In the next step, we also require that the given initial and final states a(0) = a 0 , a(T ) = a T take a part in the discrete formulation. Obviously, we can use a 0 = a 0 , a N = a T as constraints in a straightforward way. The discrete forces at the first time node can be incorporated into the formulation by comparing the coefficients of variations δa j (0) and δa 0 j . To be specific, we can set through Eq. (27) by comparison with Eq. (17).
Remark 3 When the fractional order α equals one, optimal control problem (29) we consider will reduce to the integer one in Ref. [48], where the similar results in the sense of ordinary differential equations are a particular case of this paper. Now we have completed the general framework of optimal control for fractional Birkhoffian systems based on variational integrators. The direct discretization of the fractional Pfaff-Birkhoff-d'Alembert principle gives rise to the discrete fractional forced Birkhoff equations. After the approximation of the cost functional, the constraints on the configuration, boundary conditions (28) including the discrete forces on the first time node, together with the discrete fractional forced Birkhoff equations are regarded as the equality constraints of the discrete objective functional. Thus, original optimal control problem (20) is converted into a discrete optimal control problem, or a nonlinear optimization problem with constraints in the forms of algebraic equations (29), which can be solved by existing algorithms.

Example
Consider the following equation describing the dynamics of a linear damped oscillator where γ > 0 is a given constant. We study its fractional version starting from the following fractional Pfaff action integral with the boundary conditions The resulting fractional Birkhoff equations from Eq. (31) are as follows: which is equivalent to system (30) when α = 1 if we set a 1 = x, a 2 =ȧ 1 .
Suppose that the initial state of system (32) is given by a 1 (0) = a 2 (0) = 1, the final state at T = 10 is expected to be a 1 (T ) = a 2 (T ) = 0 under the influence of the control force u(t). Following the arguments in Ref. [48], we take F = (F 1 , F 2 ) = (e γ t u(t), 0), where u is the realistic control force to be optimized such that we can realize the fractional forced Birkhoffian representation of the fractional forced dynamical system.
It is required that the control effort is minimized during the control procedure. To include the control force in the formulation, we extend fractional Pfaff action integral (31) and consider the following fractional Pfaff-Birkhoff-d'Alembert principle: Substituting the exact form of R 1 , R 2 , B into problem (20), we can obtain the following optimal control problem: To solve problem (33) numerically, the strategy listed in the beginning of Sect. 4 is performed. On every sub-interval [kh, (k + 1)h], the virtual work is approximated by So we take to be the left and right discrete forces.
Finally, we obtain the discrete optimal control problem as follows: where k = 1, 2, . . . , N − 1 and we set t −1 = −h. We solve above optimization problem (35) by the fmincon function in MATLAB, which finds the local minimum that satisfies the constraints in all the following numerical tests. In each running, we set γ = 0.05 and all initial guesses are set to 1 in the calling of fmincon. To deeply investigate how the position/velocity of the system and control effort could be influenced by different fractional orders, we also carry out a series of figures for α = 0.25, 0.5, 0.75, 0.99.
We can see from Figs. 1 and 2 that both the position and velocity of the fractional linear damped oscillator change regularly. The curves are similar to the case without imposing the control forces, which show obvious decay trends and periodic phenomena in a smooth way. In all the cases, the final states are steered to almost 0 both in variables a 1 and a 2 . This means that the discrete controls we exert are actually working as intended.
In order to examine and compare the discrete controls in different situations, we also plot the discrete  Fig. 3 for different fractional order. The control efforts or the discrete cost functionals in problem (35) are 0.0076, 0.0415, 0.1543, 0.3071 for α = 0.25, 0.50, 0.75, 0.99, respectively, which has an interesting tendency that the control efforts increase as fractional order increases. The discrete control u k +u k+1 Fig. 3 almost varies continuously. This feature is favorable to the implementation of the control in practice.
In Fig. 4, we depict the trajectories in phase space for α = 0.25, 0.5, 0.75, 0.99. The trajectory becomes smoother when α increases from 0.25 to 0.99. To analyze the periodic behavior of the dynamical systems in different fractional orders, we can see that the number of circles in a fixed time interval T = 10 decreases as α increases. This phenomenon can also be confirmed in Figs. 1 and 2 by the fact that peaks and troughs are more closely spaced in lower fractional order, which indicates that lower fractional-order systems shall oscillate faster than higher ones.
We also present some numerical experiments to further investigate the computed final position and velocity as indicated in Figs. 5 and 6, respectively. The obtained numerical results coincide closely with the given final states (a 1 (10) = a 2 (10) = 0) at an extremely small scale of at most 10 −21 , which verifies the validity and effectiveness of our method. Figure 7 shows the control effort using a different number of nodes in different fractional orders. The more closely α approaches 1, the more closely the control effort stays at the same order of magnitude. For α = 0.25, 0.50, 0.75, the discrete control effort decreases as number of nodes increases. When α = 0.99, the discrete control effort in a different number of nodes can be regarded as a constant in the sense that the magnitude of difference stays at a level of 10 −2 . This is consistent with the results in Ref. [48]. Particularly, for the case of α = 1, the discrete optimization problem formulated in (35) is exactly the one treated in Ref. [48], and we can obtain exactly the same numerical results and would not analyze it here.

Conclusions
In this paper, we develop a general numerical scheme for the optimal control of fractional Birkhoffian systems in terms of Riemann-Liouville fractional derivatives. Following the strategy of variational discretization, the discrete fractional Pfaff-Birkhoff-d'Alembert principle is used to derive the discrete fractional forced Birkhoff equations, which can be served as equality constraints of the optimization problem. The resulting discrete optimal control problem is a quadratic optimization problem with algebraic constraints, which can be solved easily by existing algorithms. This kind of treatment can avoid the discretization of the right fractional derivative that appeared in the fractional forced Birkhoff equations. It costs less computing resources and is much easier to realize by programming. In the illustrated example, we make a few numerical experiments to investigate the dynamical behaviors of a fractional linear damped oscillator under control. The control efforts have an interesting tendency to increase as the fractional order increases. The discrete control also varies continuously in all cases, which is a favorable feature in the implementation in practice. The computed final states are obtained as expected with a relatively small error on an extremely small scale, which shows the feasibility and effectiveness of the proposed scheme.