Modelling and simulation of nabla fractional order systems
using vector fitting method
Abstract
The paper focuses on the numerical approximation of discrete fractional order systems. First of all, an infinite dimensional frequency distribution model of discrete fractional order system is given. Then, resorting the nabla discrete Laplace transform, the rationality of the finite dimensional frequency distribution model approaching the infinite one is illustrated. Based on this, an original algorithm to estimate the parameters of the approximate model is proposed with the help of vector fitting method. What is more, this approximation method is extended to the approximation of the general discrete fractional order systems. Numerical examples are carried out to show the applicability and flexibility of the introduced methodology.
–Demonstration–
Modelling and simulation of nabla fractional order systems
using vector fitting method
Yiheng Wei, Jiachang Wang, Peter W Tse, Yong Wang
Key Words: Discrete fractional calculus, nabla Laplace transform, frequency distributed model, vector fitting method.
section1[Introduction]Introduction Fractional calculus is a natural extension of classical calculus, whose inception can be traced back to 300 years ago. It is well known that fractional calculus has been widely applied to system modelling [1, 2], stability analysis [3, 4] and controller synthesis [5, 6], etc. Due to the great efforts devoted by researchers, a large number of valuable results have been reported on fractional calculus. For the details of the most recent advances, one can refer to some excellent papers [7, 8, 9] and the references therein.
Despite fractional calculus is the generalization of the traditional integer order calculus, it has complex definition which brings it long memory characteristic [10]. Because of this, it is difficult to simulate fractional order systems and implement fractional order controllers [11]. Besides, fractional order systems are infinite dimensional in nature [12], which makes it challenging to get their analytical time domain response of fractional order systems especially for the nonlinear case under nonzero initial conditions. Therefore, numerical approximation comes into being and becomes essential and challenging. Fundamentally speaking, the core task of numerical approximation is to approximate the fractional calculus operators. Following this idea, many valuable results have been produced in the solution of continuous time fractional order systems [13, 14, 15, 16, 17].
It is worth mentioning that the most famous work in numerical approximation is the Oustaloup method [13]. The main idea of this method is to use a series of polylines with the slope of and to approximate the phasefrequency curve that is a straight line with the slope of . However, since the designed poles will be different with different differential orders, the dimension of the approximate model with multiple differential operator will be very large, which is not conducive to practical use. In order to reduce the dimension of the approximation model, the socalled fixedpole approximation method was proposed in [15]. In this method, the poles will be fixed for each integral order , which will bring a smaller dimensional approximate model without loss of the precision. Although the aforementioned methods have been improved, all of them are still based on curve fitting method. Since the actual curve is not a broken line and the poles are limited to real numbers, these approximation methods are flawed. What is more, all of them focused on the continuous time system and none of the discrete time system was considered.
With the rapid development of computer application, it is usually to convert a continuous time system into the discrete one by sampling. Motivated by this, some scholars have begun to pay attention to the discrete fractional calculus and the discrete fractional order system [18, 19, 20, 21]. Similarly, the numerical implementation of discrete fractional order systems is a thorny and task. To solve this problem, two discretization schemes were presented for continuous fractional differentiators and integrators [22]. With the leastsquare based methods, the discretetime integerorder IIR filter are constructed to fit the behavior in the time and frequency domains [23]. By handling discrete fractional order systems directly, the numerical scheme is given in the matrix form [24]. However, this method generally does not perform effectively when the memory length changes and higher computational complexity will be brought. By means of the singular value decompositionoriginated balanced truncation method, an approach was proposed to approximate the linear time invariant discrete time fractional order systems [25]. However, it is not convenient to select the interested frequency range and approximate degree by this approach. With further research, it has found that this issue is fundamentally an approximation problem of fractional sum operator by nabla discrete Laplace transform [26]. Inspired by this, the purpose of this paper focuses on developing a new strategy that not only approximates fractional order systems with high accuracy, but also differs from traditional curve fitting methods.
Before end up the Introduction Section, the main contributions of this paper are summarized as follows.

The inverse nabla Laplace transform is defined and proven for the first time. After wards, the method to calculate from its nabla Laplace transform is discussed in detail.

Different from the general numerical approximation methods based on curve fitting method, it is transformed into a problem of nonlinear parameter identification, and the parameters are determined by vector fitting method [27].

