A Numerical Method for a System of Fractional Di ﬀ erential-Algebraic Equations Based on Sliding Mode Control

: Fractional calculus is widely used in engineering ﬁelds. In complex mechanical systems, multi-body dynamics can be modelled by fractional di ﬀ erential-algebraic equations when considering the fractional constitutive relations of some materials. In recent years, there have been a few works about the numerical method of the fractional di ﬀ erential-algebraic equations. However, most of the methods cannot be directly applied in the equations of dynamic systems. This paper presents a numerical algorithm of fractional di ﬀ erential-algebraic equations based on the theory of sliding mode control and the fractional calculus deﬁnition of Grünwald–Letnikov. The algebraic equation is considered as the sliding mode surface. The validity of the present method is veriﬁed by comparing with an example with exact solutions. The accuracy and e ﬃ ciency of the present method are studied. It is found that the present method has very high accuracy and low time consumption. The e ﬀ ect of violation corrections on the accuracy is investigated for di ﬀ erent time steps.


Introduction
Fractional calculus (FC) means non-integer order derivative or integral of the variable function, which as an extension of the classical calculus theory, has become an important branch of differential equations [1]. It has been confirmed that many physical phenomena can be described more accurately by using fractional calculus than classical integer-order derivative [2]. Recently, the research and application of fractional calculus theory have become more and more popular in several fields of mathematics, physics, mechanics, control, electrochemistry and engineering, etc. For example, fractional calculus has been successfully applied in viscoelastic material [3][4][5], heat conduction [6], electrical spectroscopy impedance [7] and control strategies [8][9][10][11].
Fractional differential-algebraic equations (FDAEs) are used to describe many physical and engineering problems, which consist of fractional differential equations and algebraic equations. In complex mechanical systems, multi-body dynamics are usually described by differential-algebraic equations (DAEs) [12][13][14]. However, if some physical characteristics of viscoelastic damping materials are modeled by the fractional constitutive relationship, a group of FDAEs are set up for the system. In control theory, if a fractional chaotic system is controlled by the sliding mode method, the control equations of motion are FDAEs and the sliding surface corresponds to the algebraic equations. FDAEs are also applied in other fields such as electrochemistry, fractional controller, biochemistry and so engineering background and is formally an integer-order differential-algebraic equation with fractional derivatives added. Consider an affine nonlinear fractional differential-algebraic equation, given by . x(t) = Ax(t) + B 1 x (α) (t) + f(x, t) + B 2 w(t) s(x) = 0 (1) where x = (x 1 , x 2 , . . . , x n ) T is the state vector, x (α) (t) is the fractional derivatives of order α, f(x, t) is the nonlinear part, w(t) is the external input, s(x) is the algebraic equation (constraint equation), and A, B 1 and B 2 are the coefficient matrices. Equation (1) is the state equation, which is widely representative in dynamic systems. Because many dynamic systems are second-order differential equations with multiple degrees of freedom, they need to be converted into the form of state equations for analysis and numerical calculation, while the independence of state variables can still be ensured when they are converted into integer type state equations. Here, we set the range of α from 0 to 1. In fact, in the fractional-order dynamic equation or control equation, the range of fractional order is usually related to the physical characteristics of the research object. For example, the order of viscoelastic material is between 0 and 1; when the order tends towards 0, it is linear stiffness; when the order tends towards 1, it is linear damping; and the order of the control parameter may be greater than 1. However, when α is greater than 1, it is easy to reduce the order of the fractional-order equation by using the state equation, and limit α to 0 to 1. Among various definitions of fractional derivatives, the Riemann-Liouville definition is more reasonable in the dynamic description of natural science problems. However, there is a certain correlation between the discrete Grünwald-Letnikov definition and the continuous Riemann-Liouville definition. The Grünwald-Letnikov definition will degenerate to the Riemann-Liouville definition when time step tends towards 0. In this paper, we choose the Grünwald-Letnikov definition to solve the above equation, which is the most straightforward definition and convenient for numerical calculation [24]. According to the Grünwald-Letnikov definition [24], x (α) (t) can be expressed as where h is the time step, Γ() is the Gamma function and C(x, t) is given by

