Modelling and simulation of nabla fractional dynamic systems with nonzero initial conditions

TL;DR

This paper develops a numerical approach for nabla fractional dynamic systems with nonzero initial conditions, introducing inverse nabla Laplace transform, frequency distributed models, and parameter estimation algorithms, validated through numerical examples.

Contribution

It introduces a novel methodology combining inverse nabla Laplace transform, frequency distributed modeling, and vector fitting for fractional systems with nonzero initial conditions.

Findings

Effective approximation of nabla fractional systems demonstrated

Parameter estimation algorithm shows high accuracy

Method applicable to general systems beyond sum operators

Abstract

The paper focuses on the numerical approximation of nabla fractional order systems with the conditions of nonzero initial instant and nonzero initial state. First, the inverse nabla Laplace transform is developed and the equivalent infinite dimensional frequency distributed models of discrete fractional order system are introduced. Then, resorting the nabla Laplace transform, the rationality of the finite dimensional frequency distributed model approaching the infinite one is illuminated. Based on this, an original algorithm to estimate the parameters of the approximate model is proposed with the help of vector fitting method. Additionally, the applicable object is extended from a sum operator to a general system. Three numerical examples are performed to illustrate the applicability and flexibility of the introduced methodology.

Tables

Algorithm 1: The numerical approximation of fractional sum operator

Prior Information: the fractional order $\alpha$ and the discrete frequency points $s_l,\left(l=1,2,\cdots ,L\right)$

Input: the total iterations number $T$ and the degree of approximation model $N+1$ $\left(L>2N+2\right)$

Output: the poles $-{\omega }_k$ and the residues $c_i$ $\left(i=0,1,\cdots ,N\right)$

Step 1: Let the number of iterations $t=0$ and choose initial poles $-p_i$.

Step 2: Use (39) to calculate $y\left(s_l\right)$ and $\varphi \left(s_l\right)$.

Step 3: Use (41) to calculate $\hat{\theta }$, and get ${\mu }_i$ and ${\lambda }_i$.

Step 4: Use (35) to calculate the zeros $-{\omega }_i$ of $h\left(s\right)$ based on $p_i$ and ${\lambda }_i$.

Step 5: If $t\ge T$, take the zeros $-{\omega }_i$ obtained in Step 4 as the poles of ${\hat{S}}_{\alpha }\left(s\right)$

and use the least squares algorithm to find the residues $c_i$ of ${\hat{S}}_{\alpha }\left(s\right)$.

Otherwise, return to Step 2 and let $p_i={\omega }_i,t=t+1$.

Step 6: Finally, put ${\omega }_i$ and $c_i$ as the parameters of ${\hat{S}}_{\alpha }\left(s\right)$.

Table 2. Table 1: The error with different approximation order

$J$	$\alpha =0.1$	$\alpha =0.2$	$\alpha =0.3$	$\alpha =0.4$	$\alpha =0.5$	$\alpha =0.6$	$\alpha =0.7$	$\alpha =0.8$	$\alpha =0.9$
$N=5$	$0.6839$	$1.2615$	$1.6254$	$1.8341$	$1.9514$	$1.9669$	$1.8500$	$1.5517$	$1.0079$
$N=10$	$0.0431$	$0.0882$	$0.1409$	$0.1780$	$0.2140$	$0.2725$	$0.2497$	$0.3315$	$0.2095$
$N=15$	$0.0090$	$0.0261$	$0.0425$	$0.0997$	$0.1073$	$0.1302$	$0.1868$	$0.1949$	$0.1785$
$N=20$	$0.0057$	$0.0407$	$0.0284$	$0.0333$	$0.0610$	$0.0880$	$0.1611$	$0.1360$	$0.1511$

Table 3. Table 2: The error with different number of iterations

$J$	$T=3$	$T=6$	$T=9$	$T=11$	$T=12$	$T=15$	$T=16$	$T=18$	$T=21$
$N=5$	$10.7382$	$1.9514$	$1.9414$	$2.1431$	$2.2061$	$2.2960$	$2.3080$	$2.3207$	$2.3272$
$N=10$	$4.9554$	$0.2140$	$0.1003$	$0.1202$	$0.1312$	$0.1549$	$0.1597$	$0.1663$	$0.1712$
$N=15$	$4.0491$	$0.1073$	$0.0105$	$0.0052$	$0.0051$	$0.0069$	$0.0074$	$0.0082$	$0.0088$
$N=20$	$3.1467$	$0.0610$	$0.0053$	$0.0013$	$7.73e-4$	$2.67e-4$	$2.66e-4$	$3.23e-4$	$4.03e-4$

Equations

$${\textstyle{}_{a}\nabla_{k}^{-\alpha}f(k)\triangleq\sum\nolimits_{i=0}^{k-a-1}(-1)^{i}\bigl{(}\begin{smallmatrix}-\alpha\\ i\end{smallmatrix}\bigr{)}f(k-i),}$$

