This work is licensed under the Creative Commons Attribution 4.0 International License.
Introduction
For a long time, we have focused on the mathematical theory of fractional calculus (FOC). It was not until the mid and late 20th century that FOC theory was introduced and applied in some fields such as electrochemistry, rheology, signal processing, and power transmission theory. Recently, the research of fractional order systems has aroused the interest of many scholars [1]. In the field of automatic control, a new branch of fractional control has emerged.
The research of fractional order control should cover at least three aspects: (1) The establishment of a fractional-order system (FOS) model and its analysis based on describing the fractional-order object more accurately and concisely. (2) The fractional-order control strategy is selected to obtain better control performance. (3) Apply fractional arithmetic to process signals, data, etc. The above work is based on the calculation of FOC and FOS.
Based on the Tustin transformation theory and the analysis of the characteristics of the discrete generating function formula used for fractional operators, and using the Maclaurin series expansion of the binomial power function, this paper proposes an improved method based on the power series expansion (PSE) and the Tustin transformation. Fractional calculation method [2]. We apply it to the solution of the linear fractional system and derive the corresponding recurrence formula. The performance of the algorithm's operation speed and accuracy is analyzed through specific example simulations and comprehensive comparisons with several conventional methods. The algorithm proposed in this paper is suitable for the fractional operation of arbitrary signals, including discrete data sequences whose mathematical model is unknown and the solution of linear systems.
Algorithm proposed
Analysis of commonly used algorithms and their limitations
The essence of the numerical method is to perform a discrete approximation of FDEs to obtain corresponding approximate numerical solutions [3]. The core problem is the discretization of the FOC operator. The basic idea of discretization is to use grid function f(nh) and generating function ω(ξ−1) to approximate the FOC of function f(t). In control theory, the sampling period T can be used to replace the lattice length h, z instead, and the function f(t) can be transformed into z the transformation of ξ the sequence f(nh). This realizes the discretization approximation. We will write it as
{D^{ \pm a}}(z) = (\omega ({z^{ - 1}}{))^{ \pm a}} = {s^{ \pm a}}
Corresponding to different methods, different generating functions can be obtained, the more commonly used ones are:
Since the corresponding generating function is irrational, it is not convenient to be directly used for solving. We can approximate it rationally [4]. Commonly used methods are PSE and continued fraction expansion (CFE). The combination of the above-mentioned different methods can get different forms of FOC operator discrete algorithm (PSE-Tustin, CFE-Tustin). Of course, no matter whether it is the PSE method or the CFE method, the solution expression obtained is still an infinite number of items. In its practical application, it can be considered to approximate its finite term, that is, the short memory method (SMP).
The research-based on traditional integer-order systems is relatively mature, and it is a basic research idea and method to approximate FOSs with integer-order systems [5]. So the problem is transformed into using standard integer-order operators to approximate fractional-order operators. This also overcomes the difficulty of not allowing direct fractional operator operations to calculate and simulate existing software. The fractional operator Sa in the complex frequency domain is generally an irrational number. The irrational transfer function can be transformed into a rational transfer function by approximating it with a rational function. Because it is easy to get the frequency domain characteristics of the FOS by the complex frequency domain calculus operator, we can approximate Sa the continuous domain integer order transfer function structure operator through the frequency domain fitting method. The approximation process can also be achieved employing stable optimal rational function fitting.
In the research, it is found that the integral transformation formula of the analytical algorithm is complicated, the calculation amount is large, and the dimension of the state space method is often too high. Models after rational approximation are often inconvenient for the solution and further analysis because of their high order and complexity. The existing analytical algorithms are all aimed at special forms of equations, but many equations often do not have analytical solutions [6]. On the other hand, the application premise of the analytical method is that the prototype function is known, and the problem of calculus for the known data is often encountered in practical engineering applications. Such as system parameter identification, signal analysis, etc. When applying PSE or CFE to the Eqs. (2) to (4) to realize the discretization approximation, it is necessary to call the related functions in Matlab or Maple. The calculation speed is very slow, and the resulting form cannot be directly called in programming. But modern control is based on the computer as the main realization tool. The meaning of its ‘realization’ is not limited to the ‘realization’ of its analysis and design process, but more importantly, it is put into actual operation. There are obvious obstacles to the above methods.
Tustin transformation
It is well known that the mapping relationship between the s plane in the continuous domain and the z plane in the discrete domain is z = exp(Ts). The left half-plane of z corresponds to the entire unit circle of the plane z. If Da f (t) is regarded as the response of the signal described by the function f(t), the Euler backward difference method only maps the left half-plane of z to the circle with the center of (0.5, 0) and the radius of 0.5 in the z plane [7]. This will inevitably lead to response distortion and error. Equation (3) Tustin transform will also be distorted but smaller than Eq. (2). The Al-Alaoui method is a compound operator of the former two, and the distortion is somewhere between. To reduce the error, we need to adopt a smaller sampling period T. This will inevitably increase the amount of calculation.
In the research based on the numerical example simulation, it is found that the Tustin method and Al-Alaoui method with the same number of expansion terms or orders are better than the Euler method in terms of accuracy no matter whether PSE or CFE is used. The Tustin method is slightly better than the Al-Alaoui method [8]. Therefore, from the perspective of response distortion or error, Tustin transform is a better choice.
Convergence and recursive realization
The convergence of the algorithm is one of the important issues that must be considered. Due to the possible limitations of the algorithm itself, when it is applied to the approximation of fractional-order operators, the system's output will diverge due to improper selection of the order (Figure 1). The fundamental reason is that the algorithm cannot guarantee convergence, or the convergence is not good [9]. Theoretically, because FOC has a unique attenuation memory characteristic, it must be shown as the uniform attenuation of the weight coefficients of the power terms of z−1 when calculating. However, they are not all well-reflected when using Matlab or Maple to expand and approximate Eqs. (2) to (4) as a whole. Taking α = 0.5 as an example, when the Taylor function in Matlab is called to expand the power series of Eq. (3) at z−1 = 0, the coefficients of the first 50 terms of z−1 are shown in Figure 2. It can be seen that the coefficient oscillates and the jump degree is large, and when the finite term is approximated, it will inevitably lose the information with a large weight. But the combination of formula (2) and PSE is an exception (Figure 3), and the convergence speed is very fast. Another significant advantage of this combination is that the transformation formula is a binomial power function. The weight coefficients of the expansion of the power series can be obtained quickly and conveniently. This is difficult to achieve in other combinations.
Fig. 1
Output divergence of CFE-Euler simulation model
Fig. 2
Maclaurin series expansion coefficient of Tustin transform
Fig. 3
Maclaurin series expansion coefficient of binomial power function
Tustin transform formula (3) is a fractional power function, and its numerator and denominator have the same structure as formula (2). Suppose the numerator and denominator are approximated as polynomials, and the two are regarded as the numerator and denominator of the pulse transfer function of a filter. In that case, this filter is a discrete FOC filter [10]. Furthermore, we can use the Maclaurin series expansion of the binomial power function to ensure convergence, and the weight coefficient of the expansion term can be quickly and recursively obtained.
We carry out the Maclaurin series expansion of the numerator and denominator of the function
f(x) = {\left( {(1 - x)/(1 + x)} \right)^{{1 \over 2}}}
of the structure of formula (3). The comparison between the approximate result and the theoretical curve of the original function is shown in Figure 4. It can be seen that a fairly high fitting accuracy can be achieved by taking 20 expansions.
Based on the above analysis, this paper proposes an improved FOC algorithm based on PSE and Tustin transform and applies it to the solution of linear fractional systems.
A recursive algorithm for FOC
We assume that the FOC of function f(t) is Da f (t). Where a is the order. Z transformation is denoted as F(z), DF(z) respectively. Considering Eqs. (1) and (3), we can get
{{DF(z)} \over {F(z)}} = {\left( {\omega \left( {{z^{ - 1}}} \right)} \right)^a} = {\left( {{2 \over T}\left( {{{1 - {z^{ - 1}}} \over {1 + {z^{ - 1}}}}} \right)} \right)^a} = {\left( {{2 \over T}} \right)^a}{{{{(1 - {z^{ - 1}})}^a}} \over {{{\left( {1 + {z^{ - 1}}} \right)}^a}}}
But
{(1 + {z^{ - 1}})^a}DF(z) = {\left( {{2 \over T}} \right)^a}{(1 - {z^{ - 1}})^a}F(z)
We carry out the Maclaurin series expansion of the binomial power function (1 − z−1)a (take n terms to approximate) and write it in matrix form
{\left( {1 - {z^{ - 1}}} \right)^a} \approx \sum\limits_{j = 0}^n {( - 1)^j}\left( {\matrix{ a \cr j \cr } } \right){z^{ - j}} = \left[ {{p_0},{p_1}, \cdots , \cdots } \right]{\left[ {1,{z^{ - 1}} \cdots {z^{ - j}} \cdots } \right]^T} = PZ\left( {\matrix{ a \cr j \cr } } \right)
is the binomial coefficient. P is the weight coefficient matrix. It is not difficult to get the recurrence relationship of the element
{p_j} = ( - {1)^j}\left( {\matrix{ a \cr j \cr } } \right)
of P by mathematical induction
{p_0} = 1,{p_j} = \left( {1 - {{a + 1} \over j}} \right){p_{j - 1}}
The same can be obtained
\matrix{ {{{\left( {1 - {z^{ - 1}}} \right)}^a} \approx \left[ {{q_0},{q_1}, \cdots ,{q_j}, \cdots } \right]} \hfill \cr {\left[ {1,{z^{ - 1}} \cdots {z^{ - j}} \cdots } \right] = QZ} \hfill \cr } {q_0} = 1,\quad {q_j} = \left( { - 1 + {{a + 1} \over j}} \right){q_{j - 1}}
We assume that the discrete sequences of f(t) and Daf (t) are f(kT) and Df (kT). Abbreviated as fk and dfk. k = 0, 1, 2,⋯ ,n. We substitute formula (7) and formula (9) into formula (6) to obtain d fk + q1d fk−1 + ⋯ + qk−1d f1 + qkd f0 = (2/T)a (fk + p1fk−1 + ⋯ + pk−1f1 + pk f0). Then
d{f_k} = (2/T{)^a}\left( {{f_k} + {p_1}{f_{k - 1}} + \cdots + {p_{k - 1}}{f_1} + {p_k}{f_0}} \right) - \left( {{q_1}d{f_{k - 1}} + \cdots + {q_{k - 1}}d{f_1} + {q_k}d{f_0}} \right)Formula (11) is the FOC solution formula of f(t). The specific recurrence steps are as follows.
Step 1: Discrete f(t) to obtain the sequence fk of f(t). And get matrix P, Q by recursive formula (8) and formula (10).
Step 2: The initial value of the parameter P, Q and d f0 is generally controlled as d f0 = 0.
Step 3: Take k = 1 and find d f1. d f1 = (2/T)a(f1 + p1f0) − q1d f0.
Step 4: Let k = 2 and find d f2: d f2 = (2/T)a (f2 + p1f1 + p2f0) − (q1df1 + q2df0).
Step 5: df3, df4, df5,⋯ ,dfn can be obtained by formula (11).
Solving fractional linear systems
First consider the general form of the differential equation of a special fractional-order linear system:
{a_1}{D^{a1}}y(t) + {a_2}{D^{a2}}y(t) + \cdots + {a_n}{D^{an}}y(t) = u(t)u(t) can be formed by a certain function and its fractional derivative. Considering formula (6), formula (7), formula (9), and formula (11), the approximate value daiyk of Daiy(t) can be obtained. We rewrite it as
{d^{{a_i}}}y(t) \approx {d^{{a_i}}}{y_k} = {\left( {{2 \over T}} \right)^{{a_i}}}\left[ {{y_k} + \sum\limits_{j = 1}^k p{j^{({a_i})}}{y_{k - j}}} \right] - \sum\limits_{j = 1}^k q{j^{({a_i})}}{d^{{a_i}}}{y_{k - j}}pj(ai) and qj(ai) can be directly recursively obtained by formula (8) and formula (10). We substitute it into Eq. (12) and arrange it as
\sum\limits_{i = 1}^n {a_i}\left\{ {{{\left( {{2 \over T}} \right)}^{{a_i}}}\left[ {{y_k} + \sum\limits_{j = 1}^k p{j^{({a_i})}}{y_{k - j}}} \right] - \sum\limits_{j = 1}^k q{j^{({a_i})}}{d^{{a_i}}}{y_{k - j}}} \right\} = {u_k}
. Then
{y_k} = {1 \over {\sum\limits_{i = 1}^n {a_i}{{\left( {{2 \over T}} \right)}^{{a_i}}}}}\left[ {{u_k} - \sum\limits_{i = 1}^n {a_i}{{\left( {{2 \over T}} \right)}^{{a_i}}}\sum\limits_{j = 1}^k q{j^{({a_i})}}{y_{k - j}} + \sum\limits_{i = 1}^n {a_i}\sum\limits_{j = 1}^k q{j^{({a_i})}}d{j^{({a_i})}}{y_{k - j}}} \right]Equations (13) and (14) are recursive algorithm formulas that constitute the output of the fractional linear system. Specific steps are as follows.
Step 1: Calculate the sequence uk of the input signal u(t) by the FOC recursive algorithm given in this article.
Step 2: We take the initial value: y0 = 0, daiy0 = 0.
Step 3: We make k = 1. Obtain y1 from Eq. (14). Substitute formula (13) to obtain the series daiy1.
Step 4: Let k = 2, use the previous result to obtain y2 from Eq. (14), and substitute Eq. (13) to get the series daiy2.
Step 5: Use the results of the previous steps to calculate y3, y4, y5,⋯yk by analogy.
Simulation analysis of algorithm accuracy
First, to verify the effectiveness of the FOC recursive solution algorithm proposed in this paper, here is a simple simulation curve of the 0.5-order derivative, 0.5-order integral, and the first-order integral of the sine function f(t) = sint (Figure 5).
Fig. 5
FOC curve of function sint. FOC, fractional calculus
It can be seen that the algorithm curve in this paper is the same as the result of the conventional PSE expansion approximation of the Tustin transform operator as a whole. This is almost the same as the curve obtained by the most direct and simple G-L definition discrete method [11]. The G-L method is almost identical to the theory obtained by applying the Cauchy definition to remove the area around t = 0+. Therefore, the results in Figure 5 show that the algorithm in this paper is feasible.
The unit step response is obtained through a specific FOS model. Based on the accuracy and speed, a comprehensive analysis and comparison between the method in this paper and the commonly used algorithms are made. The calculation and simulation tool is Matlab7. Let the FDEs of the simulation example be
{D^2}y(t) + 0.5{D^{0.5}}y(t) + y(t) = {D^{0.5}}u(t) + u(t)
For a class of FOS transfer function models that partial fractions can expand, the problem of obtaining theoretical numerical solutions has been well solved. We apply it to this example to easily find the analytical solution of its unit step output as
{y_{step}}(t) = {L^{ - 1}}\left[ {G(s)U(s)} \right] = {L^{ - 1}}\left[ {\sum\limits_{k = 1}^4 {{{c_k}{s^{ - 1}}} \over {{s^{0.5}} - {\lambda _k}}}} \right] = \sum\limits_{k = 1}^4 {c_k}{t^{0.5}}{E_{0.5,0.5 + 1}}\left( {{\lambda _k}{t^{0.5}}} \right)ck, λk is the pole and residue of G (s0.5) respectively. Take the error function
E = \left\| {y\mathop {step}\limits^* - ystep} \right\|1/\left\| {y\mathop {step}\limits^* } \right\|1y\mathop {step}\limits^*
is the comparison reference value, which is calculated by Eq. (16). ystep is the compared value. Considering that
y\mathop {step}\limits^*
reaches the steady-state value in about 20s, the length of the comparison data is 200 (step T = 0.1s). The precision index function is J = E*/E. E* is the reference value. Here we take the E value of the full memory G-L definition method (T = 0.1s). The index value has standard unity meaning. The simulation curves of each algorithm are shown in Figures 6–8. Table 1 shows the algorithm error and accuracy index values.
Fig. 6
Comparison of algorithm curve and theoretical solution curve in this paper
Fig. 7
The comparison curve between the full memory algorithm in this paper and the G-L definition @@LISAN method
Fig. 8
Comparison curve between the algorithm in this paper and the conventional PSE-Tustin method. PSE, power series expansion
Algorithm error and accuracy indicators
Algorithm name
Error
Accuracy index
This paper algorithm (memory length: N)
N = 20
0.0649
13.3
N = 50
0.0301
28.6
N = 100
0.0161
53.5
Full memory
0.012
370.1
G-L definition discrete method (full memory)
T = 0.1s
0.862
1
T = 0.05s
0.0557
15.4
T = 0.01s
0.0194
44.5
T = 0.001s
0.0052
165.8
Tustin PSE overall expansion
50
0.0683
12.6
100
0.0397
21.7
200
0.0123
70.1
50
0.0575
15
Al-Alaoui PSE overall expansion
100
0.0543
15.9
200
0.0536
16.1
PSE, power series expansion
The error is larger for step T = 0.1s (E = 0.862), and the smaller T is, the closer to the theoretical curve. Although high enough accuracy can be achieved when T = 0.001s is used, it is very difficult and too slow in terms of engineering realization. Considering the actual application of the project, the sampling period in the following simulations is selected to be easy to implement in the project.
Figure 6 shows that using the recursive discrete approximation algorithm given in this article to obtain the system's step response can obtain a relatively satisfactory result. The response curve of the memory length N = 100 almost coincides with the theoretical solution curve. Research shows (see Table 1) that the algorithm in this paper with an overall error of 10-2 orders of magnitude and full memory can achieve fairly high accuracy. The accuracy when N = 20 is equivalent to the result of the smaller step size T = 0.05s of the full memory G-L method. This verifies the conclusion that the response distortion of the Tustin discrete operator mentioned in the previous analysis is smaller than the Euler discrete operator. But when N > 100, the effect of accuracy improvement is not obvious [12]. This shows and verifies that the coefficients of the binomial power function can converge quickly. The above-mentioned excellent characteristics can be more clearly reflected in the comparison curve between the full memory algorithm in this paper and the G-L definition discrete method shown in Figure 7. The algorithm in this paper quickly converges to the theoretical value at the initial moment, and the accuracy advantage is also very significant.
It can be seen from Figure 8 that the initial response of the comparison curve between the algorithm in this paper and the conventional PSE-Tustin method is not much different. After 5 s, the difference is obvious, but this phenomenon is not reflected in Figure 5. The reason may be that different FOC discrete approximation methods exhibit different accuracy due to the difference in signal behavior. For this reason, we apply the fractional-order Al-Alaoui discrete operator to this example using the conventional PSE method and compare and analyze the results to show (Figure 9 and Table 1). The response curve shape of PSE-Al-Alaoui is similar to that of PSE-Tustin. In addition, the difference between the number of expansion terms of PSE-Tustin and the obvious influence on the accuracy is that the number of expansion terms has a less obvious influence on the accuracy. This shows that the convergence of the algorithm is better than that of PSE-Tustin. Overall, the results of the algorithm in this paper are the best. The PSE-Tustin method is better than the PSE-Al-Alaoui method.
Fig. 9
Step response of CFE-Al-Alaoui discrete method. CFE, continued fraction expansion
Concluding remarks
This paper proposes an improved PSE-Tustin recursive algorithm. This algorithm is suitable for the fractional operation of arbitrary signals, including discrete data sequences whose mathematical models are unknown and the solution of linear systems. A simulation example verifies the effectiveness of the algorithm. We show the characteristics and advantages of the algorithm in this paper by comparing the accuracy with several commonly used fractional-order algorithms.