Transformation Based on Sliding Mode Control Theory
Sliding mode control is a variable structure control method, which can switch the structure of the nonlinear system from one structure to another along the sliding mode surface and consequently make the system asymptotically stable. To solve Equation (1), an equivalent conversion is performed for FDAEs according to the sliding mode control theory. Consider the algebraic equation s(x) = 0 as the sliding mode surface, and then convert Equation (1) into fractional differential equation, given by .
where u is the control vector added to Equation (1) and g(x) is the control matrix. g(x)u is the control term, which ensures that the state is on or close to the sliding mode surface. Typically, Equation (4) can be integrated directly by numerical methods. Meanwhile, the control matrix g(x) must satisfy the following two conditions.

1.
Based on nonlinear control theory, g(x) should satisfy the necessary and sufficient condition of the controllability for the nonlinear system Equation (4), namely rank g ad f g ad 2 f g · · · ad n−2 f g ad n−1 f g = n (5) where, ad f ·g is Lie bracket operation of vectors f = .
x and g.

2.
In order to obtain control vector u, J s g(x) should be invertible, where J s is the Jacobi matrix of s(x).
According to the theory of sliding mode control, we need to design the reaching law and the control law. The reaching law approach is adopted to make the state move towards and reach the surface s(x) = 0, which is given by [25] . s = −εs 2 sgn(s) − λs (6) where ε and λ are diagonal matrices with positive elements and sgn is the sign function. The control law can be directly solved from the reaching law. Using Equation (6), the derivative of s with respect to time can be expressed as Substituting Equation (4) into Equation (7), the control vector u can be obtained as Note that only if J s ·g(x) is invertible can u be obtained.

The Numerical Algorithm
In this section, a numerical algorithm is proposed to solve fractional differential equation Equation (4). Using discrete time domain and the Grünwald-Letnikov definition Equation (2), an approximation differential equation is obtained as a series, given by .
where A = A + B 1 /h α and t k = kh. By solving the above equation, the approximation solution of Equation (9) is obtained, In order to directly integrate the second term on the right side of the above equation, C, w and f need to be transformed appropriately. For t ∈ [t k , t k+1 ], C and w can be expressed as first-order Newton interpolation polynomials, given by The nonlinear part f(x, t) is expanded to Taylor series at t = kh and x = x(t k ), given by where Finally, using Equations (11)-(13), Equation (10) can be expressed as where the parameters H 0 , H 1 and T j are

Violation Corrections
The existence of numerical integration error results in violation, and hence, the algebraic equation will do not hold. In order to ensure the accuracy of calculation, it is necessary to modify the results of numerical calculation.
Assuming that the value of integration in step k is x k , the constraint violation will occur due to the integration error, that is s(x k ) 0. In order to improve the calculation accuracy, the value of x k needs to be corrected. Let the corrected term be ∆x k , and then the corrected value and the algebraic equation arê Expanding the above equation in Taylor series at x k up to first order, yields Because the Jacobi matrix J s is usually not a square matrix, the solution of the above equation is not unique. Since the constraint equations are independent of each other, the generalized inverse of matrix J s is The minimum norm solution is Therefore, the corrected value can be written aŝ

Results
In this section, we first use the present method to calculate an example with an exact solution to verify its validity, and then examine the accuracy and efficiency of the present method by comparing it with other methods. The effect of violation corrections is discussed at last.

