Huimin Zhao^{1} , Wu Deng^{2} , Xinhua Yang^{3} , Yu Xue^{4}
^{1, 2, 3}Software Institute, Dalian Jiaotong University, Dalian 116028, China
^{1, 2}The State Key Laboratory of Mechanical Transmissions, Chongqing University, Chongqing 400044, China
^{1, 2}Traction Power State Key Laboratory of Southwest Jiaotong University, Chengdu 610031, China
^{1, 2, 3}Dalian Key Laboratory of Welded Structures and Its Intelligent Manufacturing Technology(IMT) of Rail Transportation Equipment, Dalian Jiaotong University, Dalian 116028, China
^{4}Nanjing University of Information Science and Technology, Nanjing, 210044, China
^{2}Corresponding author
Journal of Vibroengineering, Vol. 18, Issue 5, 2016, p. 34123424.
https://doi.org/10.21595/jve.2016.16966
Received 12 March 2016; received in revised form 16 June 2016; accepted 14 July 2016; published 15 August 2016
Copyright © 2016 JVE International Ltd.
This paper presents a rational approximation method for fractional calculus operators ${s}^{\pm \mathrm{\alpha}}$in the given frequency range and the error, which is based on the best rational approximation definition. The fractional integral operator is selected as an example to describe the construction of the rational approximation functions. An application case (${s}^{0.3}$) is used to illustrate the effectiveness of the proposed method. The obtained approximation function in the frequency domain is a best rational approximation function, which can further improve the accuracy of the approximation without increasing the orders. On the basis of the presented rational approximation method, a rational approximation equation of a fractionalorder PID controller is obtained. Finally, the method for analyzing the optimization and frequency characteristics of the fractionalorder controller is implemented to demonstrate the good frequency characteristic and best structure. The results from theoretical analysis and experimental verification show that the proposed method provides a new design idea for the effective application of the fractionalorder PID controller in engineering.
Keywords: best rational approximation definition, fractional calculus operators, fractionalorder PID controller, model optimization.
In recent years, fractionalorder systems and controllers have attracted much attention [1]. Fractionalorder systems, as a generalization of classical integerorder systems, can be defined by the mentioned integral and derivative operators [25]. Compared with integerorder systems, fractionalorder systems can better describe the dynamics of real systems, such as electrodeelectrolyte polarization, electromagnetic waves, diffusion equations and so on [615]. Because the conventional proportionintegralderivative(PID) controller is the most widely used control type in industry, the fractionalorder PID controller has been considerably researched [16, 17]. The fractionalorder PID controllers, as a generalization of traditional integerorder PID controllers, are used to enhance the performance of control systems.
The differential equation of a fractionalorder PID controller is described as:
where ${K}_{p}$, ${K}_{i}$ and ${K}_{d}$ represent the proportion, integral, and derivative respectively.
The Eq. (1) is taken the Laplace transform, and zero initial conditions are assumed, thus the following equation is obtained:
where $s$ is the variable of Laplace. When $\alpha $ is 1 and $\beta $ is 1, the PID is a particular case.
As is known, when the integerorder systems have a limited memory, the fractionalorder systems have an unlimited memory. There mainly has the frequency domain based on the RiemannLiouville definition and GrunwaldLetnikov definition and the time domain based on Caputo definition. Tseng [18] investigated the design of a fractionalorder FIR differentiator. Chen et al. [19] presented a direct recursive discretization of the Tustin operator and direct discretization method by using the AlAlaoui operator. Chen and Vinagre [20] presented a new infinite impulse response type digital fractionalorder differentiator by using a new family of firstorder digital differentiators. Podlubny et al. [21] proposed several possible analog methods for obtaining suitable rational approximations. Krishan [22] reviewed the analysis, design and applications of analog and digital differentiators and integrators of fractionalorder and proposed time and frequency domain analysis method to realize a fractance device. Oustaloup et al. [23] conducted the synthesis of differentiators with integer, noninteger, real or complex orders, and the action is limited to any given frequency bandwidth. Charef [24] presented an analogue realization of fractionalorder integrator, differentiator and fractional PID controllers. Romero et al. [25] discussed the use of Chebyshev polynomial theory to achieve accurate discretetime approximations to the fractionalorder differentiator/integrator in terms of IIR filters. Mekhnache [26] presented a method for designing fullband and nonfullband IIR digital integrators with linear phase to approximate fractionalorder integrators with IIR filters. Tseng [27] presented the design problems of the digital Feller fractionalorder integrator. Benkhettou et al. [28] introduced a general notion of fractional derivative to define functions on arbitrary time scales and develop the basic tools for timescale fractional calculus.
However, these methods could not address the best rational approximation. In fact, if the obtained approximation function in the frequency domain is the best rational approximation function, then it can further improve the accuracy of approximation without increasing the orders. In this paper, a rational approximation method of fractional calculus operator ${s}^{\pm \alpha}$ ($0<\alpha <1$) based on the best rational approximation definition is presented. A design and optimization method of a fractionalorder PID controller is also implemented.
The approximation is a subset of rational functions on $\left[a,b\right]$, and it is the best approximation of a given irrational function $f\in [a,b]$ under the normal $\Vert f\Vert =\underset{a\le x\le b}{\mathrm{m}\mathrm{a}\mathrm{x}}\leftf\right(x\left)\right$. For fractional calculus operators $f={s}^{\pm \alpha}$, the best rational approximation function definition can be given.
Definition 1. Let $f={s}^{\pm \alpha}$, $0<\alpha <1$. For any point ${\omega}_{i}$ in the frequency range $[a,b]$, ${M}_{{R}_{mn}}^{\left(i\right)}$ is the amplitude of rational approximation function ${R}_{mn}\left(s\right)$ at point ${\omega}_{i}$ and ${M}_{f}^{\left(i\right)}$ is the amplitude of $f$ at point ${\omega}_{i}$. The best approximation function of $f$ in ${R}_{mn}[a,b]$ can be defined as:
where ${\omega}_{1}=a$, ${\omega}_{N}=b$ and ${R}_{mn}[a,b]$ is a rational function set:
This is the maximum amplitude error between ${R}_{mn}\left(s\right)$ and $f$. If the function $E\left({R}_{mn}^{\mathrm{*}}\right(s\left)\right)\in {R}_{mn}[a,b]$. It satisfies $E\left({R}_{mn}^{\mathrm{*}}\right(s\left)\right)=\mathrm{\Delta}(f;{R}_{mn}[a,b\left]\right)$.
Then ${R}_{mn}^{\mathrm{*}}\left(s\right)$ is called the best rational approximation function of $f$.
In Definition 1, the notation ${R}_{mn}\left(s\right)$ means a rational function to be expressed as follow:
${R}_{mn}[a,b]$ is a set, which is composed of ${R}_{mn}\left(s\right)$. ${R}_{mn}[a,b]$ is a rational function class.
The construction steps of the fractionalorder integral operator are shown in Fig. 1.
Fig. 1. Logarithmic amplitude frequency characteristic curve and the asymptote
To improve the approximation accuracy, the approximation range $[a,b]$ contains the interested approximation band [${\omega}_{a},{\omega}_{b}$].
The construction method is based on the minimum phase system. It can ensure that the phase frequency and amplitude frequency characteristics of the rational approximation function are consistent. The constructed rational approximation function is composed of n single real negative zeros and m single real negative poles. In order to achieve better approximation characteristics, let $nm=\left\{\mathrm{0,1}\right\}$. The fractionalorder integral operator is represented:
When ${k}_{0}=$1, ${G}_{0}\left(s\right)$ is a pure fractionalorder integral operator. The [${\omega}_{a},{\omega}_{b}$] is the interested approximation frequency band and ${\omega}_{c}$ is the crossing frequency of the logarithm amplitude frequency characteristic, ${\omega}_{a}<{\omega}_{c}<{\omega}_{b}$.
First, the curve of Eq. (3) is shown by using line ${L}_{0}$ in Fig. 1. The other asymptote line ${L}_{1}$ and line ${L}_{2}$ are parallel to line ${L}_{0}$. The amplitude error from ${L}_{1}$ and line ${L}_{2}$ to ${L}_{0}$ is $\epsilon $ (dB), which is the maximum permission error of the amplitude frequency characteristic. Line ${L}_{1}$ and line ${L}_{2}$ are given by:
Let $s=j\omega $; the slope of the logarithm amplitude frequency characteristic curve is –20$\alpha $ in Eq. (3) and (4). The following equations are obtained:
The relationship among ${L}_{0}\left(\omega \right)$, ${L}_{1}\left(\omega \right)$ and ${L}_{2}\left(\omega \right)$ is described as:
From Eq. (5) and (6), the following equation is obtained:
Next, the detailed steps for obtaining best rational approximation function of fractionalorder integral operator ${s}^{\alpha}$ are described:
Step 1. Select the first turnover frequency ${\omega}_{1}$.
Step 2. Draw a line segment with slope –20 dB/dec from the first turnover frequency ${\omega}_{1}$ on line ${L}_{1}$; this line segment can be represented as ${h}_{1}/s$, which intersects line ${L}_{2}$ at frequency ${\omega}_{1}^{\text{'}}$. For the frequency point ${\omega}_{1}$, the following expression is obtained:
From Eqs. (5) and (8), the following equation is obtained: ${h}_{1}={k}_{1}{{\omega}_{1}}^{1\mathrm{\alpha}}$.
For the frequency point ${\omega}_{1}^{\text{'}}$, the following equation is obtained:
From Eqs. (5) and (9), the following equation is obtained: ${\omega}_{1}^{\text{'}}=({h}_{1}/{k}_{2}{)}^{\frac{1}{1\alpha}}$.
Step 3. ${\omega}_{2}$ is obtained by crossing line ${L}_{2}$ with line ${L}_{1}$. ${\omega}_{2}$ meets the following equation:
From Eq. (10), the following equation is obtained: ${\omega}_{2}={\left({k}_{1}/{k}_{2}\right)}^{\frac{1}{\alpha}}{\omega}_{1}^{\text{'}}$.
The obtained ${\omega}_{1}^{\text{'}}$ and ${\omega}_{2}$ are the turnover frequencies of the firstorder lagging link and secondorder leading link respectively.
The previous steps are repeated until the turnover frequency ${\omega}_{1}^{\text{'}}$ of the $i$th order 1agging link and turnover frequency ${\omega}_{i+1}$ of the $(i+1)$th order leading link are obtained. The recursive calculative formulas are obtained:
where ${h}_{i}={k}_{1}{\omega}_{i}^{1\alpha}$, $i=\mathrm{}$1, 2,….
If ${\omega}_{1}^{\text{'}}$ and ${\omega}_{i+1}$ are far greater than ${\omega}_{b}$, they are removed and step 3 is terminated. At this time, an approximation function of ${G}_{0}\left(s\right)$ is obtained.
Step 4. In the approximation range [$a,b$], $N$ points with the same distance are selected, the difference between the value of ${G}_{0}^{i}\left(s\right)$ and value of ${R}_{mn}^{i}\left(s\right)$ ($i=$1, 2,…, $N$) is calculated, and the maximum difference is selected as $E\left({R}_{{\mathrm{}}_{mn}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\right(s\left)\right)$:
Step 5. Repeat step 14 for different values of ${\omega}_{1}$ to obtain $M$ rational approximation functions of ${G}_{0}\left(s\right)$.
Step 6. The minimum value of $E\left({R}_{mn}^{\mathrm{*}}\right(s\left)\right)$ is selected from ${E}_{k}\left({R}_{mn}^{\mathrm{m}\mathrm{a}\mathrm{x}}\left(s\right)\right)$ ($k=$1, 2,…, $M$).
Thus, $E\left({R}_{mn}^{\mathrm{*}}\right(s\left)\right)=\mathrm{\Delta}\left({G}_{0}\right(s);{R}_{mn}[a,b\left]\right)$, $s=j\omega $.
After all the steps, the approximation function of ${G}_{0}\left(s\right)$ is obtained:
where $K={k}_{1}/{\omega}_{1}^{\alpha}$.
Eq. (13) can be written as a general rational fractional expression:
where $nm=\left\{\mathrm{0,1}\right\}$.
To verify the effectiveness of the proposed method, ${s}^{0.3}$ is approximated separately by Oustaloup’s method and the proposed method, and the results of approximation are compared. The approximation frequency range $[a,b]$ is [1, 10000]. when the approaching order $N$ is 3, the approximation results are obtained as follows:
At the same time, the approximation frequency range $[a,b]$ is [1, 10000], and the maximum approximation error of the amplitude frequency characteristic curve is $\epsilon =$1.5 (dB); thus, the best rational approximation function of ${s}^{0.3}$ is written as:
The comparison of the Bode plots of the above two methods is shown in Fig. 2.
Fig. 2. Comparison result of the two methods
a) Amplitude frequency characteristic and error curve
b) Phase frequency characteristic and error curve
As can be seen from Fig. 2(a), on the basis of analyzing the amplitude frequency characteristic and error curve, the result indicates that the proposed method better approximates the ideal curve in the approximation range. At the same time, this method compensates for the deficiency of the approximation effect of Oustaloup’s method in the near boundary point of the approximation interval. As can be seen from Fig. 2(b), on the basis of analyzing the phase frequency characteristic and error curve, the result indicates that the proposed method approximates the ideal curve better than Oustaloup’s method in the frequency range [1, 10]. The two methods deliver better approximation effects in the frequency range [10, 450], have certain approximation errors in the frequency range [450, 10000] and require further study.
In summary, on the basis of analyzing the amplitude frequency characteristic, phase frequency characteristic and error curve, the proposed method has a better approximation effect and compensates for the deficiency of the approximation effect of Oustaloup’s method to a certain extent.
Because this paper involves only the design idea for one fractionalorder controller, it is independent of the values of the parameters. Thus, the selection parameters are not discussed, and the determined parameter values are randomly selected.
First, the fractionalorder $P{I}^{\alpha}{D}^{\beta}$ controller in the frequency domain is given:
The design steps of the fractionalorder $P{I}^{\alpha}{D}^{\beta}$ controller can be described.
Step 1. The values of the parameters (${K}_{p},{K}_{i},{K}_{d},\alpha ,\beta $) are initialized.
Step 2. The best rational approximation functions of fractionalorder integral operators ${s}^{\alpha}$ and ${s}^{\beta}$ are calculated by using the proposed method.
Step 3. The best rational approximation functions of ${s}^{\alpha}$ and ${s}^{\beta}$ are substituted into Eq. (17), and a new equation ${G}_{c1}\left(s\right)$ is obtained.
As can be seen from Eq. (17), the fractionalorder controller includes five parameters: ${K}_{p}$, ${K}_{i}$, ${K}_{d}$, $\alpha $, $\beta $, $0<\alpha ,\beta <1$.
In this paper, let ${K}_{p}=1$, ${K}_{i}=1$, ${K}_{d}=1$, $\alpha =0.3$, $\beta =0.6.$
Five parameters of controller are substituted into Eq. (17), and Eq. (18) is obtained:
Because the design is based on the best rational approximation method of fractionalorder integral operators and the corresponding fractionalorder integral operator is equivalent to a rational function, Eq. (18) can be rewritten as:
where ${s}^{\alpha}=A\left(s\right)/B\left(s\right)$, ${s}^{\beta}=C\left(s\right)/D\left(s\right)$.
Let the approximation frequency band [${\omega}_{a}$, ${\omega}_{b}$] is [1,1000] and the maximum approximation error of amplitude frequency characteristic $\epsilon $ is 4 (dB). The best rational approximation functions of fractionalorder integral operators ${s}^{0.3}$ and ${s}^{0.6}$ are calculated by the proposed method. Thus the best rational approximation equations of ${s}^{0.3}$ and ${s}^{0.6}$ are obtained:
Eqs. (20) and (21) are substituted into Eq. (19). Thus Eq. (22) is obtained:
Eq. (22) is the best rational approximation transfer function of fractionalorder controller. Eq. (22) is transformed into Eq. (23) by the coefficient transformation:
Eq. (23) is used to describe the integerorder transfer function for the corresponding fractionalorder controller. The amplification factor of the integerorder transfer function is ${K}_{k}=$541.8. Its zeros and poles are shown in Table 1.
Table 1. The zeros and the poles of integer order transfer function
No.

Zeros

Poles

1

–1.113

–1.111

2

–88.9

–89.24

3

–7162

–7166

4

–1701

–7.691

5

–13.9+j163.6

–357

6

–13.9j163.6

To optimize and analyze the logarithmic frequency characteristics of transfer function, $j\omega $ is used to replace $s$ in Eq. (23), and Eq. (23) is transformed into Eq. (24):
Next, the frequency characteristic of fractionalorder controller is optimized and analyzed by using Eq. (24). The optimization removes partial zeros and poles, and the amplification factor ${K}_{k}$ of the transfer function is a constant. The steps of the optimization are explained:
1) The dipoles $\left(j\frac{\omega}{7162}+1\right)$ and $1/\left(j\frac{\omega}{7166}+1\right)$ outside the frequency band are removed, and the following Eq. (25) is obtained:
2) The zeros outside the approximation frequency [1, 1000] are removed, and the following transfer function of the controller is obtained:
3) Based on Eq. (26), the dipoles $\left(j\frac{\omega}{1.113}+1\right)$ and $1/\left(j\frac{\omega}{1.111}+1\right)$ are removed and the transfer function is obtained:
4) Based on Eq. (27), the dipoles $\left(j\frac{\omega}{88.9}+1\right)$ and $1/\left(j\frac{\omega}{89.24}+1\right)$ are removed and the transfer function is obtained:
5) Based on Eq. (28), one order leading link $\left(j\frac{\omega}{1701}+1\right)$ outside the approximation frequency [1, 1000] is removed, and the transfer function is obtained:
The corresponded logarithmic frequency characteristic curves of Eqs. (24)(29) are shown in Fig. 3Fig. 8, and the result of the comparison is shown in Fig. 9.
The Bode diagram of Eq. (24) is shown in Fig. 3. As can be seen from Fig. 3, the amplitude and phase of the controller changed greatly in the frequency band [100, 300] (rad/sec), and the maximum phaselag was close to 90°. Fig. 4 shows the frequency characteristic curve obtained by removing a dipole nearby at 7000 rad/sec in Fig. 3. Because the dipole is far from the approximation band, Fig. 3 and Fig. 4 within the given frequency band do not change. Fig. 5 is obtained by removing the one order leading link outside the approximation frequency in Fig. 4. As can be seen from Fig. 5, the amplitude and phase became smaller in the highfrequency band. Because the one order leading link is much closer to the approximation frequency band, it exerts some influence on the highfrequency, as compared with Fig. 3.
According to the principle of automatic control, the reduction of amplitude will be beneficial for suppressing highfrequency interference, and the phase is still positive. Thus, the designed controller cannot decrease the system phase angle, which affects the system stability.
Fig. 3. Bode diagram of the original rational approximation transfer function
Fig. 4. Bode diagram of the transfer function corresponding to Eq. (25)
Two dipoles are located at the frequencies 1 rad/sec and 89 rad/sec on the approximation band. Fig. 6 shows a Bode diagram obtained by removing one dipole in Fig. 4. Fig. 7 shows a Bode diagram obtained by removing two dipoles in Fig. 4. Because the distance between the pole and zero is far smaller than the distance from the origin to the pole or zero, the influence of the amplitude and phase of the transfer function can be completely offset. Thus, Fig. 6 and Fig. 7 are close to Fig. 4. As can be seen from Figs. 47, the orders of the controller decreased; thus, the complexity is reduced by removing the dipole in the obtained transfer function.
Fig. 5. Bode diagram of the transfer function corresponding to Eq. (26)
Fig. 6. Bode diagram of the transfer function corresponding to Eq. (27)
To further simplify the structure of the controller based on maintaining the performance of the controller, Eq. (29) will be regarded as the controller model, and the order of the numerator and denominator is set to 2. The corresponding Bode diagram is shown in Fig. 8.
Fig. 7. Bode diagram of the transfer function corresponding to Eq. (28)
Fig. 8. Bode diagram of the transfer function corresponding to Eq. (29)
Fig. 9. Comparison of bode plots
Fig. 9. is the comparison of the Bode diagrams. All the discussed characteristic curves are described in Fig. 9. It clearly shows that the removed dipoles will rarely affect the characteristic curves during the order reduction, but the reduction of one order leading link $\left(j\frac{\omega}{1701}+1\right)$ will obviously affect the high frequency of the characteristic curve. The final rational approximation transfer function of the fractionalorder controller is shown below:
Based on the proposed best rational approximation definition of the fractional calculus operators, the construction method for the best rational approximation of a fractionalorder integral operator is presented in this paper. The obtained approximation function in the frequency domain is the best rational approximation function, which can further improve the accuracy of approximation without increasing orders. The design and optimizing process of a fractionalorder PID controller are also discussed. The final obtained fractionalorder controller presents a good frequency characteristic, the best structure, easy realization and so on. This method is convenient, targeted to determining the controller model according to the given frequency band, and can be easily implemented with the program. This method can also provide better engineering applications of the fractionalorder controller.