The proposed approximation strategy is applied to simulate discrete fractional order systems and the resulting finite dimensional state space model is not only helpful for the nonlinear case with nonzero initial conditions.
The remainder of this paper is organized as follows. Section 2 presents preliminaries, including some fundamental knowledge of nabla discrete fractional calculus and nabla discrete Laplace transform. Section 3 proposes the numerical approximation method for fractional sum operators and discrete fractional order systems. To illustrate the validity of proposed approach, three illustrative examples are contained in Section 4. Finally, some conclusions are drawn in Section 5.
section1[Preliminaries]Preliminaries In this section, the fundamental knowledge for nabla fractional calculus and the preparatory problem statement [20] will be provided.
The th nabla fractional sum of a function : can be defined by
(1) 
where , , , and is the Gamma function.
To facilitate the analysis and synthesis of nabla fractional calculus, the range of sum is changed from [28] to . With this definition, the nabla Caputo fractional difference is defined by
(2) 
where , , and represents the normal th backward difference .
The nabla Laplace transform of a function : is defined as
(3) 
Assuming , then one has
(4) 
Moreover, the initial value of can be calculated
(5) 
Along this way, a general case follows
(6) 
where . Till now, it can be concluded that given , the value of can be computed. However, as increases gradually, calculating the value of one by one has proved difficult. Taking limit operation is not an easy task either. As a smooth consequence, a more practical solution is expected.
If is a closed curve rotating around the point clockwise and it also locates in the convergent domain of , one has
(7) 
When all the poles of are inside the integration line, one can define with . Afterwards, it follows and
(8) 
With the help of the following fact
(9) 
one has
(10) 
From the existence and uniqueness of nabla Laplace transform, the inverse nabla Laplace transform can be expressed as
(11) 
It is noteworthy that if the nabla Laplace transform is given, can be calculated via the newly developed (11). Sometimes, it is difficult to solve such a contour integral problem, the corresponding equivalent realization approach can be developed.
Define , then one obtains
(12) 
for . Inspired from the continuous case, the function can be called as the discrete exponential function. With the help of such a function, residue method, partial fraction expansion method or powers series expansion method might be helpful to realizing the inverse nabla Laplace transform equivalently. In other words, no matter rational or nonrational , one can obtain the desired .
However, in many situations, is not given previous or impossible to access, for example, is the state or output of a nonlinear system. As a result, an alternative and effective method is expected.
Lemma 1
Consider the fundamental nabla discrete fractional order system
(15) 
where , are the input and output sequences, respectively. With the above basic knowledge, one can get its transfer function
(16) 
which can be as a nabla fractional sum operator in the frequency domain.
With the help of the frequency distribution model theory in [29, 30, 31], the system (15) with can be expressed equivalently as
(17) 
where , is the weighting function and is the exact system state. Interestingly, the transfer function of (17) can be expressed as
(18) 
section1[Main Results]Main Results In this section, an effective algorithm is developed to approximate the nabla fractional sum operator by the concept of identification. To deal with the nonlinear coupling problem, the vector fitting method is introduced. From this, the numerical approximation of nabla fractional order systems is further investigated.
subsection2[Approximation for nabla fractional sum]Approximation for nabla fractional sum
Due to the high complexity of the related infinite dimensional problem, the system (17) usually cannot be used directly. A practical solution is to discretize the continuous frequency range with finite distributed frequency points . Thus, an approximate finite dimensional state space model can be expressed as
(19) 
where is the weighting value and is the approximate system state. Similarly, the transfer function of the above model (19) can be obtained as follows
(20) 
Therefore, the task of numerical approximation is to use to approximate . It has been pointed by the reference [29] that if the condition that , , , are satisfied, then the system (20) would approximate the system (18) with arbitrary accuracy. In this case, the finite dimensional frequency distribution model (19) can be adopted to approximate the infinite dimensional frequency distribution model (17). Therefore, the problem of operator approximation has been transformed into a system identification problem: when the order and are known, how to determine the parameters and , so as to minimize the error between and , that is
(21) 
This is apparently a nonlinear least squares problem. To solve this problem, it is necessary to introduce two auxiliary transfer functions as
(22) 
For the purpose of obtaining the parameters and , the following equation should be guaranteed
(23) 
It can be seen from (23) that the numerator of is the denominator of , and the numerator of is the numerator of . Therefore, as long as the numerators of two auxiliary functions are determined, would be obtained. Also, (23) can be easily expanded as
(24) 
where is the known parameters when introducing the auxiliary transfer functions. Because and are the given discrete frequency points, can be regarded as the frequencydomain response. Obviously, the unknown parameters in (24) are only and , and they only exist in the numerators. Thus, the original nonlinear least squares problem in (21) is successfully transformed into the following least square problem
(25) 
where , , and the corresponding criterion function is
(26) 
Taking what is the above mentioned into account, the estimated values of and are easily acquired as follows
(27) 
where ,
It should be noted that (22) implies that (27) is a precise result only when the initial poles of the auxiliary function are chosen to be very close to the exact poles . However, one could not have any prior information when choosing auxiliary functions, so it is often difficult to obtain accurate results in one calculation. Fortunately, formula (22) also suggests that after the above calculation, the zeros of can be used as the initial poles in the next step. In accordance with this iterative relationship, the more accurate results would be obtained.
Remark 1
Note that the backwards difference is relevant to and . Then the approximate model (19) cannot be transformed into the general case , and then it cannot be solved directly by some existing modules in MATLAB. However, the obtained approximate model is causal, stable, minimumphase, and suitable for a digital implementation.
subsection2[Treatment of conjugate complex numbers]Treatment of conjugate complex numbers
One of the advantages of the above approach is that the choice of the poles and residues in the approximation model can be extended to the complex field. However, since the transfer function must be strictly regular, the chosen complex poles must be complex conjugate. Besides, consider that the calculation process of the above method is iterative, the complex poles and residuals obtained during each iteration should be conjugate. Thus, it is essential to make some improvements to the method.
In (25), assume that the th and th poles are the complex conjugate, that is
(28) 
Then, the corresponding elements , , and in vector are rewritten as
(29) 
In this way, the original least squares solution must be all the real numbers. If the estimated parameters corresponding to and are and , then the original parameters and also satisfy the following conjugate relationship
(30) 
Based on the above operation, the original least squares problem in (25) has been converted to the following form
(31) 
It is obviously that (31) only involves real numbers. The calculation results will surely guarantee conjugacy. Thus, the entire calculation process proposed above can be described as Algorithm 1.
Algorithm 1: The numerical approximation of fractional sum operator 
Prior Information: the fractional and the discrete frequency points 
Input: the total iterations number and the order of approximation model 
Output: the poles and residues 
Step 1: Let the number of iterations , and choose initial poles . 
Step 2: Use (25) to calculate and . 
Step 3: Use (27) to calculate , and combine (30) to get and . 
Step 4: Calculate the zeros of based on , and take these zeros as the poles of . 
Step 5: Based on obtained in the Step 4, use the least squares algorithm to find the residues of . 
Step 6: If , then put and as the parameters of . Otherwise, return to the Step 3 and let . 
subsection2[Approximation for nabla fractional order systems]Approximation for nabla fractional order systems After approximating the nabla fractional sum operator , the presented results can be extended to nabla discrete fractional order systems
(32) 
where the parameters satisfy
By applying the frequency distribution model theory, the finite dimensional pseudo state space system in (32) can be rewritten as an infinite dimensional exact state space model
(33) 
where , , and .
As mentioned earlier, such an infinite dimensional model can not be used directly in practice. Therefore, the system (33) needs to be approximated by a finite dimensional model as follows
(34) 
where , , , , , , the parameters and are generated by the algorithm in previous section for .
The above discussion does not consider the case of nonzero initial values. When it comes to this situation, just consider how to reasonably assign the initial pseudo state to the real state . For the discussed Caputo definition, is completely distributed on the frequency point and therefore the initial state can be configured as follows
(35) 
Remark 2
This manuscript was first completed in 2017. To facilitate academic exchange, it is attached in arXiv platform. Once it was accepted by some journal, it will be removed here.
section1[Simulation Study]Simulation Study In this section, three numerical examples will be presented to demonstrate the effectiveness of the proposed method.
Example 1: Approximation for fractional sum operator
The range of discrete frequency points is set as in this section. To illustrate the effectiveness of approximate method, the following approximate error is introduced.
(36) 
With different fractional orders , the performance of fractional sum operator approximation is studied, and the results are shown in Table 1. It can be seen that the approximation error decreases gradually as the order of approximation model increase.
To illustrate the influence of the iteration number further, the approximation for is discussed with results shown in Table 2. From a macroscopic point of view, as the number of iteration increased, the approximation accuracy has been increased. However, from a local point of view, this approach to improve the accuracy is limited. An obvious conclusion is that the better approximation performance can be obtained by determining the number of iterations around the chosen approximation order . Based on this discovery, it is good enough for the performance of approximation to set the number of iteration in the subsequent simulation as , or slightly less than . For convenience, the numbers of iteration in the following examples are set to if not stated.
As can be seen from Table 1 and Table 2, the approximation accuracy will be improved with the approximate order increased. Actually, the two tables are the comparison between and . To reflect the approximation performance of the proposed approach more intuitively, the timedomain response is considered.
Fig. 1 illustrates that the results of the proposed method and the analytical solutions are basically the same, and this conclusion can be seen more intuitively from the approximation error in Fig. 2. The analytical solution of can be obtained by utilising the definition of fractional sum introduced in Section Modelling and simulation of nabla fractional order systems using vector fitting method.
Example 2: Approximation for the linear system
Consider a system that
(38) 
where the input and . With the approximate method, the analytical output and the simulated output are displayed in Fig. 3. The approximate accuracy can be described using an error function , and it is represented in Fig. 4. The relative error is inferior to , which is quite acceptable.
Example 3: Approximation for the nonlinear system
Consider a fractional order difference equation
(39) 
where the input is the sampling point on a sawtooth wave whose frequency is Hz and amplitude is 1. The outputs solved by the definition and the approximate method are displayed in Fig. 5. The approximate error is shown in Fig. 6.
There are visible small distinction at individual sampling points, but the approximate results obtained by our method can reflect the dynamic performance of the original system to a large degree.
section1[Conclusions]Conclusions In this paper, the numerical approximation for nabla fractional order systems has been investigated. According to the frequency distribution model theory, the exact infinite dimensional state space model of fractional sum operator is given and then the approximation problem is transformed into a nonlinear identification problem for the first time. By applying the vector fitting approach, an innovative and effective numerical approximation is established. Finally, three numerical examples have verified the general applicability and flexibility of the proposed results. It is believed that the developed method will play an essential role in the applications of discrete fractional calculus.
References
 [1] Victor, S., Malti, R., Garnier, H., and Oustaloup, A., 2013. “Parameter and differentiation order estimation in fractional models”. Automatica, 49(4), pp. 926–935.
 [2] Song, B., Zheng, S. Q., Tang, X. Q., and Qiao, W. J., 2017. “Fractional order modeling and nonlinear fractional order PItype control for PMLSM system”. Asian Journal of Control, 19(2), pp. 521–531.
 [3] Chen, Y. Q., Wei, Y. H., Zhou, X., and Wang, Y., 2017. “Stability for nonlinear fractional order systems: an indirect approach”. Nonlinear Dynamics, 89(2), pp. 1011–1018.
 [4] Taghavian, H., and Tavazoei, M. S., 2018. “Algebraic conditions for stability analysis of linear timeinvariant distributed order dynamic systems: a Lagrange inversion theorem approach”. Asian Journal of Control. doi: 10.1002/asjc.1780.
 [5] Yin, C., Huang, X. G., Chen, Y. Q., Dadras, S., Zhong, S. M., and Cheng, Y. H., 2017. “Fractionalorder exponential switching technique to enhance sliding mode control”. Applied Mathematical Modelling, 44, pp. 705–726.
 [6] Nangrani, S. P., and Bhat, S. S., 2018. “Fractional order controller for controlling power system dynamic behavior”. Asian Journal of Control, 20(1), pp. 403–414.
 [7] Chen, K., Tang, R. N., Li, C., and Wei, P. N., 2018. “Robust adaptive fractionalorder observer for a class of fractionalorder nonlinear systems with unknown parameters”. Nonlinear Dynamics, 94(1), pp. 415–427.
 [8] Xue, D. Y., and Bai, L., 2017. “Benchmark problems for Caputo fractionalorder ordinary differential equations”. Fractional Calculus and Applied Analysis, 20(5), pp. 1305–1312.
 [9] Sun, H. G., Zhang, Y., Baleanu, D., Chen, W., and Chen, Y. Q., 2018. “A new collection of real world applications of fractional calculus in science and engineering”. Communications in Nonlinear Science and Numerical Simulation, 64, pp. 213–231.
 [10] Wei, Y. H., Chen, Y. Q., Cheng, S. S., and Wang, Y., 2017. “A note on short memory principle of fractional calculus”. Fractional Calculus and Applied Analysis, 20(6), pp. 1382–1404.
 [11] Zhou, X., Wei, Y. H., Liang, S., and Wang, Y., 2017. “Robust fast controller design via nonlinear fractional differential equations”. ISA Transactions, 69, pp. 20–30.
 [12] Trigeassou, J. C., Maamri, N., and Oustaloup, A., 2013. “The infinite state approach: origin and necessity”. Computers Mathematics with Applications, 66(5), pp. 892–907.
 [13] Oustaloup, A., Levron, F., Mathieu, B., and Nanot, F. M., 2000. “Frequencyband complex noninteger differentiator: characterization and synthesis”. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 47(1), pp. 25–39.
 [14] Poinot, T., and Trigeassou, J. C., 2003. “A method for modelling and simulation of fractional systems”. Signal Processing, 83(11), pp. 2319–2333.
 [15] Liang, S., Peng, C., Liao, Z., and Wang, Y., 2014. “State space approximation for general fractional order dynamic systems”. International Journal of Systems Science, 45(10), pp. 2203–2212.
 [16] Wei, Y. H., Gao, Q., Peng, C., and Wang, Y., 2014. “A rational approximate method to fractional order systems”. International Journal of Control, Automation and Systems, 12(6), pp. 1180–1186.
 [17] Du, B., Wei, Y. H., Liang, S., and Wang, Y., 2017. “Rational approximation of fractional order systems by vector fitting method”. International Journal of Control, Automation and Systems, 15(1), pp. 186–195.
 [18] Yucra, E. A., Yuz, J. I., and Goodwin, G. C., 2013. “Sampling zeros of discrete models for fractional order systems”. IEEE Transactions on Automatic Control, 58(9), pp. 2383–2388.
 [19] Ostalczyk, P., 2015. Discrete Fractional Calculus: Applications in Control and Image Processing. World Scientific, Berlin.
 [20] Goodrich, C., and Peterson, A. C., 2015. Discrete Fractional Calculus. Springer, Berlin.
 [21] Wu, G. C., Baleanu, D., and Luo, W. H., 2017. “Lyapunov functions for RiemannLiouvillelike fractional difference equations”. Applied Mathematics Computation, 314, pp. 228–236.
 [22] Chen, Y. Q., and Moore, K. L., 2002. “Discretization schemes for fractionalorder differentiators and integrators”. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 49(3), pp. 363–367.
 [23] Barbosa, R. S., and Machado, J. A. T., 2006. “Implementation of discretetime fractionalorder controllers based on LS approximations”. Acta Polytechnica Hungarica, 3(4).
 [24] Sierociuk, D., Malesza, W., and Macias, M., 2016. “Numerical schemes for initialized constant and variable fractionalorder derivatives: matrix approach and its analog verification”. Journal of Vibration and Control, 22(8), pp. 2032–2044.
 [25] Stanisławski, R., Rydel, M., and Latawiec, K. J., 2017. “Modeling of discretetime fractionalorder state space systems using the balanced truncation method”. Journal of the Franklin Institute, 354(7), pp. 3008–3020.
 [26] Atıcı, F. M., and Eloe, P. W., 2009. “Discrete fractional calculus with the nabla operator”. Electronic Journal of Qualitative Theory of Differential Equations, 2009(3), pp. 1–12.
 [27] Gustavsen, B., and Semlyen, A., 1999. “Rational approximation of frequency domain responses by vector fitting”. IEEE Transactions on Power Delivery, 14(3), pp. 1052–1061.
 [28] Cheng, J. F., 2011. Theory of Fractional Difference Equations. Xiamen University Press, Xiamen.
 [29] Montseny, G., 1998. “Diffusive representation of pseudodifferential timeoperators”. In Fractional Differential Systems: Models, Methods and Applications, pp. 159–175.
 [30] Wei, Y. H., Tse, P. W., Du, B., and Wang, Y., 2016. “An innovative fixedpole numerical approximation for fractional order systems”. ISA Transactions, 62, pp. 94–102.
 [31] Wei, Y. H., Chen, Y. Q., Wang, J. C., and Wang, Y., 2018. “Analysis and description of infinite dimensional nature for nabla discrete fractional order systems”. Communications in Nonlinear Science and Numerical Simulation. Under review.