An Example with Exact Solutions
Consider a nonlinear fractional order differential-algebraic equation 3 .
Its exact solutions are According to Equations (1) and (24), this can be expressed in the matrix form, then we get Using the sliding mode control theory, we convert the fractional order differential-algebraic equation into an equivalent fractional differential equation, and obtain the sliding mode surface, given by .
Considering the controllability condition, the control matrix is chosen as g(x) = (0, 1, 0) T . The parameters of sliding mode are set to be ε = 5 and λ = 5. The control vector u is consistent with the control law, Equation (8).
We start with the initial condition x(τ) = (0.001, 0.1, 0.01) T , (τ ≤ 0). Equation (27) is solved numerically by using the present method, and the numerical solutions are studied by comparing with the exact solutions in the following.
First, the numerical solutions are obtained by the present method with time step h = 0.001. Figure 1 shows the numerical solution calculated by the present method and the exact solution obtained by Equation (25). As shown in the figure, the numerical results have a good agreement with the exact results. Figure 2 shows the percentage difference between the numerical solutions and exact solutions. As shown in Figure 2, the difference is very small and decreases monotonically with the increase in t. Note that the maximum occurs at t = 1 when 1 ≤ t ≤ 20, and its magnitude is less than 7%. The range 0 ≤ t < 1 is not taken into consideration, because the exact solution tends towards 0 when t → 0, which results in very large percentage difference, even for a small error.   Next, the effect of time step on the accuracy of the present method is studied. Figures 3-5 show the percentage difference between the numerical solutions and the exact solutions for different time steps, h = 0.001, 0.005, 0.01 and 0.02. As shown in the figures, for all cases of time steps, the percentage difference increases when t → 1, but within 5% when 5 < t < 35. For the case of t > 25, the percentage difference increases with the increase in t. Note that for all solutions, the larger the step is, the larger the error is. For the case of t > 35, the accuracy of small step is still good (when h = 0.001, the maximum error at t = 50 is about 1%), while the accuracy of large step drops sharply.  Next, the effect of time step on the accuracy of the present method is studied. Figures 3-5 show the percentage difference between the numerical solutions and the exact solutions for different time steps, h = 0.001, 0.005, 0.01 and 0.02. As shown in the figures, for all cases of time steps, the percentage difference increases when t → 1, but within 5% when 5 < t < 35. For the case of t > 25, the percentage difference increases with the increase in t. Note that for all solutions, the larger the step is, the larger the error is. For the case of t > 35, the accuracy of small step is still good (when h = 0.001, the maximum error at t = 50 is about 1%), while the accuracy of large step drops sharply.  Next, the effect of time step on the accuracy of the present method is studied. Figures 3-5 show the percentage difference between the numerical solutions and the exact solutions for different time steps, h = 0.001, 0.005, 0.01 and 0.02. As shown in the figures, for all cases of time steps, the percentage difference increases when t → 1, but within 5% when 5 < t < 35. For the case of t > 25, the percentage difference increases with the increase in t. Note that for all solutions, the larger the step is, the larger the error is. For the case of t > 35, the accuracy of small step is still good (when h = 0.001, the maximum error at t = 50 is about 1%), while the accuracy of large step drops sharply.

Comparisons with other Methods
Chen et al. [26] proposed a numerical method for solving fractional differential equations by using the Grünwald-Letnikov definition. Wang et al. [20] introduced the sliding mode control method to the FDAEs, and obtained the numerical solution by using the predictor-corrector method [23]. In this section, we aim to study the characteristics of the present method by comparing with the two numerical methods.
Consider a numerical example, given by Assume that the initial condition is x(τ) = (1, 1, 1) T , (τ ≤ 0) and the order of fractional calculus is α = 0.25. According to Equation (1), the above equation can be expressed in matrix form, and hence, the coefficient matrices and the nonlinear part are

