Modelling and simulation of nabla fractional order systemsusing vector fitting method

Modelling and simulation of nabla fractional order systems
using vector fitting method

Yiheng Wei, Jiachang Wang, Peter W Tse, Yong Wang

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.


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.

footnotetext: Manuscript received September 25, 2018.footnotetext: Y. Wei, J. Wang and Y. Wang (Γ) are with the Department of Automation, University of Science and Technology of China, Hefei, 230026, China. Email: P. W. Tse is with the Department of Systems Engineering and Engineering Management, City University of Hong Kong, Hong Kong, 999077, China.footnotetext: The work described in this paper was fully supported by the National Natural Science Foundation of China (61601431, 61573332), the Anhui Provincial Natural Science Foundation (1708085QF141), the fund of China Scholarship Council (201806345002), the Fundamental Research Funds for the Central Universities (WK2100100028) and the General Financial Grant from the China Postdoctoral Science Foundation (2016M602032).\@sect

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 phase-frequency 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 so-called fixed-pole 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 least-square based methods, the discrete-time integer-order 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 decomposition-originated 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.

  1. 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.

  2. 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].

  3. 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 non-zero 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


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


where , , and represents the normal -th backward difference .

The nabla Laplace transform of a function : is defined as


Assuming , then one has


Moreover, the initial value of can be calculated


Along this way, a general case follows


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


When all the poles of are inside the integration line, one can define with . Afterwards, it follows and


With the help of the following fact


one has


From the existence and uniqueness of nabla Laplace transform, the inverse nabla Laplace transform can be expressed as


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


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

[20] Assuming : , one can get


where , .

Consider the fundamental nabla discrete fractional order system


where , are the input and output sequences, respectively. With the above basic knowledge, one can get its transfer function


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


where , is the weighting function and is the exact system state. Interestingly, the transfer function of (17) can be expressed as


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


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


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


This is apparently a nonlinear least squares problem. To solve this problem, it is necessary to introduce two auxiliary transfer functions as


For the purpose of obtaining the parameters and , the following equation should be guaranteed


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


where is the known parameters when introducing the auxiliary transfer functions. Because and are the given discrete frequency points, can be regarded as the frequency-domain 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


where , , and the corresponding criterion function is


Taking what is the above mentioned into account, the estimated values of and are easily acquired as follows


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, minimum-phase, 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


Then, the corresponding elements , , and in vector are rewritten as


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


Based on the above operation, the original least squares problem in (25) has been converted to the following form


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


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


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


where , , , , , , the parameters and are generated by the algorithm in previous section for .

The above discussion does not consider the case of non-zero 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

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.


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.

Table 1: The error with different approximation order

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.

Table 2: The error with different number of iterations

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 time-domain response is considered.

Considering , , and the inputs


one can get the output of system (15) shown in Fig. 1.

Fig. 1: outputs of system (15) with

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.

Fig. 2: The approximation error

Example 2: Approximation for the linear system

Consider a system that


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.

Fig. 3: The outputs of system (38)
Fig. 4: the relative error

Example 3: Approximation for the nonlinear system

Consider a fractional order difference equation


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.

Fig. 5: The outputs of system (39)
Fig. 6: the relative error

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.


  • [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 PI-type 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 time-invariant 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. “Fractional-order 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 fractional-order observer for a class of fractional-order nonlinear systems with unknown parameters”. Nonlinear Dynamics, 94(1), pp. 415–427.
  • [8] Xue, D. Y., and Bai, L., 2017. “Benchmark problems for Caputo fractional-order 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. “Frequency-band 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 Riemann-Liouville-like fractional difference equations”. Applied Mathematics Computation, 314, pp. 228–236.
  • [22] Chen, Y. Q., and Moore, K. L., 2002. “Discretization schemes for fractional-order 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 discrete-time fractional-order 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 fractional-order 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 discrete-time fractional-order 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 pseudo-differential time-operators”. 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 fixed-pole 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.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description