$${}_a{\nabla }_k^{\alpha }f\left(k\right)\triangleq {}_a{\nabla }_k^{\alpha -n}{\nabla }^{n}f\left(k\right),$$

$${\mathcal{N}}_a\{f(k)\}\triangleq {\sum }_{k=1}^{+\infty }(1-s{)}^{k-1}f(k+a),s\in \mathbb{C}.$$

$$f\left(a+\kappa \right)=\underset{s\to 1}{\lim }\frac{F\left(s\right)-{\sum }_{k=1}^{\kappa -1}{\left(1-s\right)}^{k-1}f\left(a+k\right)}{{\left(1-s\right)}^{\kappa -1}},$$

$${\textstyle\begin{array}[]{rl}f\left(k\right)=&{\mathscr{N}}_{a}^{-1}\{F\left(s\right)\}\\ \triangleq&\frac{1}{{2\pi\rm{j}}}\oint_{c}{F\left(s\right){{(1-s)}^{-k+a}}{\rm{d}}s},k\in{\mathbb{N}_{a+1}},\end{array}}$$

$$\displaystyle{\textstyle\begin{array}[]{l}\frac{1}{{2\pi{\rm{j}}}}\oint_{c}{F\left(s\right){{(1-s)}^{-k+a}}{\rm{d}}s}\\ =\frac{1}{{2\pi{\rm{j}}}}\oint_{c}{\sum\nolimits_{i=1}^{+\infty}{{{(1-s)}^{i-1}}f(i+a)}{{(1-s)}^{-k+a}}{\rm{d}}s}\\ =\frac{1}{{2\pi{\rm{j}}}}\sum\nolimits_{i=1}^{+\infty}{f(i+a)}\oint_{c}{{{(1-s)}^{i-k+a-1}}{\rm{d}}s}.\end{array}}$$

$$\displaystyle{\textstyle\begin{array}[]{l}\frac{1}{{2\pi{\rm{j}}}}\oint_{c}{F\left(s\right){{(1-s)}^{-k+a}}{\rm{d}}s}\\ =\frac{1}{{2\pi{\rm{j}}}}\sum\nolimits_{i=1}^{+\infty}{f(i+a)}\int_{-\pi}^{\pi}{{{(r{{\rm e}^{-{\rm{j}}\theta}})}^{i-k+a-1}}{\rm{j}}r{{\rm e}^{-{\rm{j}}\theta}}{\rm{d}}\theta}\\ =\frac{1}{{2\pi}}\sum\nolimits_{i=1}^{+\infty}{{r^{i-k+a}}}{f(i+a)}\int_{-\pi}^{\pi}{{{\rm e}^{-{\rm{j}}\theta\left({i-k+a}\right)}}{\rm{d}}\theta}.\end{array}}$$