Comparisons with Other Methods
Chen et al. [26] proposed a numerical method for solving fractional differential equations by using the Grünwald-Letnikov definition. Wang et al. [20] introduced the sliding mode control method to the FDAEs, and obtained the numerical solution by using the predictor-corrector method [23]. In this section, we aim to study the characteristics of the present method by comparing with the two numerical methods.
Consider a numerical example, given by Assume that the initial condition is x(τ) = (1, 1, 1) T , (τ ≤ 0) and the order of fractional calculus is α = 0.25. According to Equation (1), the above equation can be expressed in matrix form, and hence, the coefficient matrices and the nonlinear part are Note that Equation (29) has no external input, hence w(t) = 0 and B 2 = 0. According to the sliding mode theory, we convert Equation (29) into fractional differential equation as Equation (4). The control matrix g(x) should satisfy two conditions as mentioned above, and then we let g(x) = (0, 1, 0) T . The sliding mode surface is s = x 2 − x 1×3 = 0, and we set the parameters of sliding mode ε = λ = 2.
First, for comparison, we calculate Equation (29) using the numerical method [26], the predictor-corrector method [20] and the present method, respectively. However, the numerical method [26] cannot solve Equation (29) directly, since there is an algebraic equation. Therefore, we convert Equation (29) by substituting the algebraic equation x 2 − x 1 x 3 = 0 into the fractional differential equations, and obtain The above equations are equivalent to Equation (29) and suitable for the numerical method [26]. Next, we compare the solutions calculated by the numerical method [26], the predictor-corrector method [20] and the present method. Figure 6 shows the numerical results of three methods with the time step h = 0.01. As shown in this figure, the curves of results are not in good agreement when t = 1 to 10. The results of the present method are basically smaller than those of other method. Figure 7 shows the numerical results of three methods with time step h = 0.001. From the figure, the results of three methods are in good agreement for the case of h = 0.001. Therefore, the validity of the present method can be proved. Note that the results of the two methods [20,26] are very close, even for different h. However, Figures 6 and 7 cannot examine the accuracy of the present method since there is no exact solutions for comparison. Equation (29) by substituting the algebraic equation x2 − x1x3 = 0 into the fractional differential equations, and obtain The above equations are equivalent to Equation (29) and suitable for the numerical method [26]. Next, we compare the solutions calculated by the numerical method [26], the predictor-corrector method [20] and the present method. Figure 6 shows the numerical results of three methods with the time step h = 0.01. As shown in this figure, the curves of results are not in good agreement when t = 1 to 10. The results of the present method are basically smaller than those of other method. Figure 7 shows the numerical results of three methods with time step h = 0.001. From the figure, the results of three methods are in good agreement for the case of h = 0.001. Therefore, the validity of the present method can be proved. Note that the results of the two methods [20,26] are very close, even for different h. However, Figures 6 and 7 cannot examine the accuracy of the present method since there is no exact solutions for comparison.    x x x = −    The above equations are equivalent to Equation (29) and suitable for the numerical method [26]. Next, we compare the solutions calculated by the numerical method [26], the predictor-corrector method [20] and the present method. Figure 6 shows the numerical results of three methods with the time step h = 0.01. As shown in this figure, the curves of results are not in good agreement when t = 1 to 10. The results of the present method are basically smaller than those of other method. Figure 7 shows the numerical results of three methods with time step h = 0.001. From the figure, the results of three methods are in good agreement for the case of h = 0.001. Therefore, the validity of the present method can be proved. Note that the results of the two methods [20,26] are very close, even for different h. However, Figures 6 and 7 cannot examine the accuracy of the present method since there is no exact solutions for comparison.   Next, in order to illustrate the accuracy of the present method, we compare its results with those obtained by extremely small step h. Figure 8 shows the relative error of the present method and the numerical method [26] comparing to the solution of the numerical method [26] with h = 0.0001. Because there is no exact solutions for Equation (29), the results calculated by the numerical method [26] with a sufficiently small step (h = 0.0001) are regarded as high-precision solutions. In this section, we only consider the value of x 1 and examine its accuracy. Note that for the both cases of h = 0.01 and 0.001, the present method has higher precision than the numerical method [26]. numerical method [26] comparing to the solution of the numerical method [26] with h = 0.0001. Because there is no exact solutions for Equation (29), the results calculated by the numerical method [26] with a sufficiently small step (h = 0.0001) are regarded as high-precision solutions. In this section, we only consider the value of x1 and examine its accuracy. Note that for the both cases of h = 0.01 and 0.001, the present method has higher precision than the numerical method [26].  [26] comparing to the numerical method [26] with h = 0.0001. Table 1 lists the maximum percentage differences of three methods for different steps h. It clearly indicates that the accuracy of the numerical method [26] and the predictor-corrector method [20] is equivalent. However, the present method is much more accurate than the other two methods. For the case of h = 0.001, the maximum percentage difference of the present method is less than 1%, while those of the other methods are about 4.5%. For the case of h = 0.0005, the accuracy of the present method is sufficiently high, and the error is less than 0.07%. Table 2 lists the time consumptions of the three methods for different time steps. Note that the numerical method [26] takes the longest time for the same case of h, which is about twice as long as the predictor-corrector method [20]. This illustrates that the computing efficiency is low when solving FDAEs using the numerical method of fractional differential equations. For all cases of h, the predictor-corrector method is the most efficient algorithm, and the present method is very close to it. In fact, the time differences between the two methods can be negligible in most cases. However, according to Table 1, the present method with h = 0.001 is much more accurate than the predictorcorrector method with h = 0.0005. It only takes 1.60 s for the present method to yield solutions with maximum error about 0.85%. From this point of view, the present method is the most efficient one, because it takes the shortest CPU time but obtains the highest precision. Percentage difference (%) Figure 8. Percentage difference of x 1 obtained by the present method and the numerical method [26] comparing to the numerical method [26] with h = 0.0001. Table 1 lists the maximum percentage differences of three methods for different steps h. It clearly indicates that the accuracy of the numerical method [26] and the predictor-corrector method [20] is equivalent. However, the present method is much more accurate than the other two methods. For the case of h = 0.001, the maximum percentage difference of the present method is less than 1%, while those of the other methods are about 4.5%. For the case of h = 0.0005, the accuracy of the present method is sufficiently high, and the error is less than 0.07%.  Table 2 lists the time consumptions of the three methods for different time steps. Note that the numerical method [26] takes the longest time for the same case of h, which is about twice as long as the predictor-corrector method [20]. This illustrates that the computing efficiency is low when solving FDAEs using the numerical method of fractional differential equations. For all cases of h, the predictor-corrector method is the most efficient algorithm, and the present method is very close to it. In fact, the time differences between the two methods can be negligible in most cases. However, according to Table 1, the present method with h = 0.001 is much more accurate than the predictor-corrector method with h = 0.0005. It only takes 1.60 s for the present method to yield solutions with maximum error about 0.85%. From this point of view, the present method is the most efficient one, because it takes the shortest CPU time but obtains the highest precision.

The Effect of Violation Corrections
The violation corrections can extremely improve the accuracy of the present method. In this section, we study the effect of violation corrections on the accuracy of the present method by using the example Equation (29). Additionally, we consider only the value of x 1 for comparison in this section. Figure 9 shows the comparison between the present method with violation correction and that without violation correction for different time steps. As shown in the figure, the percentage difference increases with increasing h. It illustrates that the violation corrections are more effective for large step h. Figure 10 shows the percentage difference of four methods with h = 0.01 comparing to the numerical method [26] with h = 0.0001. As shown in the figure, the present method has the highest precision under such conditions. Moreover, it can be seen that even without violation correction, the accuracy of the present method is higher than that of the other two methods. The predictor-corrector method [20] 0.05 1.10 3.66 96. 16 3 The present method 0.12 1.60 4.36 102.25

The Effect of Violation Corrections
The violation corrections can extremely improve the accuracy of the present method. In this section, we study the effect of violation corrections on the accuracy of the present method by using the example Equation (29). Additionally, we consider only the value of x1 for comparison in this section. Figure 9 shows the comparison between the present method with violation correction and that without violation correction for different time steps. As shown in the figure, the percentage difference increases with increasing h. It illustrates that the violation corrections are more effective for large step h. Figure 10 shows the percentage difference of four methods with h = 0.01 comparing to the numerical method [26] with h = 0.0001. As shown in the figure, the present method has the highest precision under such conditions. Moreover, it can be seen that even without violation correction, the accuracy of the present method is higher than that of the other two methods.   The predictor-corrector method [20] 0.05 1.10 3.66 96. 16 3 The present method 0.12 1.60 4.36 102.25

The Effect of Violation Corrections
The violation corrections can extremely improve the accuracy of the present method. In this section, we study the effect of violation corrections on the accuracy of the present method by using the example Equation (29). Additionally, we consider only the value of x1 for comparison in this section. Figure 9 shows the comparison between the present method with violation correction and that without violation correction for different time steps. As shown in the figure, the percentage difference increases with increasing h. It illustrates that the violation corrections are more effective for large step h. Figure 10 shows the percentage difference of four methods with h = 0.01 comparing to the numerical method [26] with h = 0.0001. As shown in the figure, the present method has the highest precision under such conditions. Moreover, it can be seen that even without violation correction, the accuracy of the present method is higher than that of the other two methods.  Percentage difference (%) Figure 10. Percentage difference of four methods with h = 0.01 comparing to the numerical method [26] with h = 0.0001.

Conclusions
This paper presents a numerical algorithm of FDAEs based on the theory of sliding mode control. The direct violation correction method is used to improve the accuracy of the present method. An example with exact solutions is given to illustrate the validity and precision of the present method. It proves that the solution technology of this paper can be applied in a wide range and possesses a very high accuracy. Additionally, the present method is compared to other two numerical methods. It is found that for the same time step, the present method is more accurate, and its time consumption of computing is close to the fastest method. The effect of violation corrections is investigated for the representative example of FDAEs and is found to be very effective for improving the accuracy of the present method.
Author Contributions: Software, L.W.; validation, Z.F. and J.X.; writing-review and editing, Y.T.; supervision, N.C. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the Natural Science Foundation of Jiangsu Province (China), grant number BK20160933 and the Natural Science Foundation of the Jiangsu Higher Education Institutions of China, grant number 17KJB120008.

Conflicts of Interest:
The authors declare no conflict of interest.