$${\textstyle\int_{-\pi}^{\pi}{{{\rm e}^{-{\rm{j}}\theta\left({i-k+a}\right)}}{\rm{d}}\theta}=\left\{\begin{array}[]{rl}2\pi&,i=k-a,\\ 0&,i\neq k-a,\end{array}\right.}$$

$$\displaystyle{\textstyle\begin{array}[]{rl}\frac{1}{{2\pi{\rm{j}}}}\oint_{c}{F\left(s\right){{(1-s)}^{-k+a}}{\rm{d}}s}=&\frac{{{r^{0}}}}{{2\pi}}f(k-a+a)2\pi\\ =&f(k).\end{array}}$$

$${\mathcal{N}}_a^{-1}\{F\left(s\right)\}=\frac{1}{2\pi \mathrm{j}}{\oint }_{c}F\left(s\right){(1-s)}^{-k+a}\mathrm{d}s,$$

$${\textstyle\left\{\begin{array}[]{rl}\nabla z(\omega,k)=&-\omega z(\omega,k)+u(k),\\ y(k)=&\int_{0}^{+\infty}\mu_{\alpha}(\omega)z(\omega,k){\rm d}\omega,\end{array}\right.}$$

$${\textstyle\left\{\begin{array}[]{rl}\nabla z(\omega,k)=&-\omega z(\omega,k)+\mu_{\alpha}(\omega)u(k),\\ y(k)=&\int_{0}^{+\infty}z(\omega,k){\rm d}\omega,\end{array}\right.}$$

$${\textstyle\left\{\begin{array}[]{rl}\nabla z(\omega,k)=&-\omega z(\omega,k)+\sqrt{\mu_{\alpha}(\omega)}u(k),\\ y(k)=&\int_{0}^{+\infty}\sqrt{\mu_{\alpha}(\omega)}z(\omega,k){\rm d}\omega,\end{array}\right.}$$

$${\textstyle\left\{\begin{array}[]{rl}\nabla z(\omega,k)=&-\omega^{2}z(\omega,k)+u(k),\\ y(k)=&\int_{-\infty}^{+\infty}\omega\mu_{\alpha}(\omega^{2})z(\omega,k){\rm d}\omega,\end{array}\right.}$$

$${\textstyle\left\{\begin{array}[]{rl}\nabla z(\omega,k)=&-\omega^{2}z(\omega,k)+\mu_{\alpha}(\omega^{2})u(k),\\ y(k)=&\int_{-\infty}^{+\infty}\omega z(\omega,k){\rm d}\omega,\end{array}\right.}$$

$${\textstyle\left\{\begin{array}[]{rl}\nabla z(\omega,k)=&-\omega^{2}z(\omega,k)+\sqrt{\mu_{\alpha}(\omega^{2})}u(k),\\ y(k)=&\int_{-\infty}^{+\infty}\omega\sqrt{\mu_{\alpha}(\omega^{2})}z(\omega,k){\rm d}\omega.\end{array}\right.}$$

$${\textstyle\left\{{\begin{array}[]{rl}{\nabla}z\left(\omega,k\right)=&-\omega z\left(\omega,k\right)+u\left(k\right),\\ y\left(k\right)=&\int_{0}^{+\infty}{{\mu_{\alpha}}\left(\omega\right)z\left(\omega,k\right){\rm d}\omega},\end{array}}\right.}$$

$${\textstyle\left\{{\begin{array}[]{rl}{\nabla}z\left(\omega,k\right)=&-\omega z\left(\omega,k\right)+{\nabla}y\left(k\right),\\ u\left(k\right)=&\int_{0}^{+\infty}{{\mu_{1-\alpha}}\left(\omega\right)z\left(\omega,k\right){\rm d}\omega},\end{array}}\right.}$$

$$\displaystyle\left\{{\begin{array}[]{rl}{\nabla}z\left(\omega,k\right)=&-\omega z\left(\omega,k\right)+u\left(k\right),\\ \sigma\left(k\right)=&\int_{0}^{+\infty}{{\mu_{\alpha-n+1}}\left(\omega\right)z\left(\omega,k\right){\rm d}\omega},\\ {\nabla^{n-1}}y\left(k\right)=&\sigma\left(k\right),\end{array}}\right.$$

$$\frac{1}{s^{\alpha }}={\int }_0^{+\infty }\frac{{\mu }_{\alpha }(\omega )}{s+\omega }\mathrm{d}\omega ,$$

$${\textstyle\left\{\begin{array}[]{rl}\nabla z(\omega_{i},k)=&-\omega_{i}z(\omega_{i},k)+u(k),\\ y(k)=&\sum\nolimits_{i=0}^{N}c_{i}z(\omega_{i},k),\end{array}\right.}$$

$${\hat{S}}_{\alpha }(s)={\sum }_{i=0}^N\frac{c_i}{s+{\omega }_i}=\frac{G_{\alpha }}{s+{\omega }_0}{\prod }_{i=1}^N\frac{s+{\bar{\omega }}_i}{s+{\omega }_i}.$$

$${\textstyle{\rm arg}\hskip 2.84526pt\min\limits_{c_{i},\omega_{i}}\sum\nolimits_{l=1}^{L}\big{|}\frac{1}{s_{l}^{\alpha}}-\sum\nolimits_{i=0}^{N}\frac{c_{i}}{s_{l}+\omega_{i}}\big{|}^{2}.}$$

$${\textstyle\left\{\begin{array}[]{rl}h(s)&\triangleq\prod\nolimits_{i=0}^{N}{\frac{{s+{{\omega}_{i}}}}{{s+{p_{i}}}}}=\sum\nolimits_{i=0}^{N}\frac{\lambda_{i}}{s+p_{i}}+1,\\ H(s)&\triangleq\frac{1}{s+p_{0}}\prod\nolimits_{i=1}^{N}{\frac{{s+{\bar{\omega}_{i}}}}{{s+{p_{i}}}}}=\sum\nolimits_{i=0}^{N}\frac{\mu_{i}}{s+p_{i}}.\end{array}\right.}$$

$${\hat{S}}_{\alpha }(s)=\frac{H(s)}{h(s)}.$$

$$S_{\alpha }(s)=\frac{H(s)}{h(s)}=\frac{{\sum }_{i=0}^N\frac{{\mu }_i}{s+p_i}}{{\sum }_{i=0}^N\frac{{\lambda }_i}{s+p_i}+1},$$

$${\textstyle\sum\nolimits_{i=0}^{N}\frac{\mu_{i}}{s+p_{i}}-\big{(}\sum\nolimits_{i=0}^{N}\frac{\lambda_{i}}{s+p_{i}}\big{)}S_{\alpha}(s)=S_{\alpha}(s),}$$

$$y(s_l)={\varphi }^{\mathrm{T}}(s_l)\theta ,$$

$$\underset{\theta }{\min }{\sum }_{l=1}^L{\left[y(s_l)-\varphi (s_l)\theta \right]}^2,$$

$$\begin{array}[]{l}y({s_{l}})={[{\mathop{\rm Re}\nolimits}\left\{{{S_{\alpha}}({s_{l}})}\right\},{\mathop{\rm Im}\nolimits}\left\{{{S_{\alpha}}({s_{l}})}\right\}]^{\rm{T}}},\\ \theta={\left[{{\mu_{0}},\ldots,{\mu_{N}},{\lambda_{0}},\ldots,{\lambda_{N}}}\right]^{\rm{T}}},\\ \phi({s_{l}})=\left[{{\mathop{\rm Re}\nolimits}\left\{{\varphi({s_{l}})}\right\},{\mathop{\rm Im}\nolimits}\left\{{\varphi({s_{l}})}\right\}}\right],\\ \varphi({s_{l}})=\big{[}{\frac{1}{{{s_{l}}+{p_{0}}}},\ldots,\frac{1}{{{s_{l}}+{p_{N}}}},\frac{{-{S_{\alpha}}({s_{l}})}}{{{s_{l}}+{p_{0}}}},\ldots,\frac{{-{S_{\alpha}}({s_{l}})}}{{{s_{l}}+{p_{N}}}}}\big{]}^{\rm{T}}.\end{array}$$

Full text

1 Introduction

Fractional calculus is a natural generalization of classical integer order 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], controller synthesis [5, 6], and optimization algorithm [7, 8], 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 [9, 10, 11] and the references therein.

Despite fractional calculus is the extension of the traditional calculus, it has complex definition which brings it long memory characteristic [12]. Therefore, it is difficult to simulate fractional order systems and implement fractional order controllers [13]. From another perspective, fractional order systems are essentially infinite dimensional [14], which makes it difficult to obtain their analytical time–domain response especially for the nonlinear case under nonzero initial conditions. Thus, numerical approximation problem 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 [15, 16, 17, 18, 19, 20, 21].

It is worth mentioning that the most famous work in numerical approximation is the Oustaloup method [15]. The main idea of this method is to use a series of polylines with the slope of $20{\rm dB/dec}$ and $0{\rm dB/dec}$ to approximate the magnitude–frequency curve that is an elegant straight line with the slope of $20\alpha{\rm dB/dec},\alpha\in(0,1)$. 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. To reduce the dimension of the approximation model, the so-called fixed-pole approximation method was proposed in [18, 19]. In this method, the poles will be fixed for each integral order $\alpha\in(0,1)$, 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. Bearing this in mind, the frequency–domain identification based method was considered [20, 21]. However, all of them focused on the continuous time system and none of the discrete time system was considered.

With the rapid development and wide application of digital computers, discrete–time systems become increasingly popular. In this process, some pioneering work have been done on discrete fractional calculus and discrete fractional order systems [22, 23, 24, 25, 26]. 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 [27]. Maione proved that the Laguerre approximation to the Tustin fractional operator $s^{-\alpha}$ was stable and minimum-phase [28]. Fractional order Euler-Frobenius polynomials was defined to characterize the asymptotic sampling zeros for fractional systems as the sampling period tends to zero [29]. An original direct discretization method was designed to produce low integer order $z$ domain transfer functions via inverse fast Fourier transform [30]. Different from the discretization of continuous time fractional order systems, the numerical scheme in the matrix form is given for discrete fractional order systems directly [31]. 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 [32]. However, it is not convenient to select the interested frequency range and the approximate degree by this approach. With further research, it has found that this issue is fundamentally an approximation problem of fractional sum operation $\nabla^{-\alpha}$ in the frequency domain of nabla Laplace transform. Inspired by this, the purpose of this paper is to develops a new strategy that not only approximates fractional order systems with high accuracy, but also differs from traditional curve fitting methods. To complete this, there are many challenges. i) It is difficult to realize a system in the frequency domain. ii) The nonlinear parameter identification problem hinders the access of approximation model. iii) The conditions on nonzero initial instant and nonzero initial state bring more difficulty further.

Before end up this section, the main contributions of this paper are summarized as follows.

1. i)

The inverse nabla Laplace transform is defined and proven for the first time, which implies the nonzero initial instant problem. Afterwards, the method to calculate $f(k)$ from its nabla Laplace transform $F(s)$ is discussed in detail. 2. ii)

Inspired by the unilateral frequency distributed model, three bilateral equivalent frequency distributed models are derived for fractional sum equation. Afterwards, the equivalent frequency distributed models for three special fractional difference equations are deduced. 3. iii)

Based on the infinite dimensional model of nabla fractional dynamic system, the approximation problem is transformed into a nonlinear parameter identification problem in frequency domain, and the needed parameters are determined by vector fitting method [33]. 4. iv)

The proposed approximation strategy is applied to simulate nabla fractional dynamic systems with nonzero initial instant and nonzero initial state. Though this method is developed from the frequency domain, it could also be applied for the nonlinear case.

The remainder of this paper is organized as follows. Section II presents preliminaries, including some fundamental knowledge of nabla discrete fractional calculus and nabla Laplace transform. Section III proposes the numerical approximation method for fractional sum operators and nabla fractional order systems. To illustrate the validity of proposed approach, three illustrative examples are contained in Section IV. Finally, some conclusions are drawn in Section V.

2 Preliminaries

In this section, the fundamental knowledge for nabla fractional calculus and the preparatory problem statement [24] will be provided.

The $\alpha$-th nabla fractional sum of a function $f$: $\mathbb{N}_{a+1}\rightarrow\mathbb{R}$ can be defined by

[TABLE]

where $\alpha>0$, $k\in\mathbb{N}_{a+1}$, $\mathbb{N}_{a+1}\triangleq\{a+1,a+2,a+3,\cdots\}$, \bigl{(}\begin{smallmatrix}p\\ q\end{smallmatrix}\bigr{)}=\frac{{\rm{\Gamma}}(p+1)}{{\rm{\Gamma}}(q+1){\rm{\Gamma}}(p-q+1)} and ${\rm\Gamma}(\cdot)$ is the Gamma function.

With this definition, the nabla Caputo fractional difference is defined by

[TABLE]

where $n-1<\alpha<n$, $n\in\mathbb{N}_{+}$, and ${\nabla^{n}}$ represents the normal $n$-th backward difference ${\nabla^{n}}f\left(k\right)\triangleq\sum\nolimits_{i=0}^{n}{{{{\left({-1}\right)}^{i}}\left({\begin{smallmatrix}{n}\\ i\end{smallmatrix}}\right)}f\left({k-i}\right)}$.

The nabla Laplace transform of a function $f$: $\mathbb{N}_{a+1}\rightarrow\mathbb{R}$ is defined as [34]

[TABLE]

Due to the existence of the parameter $a$, nonzero initial instant problem can be solved accordingly. Assuming $F(s)=\mathscr{N}_{a}\{f(k)\}$, then one has [35]

[TABLE]

where $\kappa\in\mathbb{Z}_{+}$. Till now, it can be concluded that given $F(s)$, the value of $f(k),~{}k\in\mathbb{N}_{a+1}$ can be computed. However, as $k$ increases gradually, calculating the value of $f(k)$ one by one has proved difficult. Taking limit operation is not an easy task either. As a consequence, an alternative solution is expected.

Theorem 1

If $F(s)={{\mathscr{N}}_{a}}\left\{{f(k)}\right\}$ with $f:\mathbb{N}_{a+1}\to\mathbb{R}$ and $a\in\mathbb{R}$, then the inverse nabla Laplace transform can be expressed as

[TABLE]

where $c$ is a closed curve rotating around the point $(1,{\rm j}0)$ clockwise and it also locates in the convergent domain of $F(s)$.

Proof For the given curve $c$, then one has

[TABLE]

Define $s=1-r{{\rm e}^{-{\rm{j}}\theta}}$ with $r>0$, $\theta\in(-\pi,\pi)$ and keep the integration curve $c$ is inside the convergent region of $F\left(s\right)$. Then, ${\rm{d}}s={\rm{j}}r{{\rm e}^{-{\rm{j}}\theta}}{\rm{d}}\theta$ and

[TABLE]

With the help of the following fact

[TABLE]

one has

[TABLE]

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

[TABLE]

where $k\in{\mathbb{N}_{a+1}}$. This completes the proof. $\Box$

It is noteworthy that if the nabla Laplace transform $F(s)$ is given, $f(k)$ can be calculated via Theorem 1. Sometimes, it is difficult to solve such a contour integral problem. Especially, in many situations, $F(s)$ is not given previous or impossible to access directly. For example, $f(k)$ is the state or output of a nonlinear system. Thus, a more effective and practical approach is to be proposed.

With the help of the frequency distribution model theory in [36], the system ${}_{a}\nabla_{k}^{-\alpha}u(k)=y(k)$ with $0<\alpha<1$ can be expressed equivalently as the following unilateral frequency distributed models

[TABLE]

[TABLE]

[TABLE]

where $u(k)$ is the input, $y(k)$ is the output, $z\left(\omega,k\right)$ is the state, ${\mu_{\alpha}}\left(\omega\right)=\frac{{\sin\left({\alpha\pi}\right)}}{{{\omega^{\alpha}}\pi}}$ is the weight function and $z\left(\omega,a\right)=0$ is the initial value. It is noted that different locations of the weight function give different models, such as, the output form in (19), the input form in (20) and the balance form in (21).

Likewise, different equivalent bilateral frequency distributed model of ${}_{a}\nabla_{k}^{-\alpha}u(k)=y(k)$ with $0<\alpha<1$ can be obtained as follows.

[TABLE]

[TABLE]

[TABLE]

If the fractional sum equation ${}_{a}\nabla_{k}^{-\alpha}u(k)=y(k)$ is replaced by the difference one ${}_{a}\nabla_{k}^{\alpha}y(k)=u(k)$, the corresponding equivalent model in the output form can be expressed as

[TABLE]

where $\alpha\in(0,1)$ and $z\left(\omega,a\right)=\frac{\delta\left(\omega\right)}{{\mu_{\alpha}}\left(\omega\right)}x\left(a\right)$ which is different from the zero initial value in (19)-(24).

In (25), $u(k)$ is known and $y(k)$ is to be calculated. On the contrary, when $y(k)$ is known and $u(k)$ is to be calculated, then the formula $u(k)={}_{a}\nabla_{k}^{\alpha}y(k)$ can be equivalently expressed as

[TABLE]

where $\alpha\in(0,1)$ and $z\left(\omega,a\right)=0$.

In (25), the order $\alpha\in(0,1)$. If $\alpha\in(n-1,n)$, $n\in\mathbb{Z}_{+}$, the system ${}_{a}\nabla_{k}^{\alpha}y(k)=u(k)$ can be equivalently expressed as

[TABLE]

where $z\left(\omega,a\right)=\frac{\delta\left(\omega\right)}{{\mu_{\alpha}}\left(\omega\right)}\nabla^{n-1}x\left(a\right)$.

Note that all the developed frequency distributed models lay the groundwork for the numerical implementation for nabla fractional order systems. They are infinite dimensional and cannot be used directly. However, their transfer functions are derived from the following identical equation [20]

[TABLE]

with $\alpha\in(0,1)$, $s\in\mathbb{C}\backslash\mathbb{R}_{-}$, which provides the possible to solve the numerical simulation problem.

3 Numerical Approximation Scheme

In this section, an effective algorithm is developed to approximate the nabla fractional sum operator $S_{\alpha}(s)=\frac{1}{s^{\alpha}}$ via system identification technique. 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.

3.1 For nabla fractional sum operator

To simulate the nabla fractional order system with infinite dimensional characteristic, a practical solution is to discretize the continuous frequency range $\omega\in[0,+\infty)$ with finite distributed frequency points $\omega_{0},\omega_{1},\cdots,\omega_{N}$. From this, an approximate finite dimensional state space model can be expressed as

[TABLE]

where $c_{i}=\mu_{\alpha}(\omega_{i})(\omega_{i}-\omega_{i-1})$ is the weight value and $z(\omega_{i},k)\in\mathbb{R}$ is the approximate system state. Similarly, the transfer function of the above model (32) can be obtained as follows

[TABLE]

Therefore, the task of numerical approximation is to approximate $S_{\alpha}(s)$ with $\hat{S}_{\alpha}(s)$. It has been pointed by the reference [37] that if the conditions that $N\to+\infty$, $\omega_{0}\to 0$, $\omega_{N}\to+\infty$, $\mathop{\sup}\limits_{1\leq i\leq N}\left|{{\omega_{i}}-{\omega_{i-1}}}\right|\to 0$ are satisfied, then the system (33) would approximate the nabla fractional sum operator $S_{\alpha}(s)$ with arbitrary accuracy. In this case, the finite dimensional frequency distribution model (32) can be adopted to approximate the infinite dimensional frequency distribution model (19). Therefore, the problem of operator approximation has been transformed into a system identification problem: when the order $\alpha$ and $s_{l},l=1,2,\cdots,L\in\mathbb{Z}_{+}$ are known, how to determine the parameters $\omega_{i}$ and $c_{i},i=0,1,\cdots,N\in\mathbb{N}$, so as to minimize the error between $S_{\alpha}(s)$ and $\hat{S}_{\alpha}(s)$, that is

[TABLE]

At superficial glance, this is a nonlinear identification problem. To solve this problem, it is necessary to introduce two auxiliary transfer functions as

[TABLE]

Compared with (33) and (35), it can be found that the numerator of $h(s)$ is the denominator of $\hat{S}_{\alpha}(s)$, and the numerator of $H(s)$ is the numerator of $\hat{S}_{\alpha}(s)$. Therefore, as long as the numerators of two auxiliary functions are determined, $\hat{S}_{\alpha}(s)$ would be obtained. From this, one can obtain that

[TABLE]

For the purpose of obtaining the parameters $\omega_{i}$ and $c_{i}$, the following equation should be guaranteed

[TABLE]

which can be easily expanded as

[TABLE]

where $p_{k}$ is the known parameters when introducing the auxiliary transfer functions. In general, we choose the given discrete frequency points $s_{l}={\rm j\,}\zeta_{l}$ with $\zeta_{l}\in\mathbb{R}_{+}$, $l=1,2,\cdots,L$. Then, $S_{\alpha}(s)$ can be regarded as the frequency–domain response. Obviously, the unknown parameters in (38) are only $\mu_{i}$ and $\lambda_{i}$ which only exist in the numerators. Thus, the original nonlinear least squares problem in (34) is successfully transformed into the following problem

[TABLE]

whose least square solution satisfies

[TABLE]

where

[TABLE]

Taking what is the above mentioned into account, the estimated values of $\mu_{k}$ and $\lambda_{k}$ are acquired as

[TABLE]

where $\Phi(s)=\left[{\phi(s_{1})},{\phi(s_{2})},\cdots,{\phi(s_{L})}\right]$, $Y(s)=\left[y(s_{1}),\right.$ $\left.y(s_{2}),\cdots,y(s_{L})\right]^{\rm T}.$

Note that $S_{\alpha}(s)$ is an irrational transfer function with infinite degree while $\hat{S}_{\alpha}(s)$ is a rational transfer function with finite degree $N+1$. As a result, only a least square solution can be obtained if we want to approximate $S_{\alpha}(s)$ with $\hat{S}_{\alpha}(s)$. Because $p_{i}$ will affect the bandwidth of $h(s)$ and $H(s)$, the desired case is that $p_{i}=\omega_{i}$. It can be observed from (35) that when the initial poles $-p_{i}$ of the auxiliary function are chosen as the exact poles $-\omega_{i}$ one has $h(s)=1$, $H(s)=S_{\alpha}(s)$. In this case, the best approximate performance can be achieved. 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 (35) also suggests that after the above calculation, the zeros $-\omega_{i}$ of $h(s)$ can be used as the poles $-p_{i}$ in the next cycle. In accordance with this iterative relationship, more accurate results would be obtained.

Remark 1

It should be noted that the backwards difference $\nabla z(\omega,k)=z(\omega,k)-z(\omega,k-1)$. The approximate model (32) cannot be transformed into the existing case \left\{\begin{array}[]{rl}x\left({k+1}\right)=&Ax\left(k\right)+Bu(k)\\ y(k)=&Cx\left(k\right)+Du(k)\end{array}\right., and therefore 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.

3.2 Treatment of conjugate complex numbers

Compared with the recursive schemes [15, 16, 17, 18, 19], one of the advantages of the proposed 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 $S_{\alpha}(s)$ must be strictly regular, the chosen complex poles must be complex conjugate. Besides, considering 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 (35), assume that the $(i+1)$-th and $(i+2)$-th poles are the complex conjugate, that is

[TABLE]

In this basis, the corresponding residues satisfy the following conjugate relationship

[TABLE]

Then, the corresponding elements in vector $\varphi(s_{l})$ and $\theta$ can be rewritten as

[TABLE]

and

[TABLE]

In this way, the original least squares solution $\hat{\theta}$ can be found in the real space and the conjugacy can be guaranteed firmly. The specific computing formula is also (41). If the estimated parameter pairs corresponding to $p_{i}$ and $p_{i+1}$ are $\theta_{i+1},\theta_{i+2}$ and $\theta_{2i+2},\theta_{2i+3}$, then the original complex parameters $(\mu_{i},\mu_{i+1})$ and $(\lambda_{i},\lambda_{i+1})$ can be recovered via the relationship in equation (43) and equation (45).

Based on the above operation, the problem with complex poles has been solved successfully. To sum up, the entire calculation process proposed above can be described as Algorithm 1.

3.3 For nabla fractional dynamic systems

After approximating the nabla fractional sum operator $\frac{1}{s^{\alpha}}$, the presented results can be extended to nabla discrete fractional dynamic systems

[TABLE]

where the parameters satisfy

[TABLE]

By applying the frequency distribution model theory, the finite dimensional pseudo state space system in (46) can be rewritten as an infinite dimensional exact state space model

[TABLE]

where ${\bf z}(\omega,k)=[~{}z_{1}(\omega,k),~{}z_{2}(\omega,k),~{}\cdots,~{}z_{n}(\omega,k)~{}]^{\rm T}$, ${\bf\mu}_{\bf\alpha}(\omega)={\rm diag}\{~{}\mu_{1}(\omega),~{}\mu_{2}(\omega),~{}\cdots,~{}\mu_{n}(\omega)~{}\}$, $\mu_{i}(\omega)=\frac{\sin(\alpha_{i}\pi)}{\omega^{\alpha_{i}}\pi}$ and $z_{j}\left(\omega,a\right)=\frac{\delta\left(\omega\right)}{{\mu_{\alpha_{j}}}\left(\omega\right)}x_{j}\left(a\right)$, $j=1,2,\cdots,n$.

As mentioned earlier, such an infinite dimensional model can not be used directly in practice. Therefore, the system (46) needs to be approximated by a finite dimensional model as follows

[TABLE]

where $M_{A}={\rm diag}\{A_{0},A_{1},\ldots,A_{N}\}$, $M_{B}=[B_{0},B_{1},$ $\ldots,B_{N}]^{\rm T}$, $M_{C}=[C_{0},C_{1},$ $\ldots,C_{N}]$, $A_{i}=-{\rm diag}\{\omega_{1,i},$ $\omega_{2,i},\ldots,\omega_{n,i}\}$, $B_{i}=I_{r}$, $C_{i}={\rm diag}\{c_{1,i},c_{2,i},\ldots,c_{n,i}\}$, the parameters $\omega_{j,i}$ and $c_{j,i}$ are generated by the algorithm in previous section for $G_{\alpha_{j}}(s)$.

The above discussion in (54) does not consider the case of nonzero initial state. When it comes to this situation, just consider how to reasonably assign the initial pseudo state $x(a)$ to the real state $z(\omega,a)$. For the discussed Caputo definition, $z(\omega,a)$ is completely distributed on the frequency point $\omega=0$ and therefore the initial state can be configured as follows

[TABLE]

Notably, the value of $\omega_{i}$ are obtained from identification and $\omega_{0}=0$ is no longer guaranteed. To achieve this, another approximation model should be assumed i.e., $\hat{S}_{\alpha}(s)=\frac{c_{0}}{s}+\sum\nolimits_{i=1}^{N}\frac{c_{i}}{s+\omega_{i}}=\frac{G_{\alpha}}{s}\prod\nolimits_{i=1}^{N}{\frac{{s+{{\bar{\omega}}_{i}}}}{{s+{\omega_{i}}}}}$. In other words, we use ${G_{\alpha}}\prod\nolimits_{i=1}^{N}{\frac{{s+{{\bar{\omega}}_{i}}}}{{s+{\omega_{i}}}}}$ to approximate $s^{1-\alpha}$ and then the zero pole can be guaranteed naturally.

Remark 2

Just as the previous discussion, the discrete frequency points $s_{l}$ are selected from the positive imaginary axis. For the continuous time case, the imaginary axis is the boundary to distinguish between stable region and unstable region. This selection strategy will maintain the stability well after approximating. However, for the discrete time case, the imaginary axis is only a curve in stable region and the actual boundary to distinguish between stable region and unstable region is a circle, $(x-1)^{2}+y^{2}=1$. Therefore, try to choose $s_{l}$ in this circle is our ongoing research.

Remark 3

If the adopted Caputo fractional difference is changed by the Riemann–Liouville fractional difference or the Grünwald–Letnikov fractional difference, only the initial value $z\left(\omega,a\right)$ should be replaced. If the nabla case is replaced by the delta case, similar infinite dimensional models can be derived. In this case, the numerical approximation in frequency domain can be performed again.

4 Simulation Study

In this section, three numerical examples will be presented to demonstrate the effectiveness of the proposed method.

Example 1: For fractional sum operator

The range of discrete frequency points is set as $[\omega_{l},\omega_{h}]=[0.001,1000]$ in this section. To illustrate the effectiveness of approximate method, the following approximate error is introduced.

[TABLE]

With different orders $\alpha=0.1,0.2,\cdots,0.9$, 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 $N$ increases.

To illustrate the influence of the iteration number further, the approximation for $\frac{1}{s^{0.5}}$ is discussed with results shown in Table 2. From a macroscopic point of view, as the number of iteration increased, the approximation accuracy also increases. 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 $T$ around the chosen approximation degree $N+1$. Based on this discovery, it is good enough for the performance of approximation to set the number of iteration in the subsequent simulation as $L=N$, or slightly less than $N$. For convenience, the numbers of iteration in the following examples are set to $T=8$ 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 $S_{\alpha}(s)$ and $\hat{S}_{\alpha}(s)$. To reflect the approximation performance of the proposed approach more intuitively, the time–domain response is considered. Notably, consider $\alpha=0.5$, $a=5$, $N=20$ and $y(a)=1$ in this and the following examples. Herein, with the following input (see Fig. 1)

[TABLE]

one can get the output of system ${}_{a}\nabla_{k}^{-\alpha}u(k)=y(k)$ shown in Fig. 1.

Fig. 2 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 approximate error $\varepsilon(k)$ in Fig. 3. The analytical solution of $y(k)$ can be obtained by utilising the definition of fractional sum introduced in Section 2.

Example 2: For the linear system

Consider a system

[TABLE]

where the system input $u(k)=5\sin(0.2\pi k)$ is shown as Fig. 4. With the similar method in Example 1, the analytical output and the simulated output are displayed in Fig. 5. The error between the two outputs is represented in Fig. 6. The relative error is inferior to $10^{-5}$, which is quite acceptable.

Example 3: For the nonlinear system

Consider a nonlinear fractional dynamic system

[TABLE]

where the input $u(k)$ is a sawtooth wave signal shown as Fig. 7. The period is 5 and the amplitude is 5. In general, it is difficult to get the analytical output of a nonlinear system. To verify the effectiveness of the established approximate method, the nonlinear term is selected regarding to $y(k-1)$ not $y(k)$. The outputs $y(k)$ solved by the definition and the approximate method are displayed in Fig. 8. The approximate error is provided in Fig. 9.

It is shown that the approximate error is quite small, which illustrates that the developed scheme is very effective and efficient and the approximate results can reflect the dynamic performance of the original system to a large degree. Without loss of generality, the discussed method can be adopted in the simulation of the general nonlinear case and other different nabla fractional order systems.

5 Conclusions

In this paper, the numerical approximation has been investigated for nabla fractional dynamic systems with nonzero initial instant and nonzero initial state. According to the frequency distributed model theory, the exact infinite dimensional model of fractional sum operator is recalled and then the approximation problem is transformed into a nonlinear identification problem for the first time. By applying the vector fitting approach in frequency domain, an innovative and effective numerical approximation is developed. 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.