# Rational Implementation of Fractional Calculus Operator Based on Quadratic Programming

Authors: Zhisong Xu; Mingqiu Li
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

When fractional calculus operators and models are implemented rationally, there exist some problems such as low approximation accuracy of rational approximation function, inability to specify arbitrary approximation frequency band, or poor robustness. Based on the error criterion of the least squares method, a quadratic programming method based on the frequency-domain error is proposed. In this method, the frequency-domain error minimization between the fractional operators±r and its rational approximation function is transformed into a quadratic programming problem. The results show that the construction method of the optimal rational approximation function of fractional calculus operator is effective, and the optimal rational approximation function constructed can effectively approximate the fractional calculus operator and model for the specified approximation frequency band.


---

## Body

## 1. Introduction

Fractional calculus is a mathematical problem to study integral and differential of arbitrary order. Fractional calculus has many applications, such as physical system modeling [1–3], control theory [4–6], and so on. In control theory, as fractional calculus is more and more applied to the design and analysis of controllers and filters, the research on approximation methods of fractional order systems or operator models has attracted more and more attention [7–9]. However, compared with the integer order system, fractional order system is often more complex, and its essence is infinite dimension system. Its characteristic equation is a polynomial with complex variable fractional order index, so it is difficult to find the analytic solution. Therefore, the fractional order system or operator model is usually approximated by the way of direct approximation in z-domain or by transforming the fractional system into rational transfer function.The direct approximation method [10] is to use a function of finite order to approximate fractional operators in z-domain, and different generating functions are usually applied to different discrete operators, for example, power series expansion (PSE) of Euler operator and continuous fractional expansion (CFE) of Tustin operator. It should be noted that when the CFE method is used to discretize the closed-loop continuous fractional transfer function, the model may be unstable. Some operators, such as Al-Alaoui operator [11], are obtained by trapezoidal interpolation and rectangular integration rules. In some cases, the interpolation and inversion process may lead to unstable fractional order models. The indirect discretization of fractional calculus operator includes two steps: approximation of transfer function in s domain and discretization of approximate function. French scholar Professor Oustaloup and his colleagues proposed filter of fractional order operator [12]. This kind of filter allows users to choose the frequency band and order of interest and use the integer order transfer function model to approximate the fractional order calculus operator. For those irrational systems which cannot be described by the standard form of fractional transfer functions, the fractional calculus operators can be discretized indirectly by frequency response fitting or Charef filter [13]. The fractional order system can also be approximated by the differential evolution method, and many scholars participated in the research. Among them, the optimization performance of differential evolution algorithm depends on the control parameters and mutation strategy to a large extent [14].In this paper, for the first time, we use the minimization of the frequency error between the calculus operator and the rational approximation function to approximate the fractional calculus operator. Firstly, the frequency range to be approximated is selected, and the approximating frequency band is equally divided inton equal intervals. The equal points are called characteristic points. On one characteristic point, the difference between the frequency characteristic function of the ideal calculus operator and the approximating function is solved. Then, the sum of the frequency differences at n+1 points is added. The expression of the sum only contains quadratic and first-order terms. The expression of the sum is optimized by quadratic programming. Then, the coefficients of the rational expression of the approximation function are obtained, and the optimal approximation expression is obtained. This quadratic programming problem can be solved by quadprog command in MATLAB. Compared with the previous methods, the method proposed in this paper has obvious improvement in approximation accuracy, precision, and algorithm encapsulation and can arbitrarily select the approximation frequency band, which makes the algorithm highly flexible and has great advantages.


## 2. Basic Knowledge of Fractional Calculus

### 2.1. Definition of Fractional Calculus

The general expression of fractional calculus is [15](1)Dtαft=dαdtαft,Reα>0,ft,Reα=0,∫0tfτdτ−α,Reα<0,where Dtα is a fractional differential or integral operator, 0 and T are the upper and lower limits (the default initial condition is 0) of calculus, and α can be a real number or a complex number.Different mathematicians give different definitions of fractional calculus. The three definitions of fractional calculus which are widely used in control field include Grünwald–Letnikov (GL), Riemann–Liouville (RL), and Caputo.


#### 2.1.1. Grünwald–Letnikov (GL) Definition

Suppose that the functionft has n+1 derivative in the interval 0,t. For any real number (the lower limit of integration at the initial time of 0 is 0), the α-order calculus of function ft is defined as(2)Dtαft=limh⟶0h−α∑j=0t/h−1jαjft−jh,where ⋅ represents the number of approximate recurrence terms, in which(3)αj=αα+1…α+j−1j!=α!j!α−j!,is the coefficient of recursive function.


#### 2.1.2. Riemann–Liouville (RL) Definition

RL definition is based on the definition of GL by simplifying the calculation process. For any real numbern−1<α<n, n∈N, the RL fractional order is defined as(4)Dtαft=1Γn−αddtn∫0tfτt−τα−n+1dτ,where Γ⋅ is a gamma function.RL fractional integral is defined as(5)Itαft=1Γ−α∫0t1t−τα−n+1dτ,t>0,α∈R+.


#### 2.1.3. Caputo Definition of Fractional Calculus

Caputo fractional calculus is defined as(6)Dtα=1Γm−α∫0tfmτt−τα−m+1dτ.Among them,m−1<α<m. It can be seen from equation (5) that Caputo definition requires the first m-order derivative of the function to be integrable.


## 2.1. Definition of Fractional Calculus

The general expression of fractional calculus is [15](1)Dtαft=dαdtαft,Reα>0,ft,Reα=0,∫0tfτdτ−α,Reα<0,where Dtα is a fractional differential or integral operator, 0 and T are the upper and lower limits (the default initial condition is 0) of calculus, and α can be a real number or a complex number.Different mathematicians give different definitions of fractional calculus. The three definitions of fractional calculus which are widely used in control field include Grünwald–Letnikov (GL), Riemann–Liouville (RL), and Caputo.


### 2.1.1. Grünwald–Letnikov (GL) Definition

Suppose that the functionft has n+1 derivative in the interval 0,t. For any real number (the lower limit of integration at the initial time of 0 is 0), the α-order calculus of function ft is defined as(2)Dtαft=limh⟶0h−α∑j=0t/h−1jαjft−jh,where ⋅ represents the number of approximate recurrence terms, in which(3)αj=αα+1…α+j−1j!=α!j!α−j!,is the coefficient of recursive function.


### 2.1.2. Riemann–Liouville (RL) Definition

RL definition is based on the definition of GL by simplifying the calculation process. For any real numbern−1<α<n, n∈N, the RL fractional order is defined as(4)Dtαft=1Γn−αddtn∫0tfτt−τα−n+1dτ,where Γ⋅ is a gamma function.RL fractional integral is defined as(5)Itαft=1Γ−α∫0t1t−τα−n+1dτ,t>0,α∈R+.


### 2.1.3. Caputo Definition of Fractional Calculus

Caputo fractional calculus is defined as(6)Dtα=1Γm−α∫0tfmτt−τα−m+1dτ.Among them,m−1<α<m. It can be seen from equation (5) that Caputo definition requires the first m-order derivative of the function to be integrable.


## 2.1.1. Grünwald–Letnikov (GL) Definition

Suppose that the functionft has n+1 derivative in the interval 0,t. For any real number (the lower limit of integration at the initial time of 0 is 0), the α-order calculus of function ft is defined as(2)Dtαft=limh⟶0h−α∑j=0t/h−1jαjft−jh,where ⋅ represents the number of approximate recurrence terms, in which(3)αj=αα+1…α+j−1j!=α!j!α−j!,is the coefficient of recursive function.


## 2.1.2. Riemann–Liouville (RL) Definition

RL definition is based on the definition of GL by simplifying the calculation process. For any real numbern−1<α<n, n∈N, the RL fractional order is defined as(4)Dtαft=1Γn−αddtn∫0tfτt−τα−n+1dτ,where Γ⋅ is a gamma function.RL fractional integral is defined as(5)Itαft=1Γ−α∫0t1t−τα−n+1dτ,t>0,α∈R+.


## 2.1.3. Caputo Definition of Fractional Calculus

Caputo fractional calculus is defined as(6)Dtα=1Γm−α∫0tfmτt−τα−m+1dτ.Among them,m−1<α<m. It can be seen from equation (5) that Caputo definition requires the first m-order derivative of the function to be integrable.


## 3. Construction of Frequency-Domain Error Function

The principle diagram of approximation of fractional order calculus operator by using the frequency-domain error minimization principle [16] is shown in Figure 1. Among them, r is the order of the ideal fractional calculus operator, G^s represents the best approximate rational expression of jωr, and ejω is the error between rational approximation function G^jω and ideal fractional calculus operator jωr in frequency domain.Φ⋅ means to use some optimization algorithm for error ejω, such as genetic algorithm, linear programming, model reduction method, and so on. Φet represents the optimized frequency-domain error function. The parameter of G^s is obtained by calculating the partial derivative of each parameter of Φet. Different from other methods, in this paper, Gs (fractional transfer function model to be optimized) is changed into the frequency-domain model of fractional order calculus operator jωr(sr|s=jω), reducing the time and complexity of the program.Figure 1
Schematic diagram of approximation of fractional calculus operator by frequency-domain error minimization principle.For the fractional calculus operatorsr to be approximated, it can be expanded into the following form:(7)sr=jωr=ωrejπ2r=ωrcosπ2r+jsinπ2r=ωrcosπ2r+jsinπ2r.After optimization by some optimization algorithm, the rational approximation functionG^jω is obtained. Its form is shown in the following formula:(8)G^jω=dmjωm+dm−1jωm−1+⋯+d1jω+d0cnjωn+cn−1jωn−1+⋯+c1jω+c0=NωDω(9)=Reω+jImωRe′ω+jIm′ω,where ci,dii=1,2,3… are the coefficient parameters of approximation expression G^jω. Re and Im represent the real and imaginary parts of the expressions. Within the specified approximation frequency band, the equidistant frequency sampling points are taken as ωii=1,2,3…. Suppose that l+1 characteristic frequency points are selected in the approximation frequency range. Among them,(10)Δωi=ωh−ωbl1≤i≤l+1,and at these points, the transfer function is approximated:(11)G^jωi=Reωi+jImωiRe′ωi+jIm′ωi.For each frequency pointωi:(12)sir=ωircosπ2r+jωirsinπ2r≈Reωi+jImωiRe'ωi+jIm'ωi,and at each characteristic frequency point, the frequency error is εi. Then,(13)εi=ωircosπ2r+jωirsinπ2r⋅Re′ωi+jIm′ωi−Reωi−jImωi=ωircosπ2rRe′ωi−ωirsinπ2rIm′ωi−Reωi+jωircosπ2rIm′ωi+ωirsinπ2rRe′ωi−Imωi,and b=1≤i≤l+1=h.


## 4. Transformation of Quadratic Programming

To solve the problem of minimizing the frequency-domain error, Levy proposed a method [16]. In this method, an original transfer function model of a fractional order linear time invariant continuous system is considered, as shown in the following equation:(14)Gs=b0sβ0+b1sβ1+b2sβ2+⋯+bmsβma0sα0+a1sα1+a2sα2+⋯+amsαn,and the frequency-domain error is(15)εω=Gjω−NωDω.Let(16)Eω=εωDω=GjωDjω−Nω=ReG+jImGRe′ω+jIm′ω−Reω+jImω=ReGRe′ω−Reω+jReGIm′ω+ImGRe′ω−Imω,(17)Eω2=ReGσ−ImGτ−α2+ReGτ+ImGσ−β2.In order to get the minimum value of the frequency-domain error,Eω2 needs to find the partial derivatives of parameters ci and di, respectively. Let the partial derivative be zero and solve the system of equations; we can get the minimum value of Eω2.For the above traditional optimization methods, it is necessary to calculate the partial derivative of each parameter and solve large-scale equations, which is very cumbersome and takes a long time to solve. Moreover, the programming is cumbersome, which is not conducive to the algorithm encapsulation. In order to overcome these shortcomings, an error-based quadratic programming method is proposed.For equation (13), take another treatment:(18)εi2=ωircosπ2rRe′ωi−ωirsinπ2rIm′ωi−Reωi2+ωircosπ2rIm′ωi+ωirsinπ2rRe′ωi−Imωi2,and then(19)εi2=ωi2rRe′2ωi+Im′2ωi+Re2ωi+Im2ωi+2ωirsinπ2rReωiIm′ωi−ImωiRe′ωi−2ωircosπ2rReωiRe′ωi+ImωiIm′ωi.LetX=c0c1…c2Kd0d1…d2NT, and X≠0⟶, where(20)Re′ωi=∑m=0K−1mc2mωi2m,(21)Im′ωi=∑m=1K−1m−1c2m−1ωi2m−1,(22)Reωi=∑m=0N−1md2mωi2m,(23)Imωi=∑m=0N−1m−1d2m−1ωi2m−1,(24)K,N=0,1,2…….Equation (19) is the quadratic form of the frequency-domain error obtained. The sum of the squares of the errors of each frequency sampling point in the approximate frequency band is ε; then,(25)ε=∑i=1l+1εi2,and (19)–(23) and (25) can be arranged into the following matrix form:(26)ε=XTQX,Q=0010…0⋮⋮0ω1…ωl+100−ω120…0100−ω13…−ωl+130ω1ω140…0−ω120⋮⋮…⋮0−ω13−1Kω12K0…0ω14000…0⋮⋮⋮⋮…⋮−1Nω12N000…0∗A1⋯0⋮⋱⋮0⋯Al+1∗0010…0⋮⋮0ω1…ωl+100−ω120…0100−ω13…−ωl+130ω1ω140…0−ω120⋮⋮…⋮0−ω13−1Kω12K0…0ω14000…0⋮⋮⋮⋮…⋮−1Nω12N000…0T=CACT,Ai=10−ωircosπ2rωirsinπ2r01−ωirsinπ2r−ωircosπ2r−ωircosπ2r−ωirsinπ2rωi2r0ωirsinπ2r−ωircosπ2r0ωi2r1≤i≤l+1.Obviously, the Q-matrix is only related to the approximation frequency segment, the order of the approximated fractional calculus operator, the number of sampling points, the denominator of the approximation function, and the order of the numerator.For fractional calculus operatorf=s±r,0<r<1, there is an optimal rational approximation function in s domain because under the premise of the same calculation amount and calculation time, the approximation expression obtained by a rational function with the same order of numerator and denominator is much better than that obtained by polynomial approximation method. When constructing the rational approximation expression (i.e., transfer function), m single negative real zeros and n single negative real poles are used, and n−m=0,1; thus, better approximation properties of fractional calculus operators can be obtained [17]. In this paper, the simulation example part adopts the case that the denominator of the approximation transfer function is the same as the numerator order and is 6.Thus, the approximation problem of fractional calculus operator is transformed into the quadratic programming problem of frequency error. For linear quadratic programming problems [18]:(27)minx12xTQx+fTx⟶yieldsA⋅x≤b,Aeq⋅x=beq,lb≤x≤ub,where Q is a symmetric matrix, A⋅x≤b is a linear constraint condition, Aeq⋅x=beq is equality constraint, and lb≤x≤ub is independent variable definition domain.(28)x=quadprogH,f,A,b,Aeq,beq,lb,ub,x0,options.The command in formula (28) can be used to solve the problem.


## 5. Algorithm Sensitivity Analysis

Sensitivity analysis refers to the sensitivity analysis of the system or surrounding things due to changes in the surrounding conditions. If the market conditions change, the value will change. This paper mainly discusses the influence of the parameter vectorX=c0c1…c2Kd0d1…d2NT on the approximation effect.Parameters such as the approximation frequency segment, the order of the approximated fractional calculus operator, the number of sampling points, the denominator of the approximation function, and the order of the numerator are set in advance, so they are not sensitive to the approximation algorithm. The only sensitive quantity isX. Next, we analyze the sensitivity caused by the change of parameters in the X vector.Firstly, the value of nonsensitive quantity is determined: assuming that the frequency band of approximation is1,100 rad/s, the value of γ is 0.5, the frequency sample point is taken as 40, and the denominator of the approximation transfer function is the same as the numerator order and is 6. From the previous analysis, the fitting effect of this method in the initial part of the approximation frequency band is not ideal, so the sensitivity analysis of the initial frequency band is carried out.Approximation function is(29)G1=0.05049s6+20.03s5+1663s4+33590s3+100000s2+0.001s+0.0010.001s6+1.269s5+217.1s4+9023s3+79840s2+36930s+100000.At the initial approximation frequency band, the rational approximation function can be reduced to(30)G1′s=33590s3+100000s279840s2+36930s+100000.The main parameters that affect the system (the quantity that determines the sensitivity) arec0,c1,c2,d2,d3.Whenc0,c1,c2,d2,d3 increases or decreases by 10% alone, the difference between the approximation phase angle generated by the rational approximation function of the system and the ideal phase angle is shown in Tables 1 and 2.Table 1
The change of each index after the parameters increase by 10% (α′ is the approximate phase angle of the parameters’ changed value at 2 rad/s and α−α′ is the difference between the phase angle of ideal fractional differential operator and the approximate phase angle after the parameters change).
The original costIncrease by 10 percentα′ (°)α−α′ (°)c010000011000036.38.7c1369304062338.36.7c2798408782436.18.9d210000011000033.511.5d3335903694938.86.2Table 2
The change of each index after the parameters decrease by 10% (α′ is the approximate phase angle of the parameters’ changed value at 2 rad/s and α−α′ is the difference between the phase angle of ideal fractional differential operator and the approximate phase angle after the parameters change).
The original costDecrease by 10 percentα′ (°)α−α′ (°)c01000009000036.28.8c1369303323734.310.7c2798407185636.18.9d21000009000039.35.7d3335903023133.411.6Investigate the sensitivity index when each parameter changes; the following can be obtained:(31)δ=α−α′±10%⋅α′.Therefore, we can derive the sensitivity when a single weight parameter changes, as shown in Table3.Table 3
Sensitivity of algorithm with single parameter change.
δ1(increase by 10%)δ2(decrease by 10%)c02.3972.431c11.7493.120c22.3242.465d23.4331.450d31.6193.473It can be seen from Table3 that when the parameters that affect the sensitivity of the algorithm increase by 10%, the sensitivity of parameter d3 is the lowest; the sensitivity of parameter d2 is the highest; when the parameters that affect the sensitivity of the algorithm decrease by 10%, the sensitivity of parameter d2 is the lowest and the sensitivity of parameter d3 is the highest.Table 4
Phase angle of characteristic frequency in approximation frequency band of different rational implementation methods.
10 (rad/s)40 (rad/s)70 (rad/s)100 (rad/s)x¯ (sample mean)S (standard deviation)Ideal45.0°45.0°45.0°45.0°45.0°0°Oustaloup19.7°45.0°45.0°19.6°32.3250°14.6359°This paper45.6°41.23°45.1°45.0°44.2325°2.0188°


## 6. Algorithm Complexity Analysis

### 6.1. Time Complexity

The time complexity of the proposed algorithm is dominated by two loop statements, which arefor i = 1:nA(:,:,i) = [0 0 1 0; 0 0 0 Omega(i); 0 0 -Omega(i)^2 0; 0 0 0 -Omega(i)^3; 0 0 Omega(i)^4 0; 0 0 0 Omega(i)^5; 0 0 -Omega(i)^6 0; 1 0 0 0; 0 Omega(i) 0 0; -Omega(i)^2 0 0 0; 0 -Omega(i)^3 0 0; Omega(i)^4 0 0 0; 0 Omega(i)^5 0 0; -Omega(i)^6 0 0 0];endandfor i = 1:nC(:,:,i) = [1 0 -Omega(i)^r∗cos(pi∗r/2) Omega(i)^r∗sin(pi∗r/2); 0 1 -Omega(i)^r∗sin(pi∗r/2) -Omega(i)^r∗cos(pi∗r/2); -Omega(i)^r∗cos(pi∗r/2) -Omega(i)^r∗sin(pi∗r/2) Omega(i)^(2∗r) 0; Omega(i)^r∗sin(pi∗r/2) -Omega(i)^r∗cos(pi∗r/2) 0 Omega(i)^(2∗r)];endThe number of execution steps of the two loop statements isn and n, respectively, and the remaining assignment and execution statements of the algorithm are 15 sentences. Therefore, the total number of execution steps of the algorithm is 2n + 15. Therefore, the complexity of the proposed algorithm is O(n).


### 6.2. Spatial Complexity

This algorithm does not use recursive algorithm, and the number of temporary variables in the algorithm has nothing to do with the problem sizen, so the space complexity is O(1).


## 6.1. Time Complexity

The time complexity of the proposed algorithm is dominated by two loop statements, which arefor i = 1:nA(:,:,i) = [0 0 1 0; 0 0 0 Omega(i); 0 0 -Omega(i)^2 0; 0 0 0 -Omega(i)^3; 0 0 Omega(i)^4 0; 0 0 0 Omega(i)^5; 0 0 -Omega(i)^6 0; 1 0 0 0; 0 Omega(i) 0 0; -Omega(i)^2 0 0 0; 0 -Omega(i)^3 0 0; Omega(i)^4 0 0 0; 0 Omega(i)^5 0 0; -Omega(i)^6 0 0 0];endandfor i = 1:nC(:,:,i) = [1 0 -Omega(i)^r∗cos(pi∗r/2) Omega(i)^r∗sin(pi∗r/2); 0 1 -Omega(i)^r∗sin(pi∗r/2) -Omega(i)^r∗cos(pi∗r/2); -Omega(i)^r∗cos(pi∗r/2) -Omega(i)^r∗sin(pi∗r/2) Omega(i)^(2∗r) 0; Omega(i)^r∗sin(pi∗r/2) -Omega(i)^r∗cos(pi∗r/2) 0 Omega(i)^(2∗r)];endThe number of execution steps of the two loop statements isn and n, respectively, and the remaining assignment and execution statements of the algorithm are 15 sentences. Therefore, the total number of execution steps of the algorithm is 2n + 15. Therefore, the complexity of the proposed algorithm is O(n).


## 6.2. Spatial Complexity

This algorithm does not use recursive algorithm, and the number of temporary variables in the algorithm has nothing to do with the problem sizen, so the space complexity is O(1).


## 7. Simulation Experiment and Result Analysis

### 7.1. Pure Fractional Differential Operators

The expression of fractional differential operator isGs=sγ0<γ<1. In engineering, the Oustaloup filter is often used to approximate the fractional calculus operator in a given frequency band, assuming that the frequency range of interest is ωb,ωh; then, a set of broken lines can be used to approximate the linear characteristics of fractional calculus. Based on this idea, French scholar Professor Oustaloup [12] put forward the design method of CRONE filter. All these broken lines are generated alternately by zeros and poles of integer order so that the slope of the amplitude frequency characteristic asymptote changes alternately between 0 dB/dec and −20 dB/dec, and the frequency-domain response itself will be very close to an oblique line.The standard form of Oustaloup filter is shown in formula (29):(32)Gs=K∏k=1Ns+ωk′s+ωk,where k = 1,2, …, N, and zeros ωk′, poles ωk, and gain K can be calculated as follows:(33)ωk′=ωbωu2k−1−γ/N,(34)ωk=ωbωu2k−1+γ/N,(35)K=ωhγ,(36)ωu=ωhωb,where ωu is called unit gain frequency or transition frequency. Assuming that the frequency band of approximation is 10,100 rad/s, the value of γ is 0.5, and the number of iterations is 5; then, using the Oustaloup filter method [12], the approximation function can be obtained as follows:(37)GOustaloups=10s5+1726s4+107800s3+3039000s2+38650000s+177800000s5+217.4s4+17090s3+606300s2+9709000s+56230000.Using the quadratic programming method based on frequency-domain error, the frequency sample point is taken as 40, and the new approximation function is obtained as follows:(38)Gnews=0.05053s6+20.08s5+1668s4+33730s3+100000s2+0.001s+0.0010.001s6+1.271s5+217.7s4+9057s3+80130s2+35630s+100000.The amplitude phase characteristic curve is obtained by MATLAB simulation.It can be seen from Figure2 that the quadratic programming method based on the frequency-domain error can approach the pure fractional differential operator well in the specified frequency band. For the approximation of phase angle, it can be seen from Table 1 that the approximation effect of the quadratic programming method based on the frequency-domain error is better than that of the Oustaloup method, and it is closer to the theoretical value.Figure 2
Ideal Bode graph ofs0.5 and Bode graph obtained by two approximation methods.According to the sample standard deviation, we can measure the dispersion of the sample data. The formula is shown in (36):(39)s=1N−1∑i=1Nxi−x¯2,where N represents the sample population, xi represents the sample, x¯ represents the sample mean. According to the results of sample standard deviation in Table 4, the proposed method has smaller standard deviation, more concentrated data, and better approximation effect. The mean value is also closer to the ideal value.Table 5
Phase angle of characteristic frequency in approximation frequency band of different rational implementation methods.
2 rad/s34 rad/s67 rad/s100 rad/sx¯ (sample mean)S (standard deviation)Ideal−45.0°−45.0°−45.0°−45.0°−45.0°0°Oustaloup−32.4°−35.0°−27.8°−22.2°−29.35°5.6199°This paper−45.1°−45.0°−45.0°−45.0°−45.025°0.05°With the increase of the number of sample points, the quadratic programming method based on frequency error can obtain the amplitude frequency characteristic curve which basically coincides with the theoretical curve. In general, the quadratic programming method based on frequency error is better to approximate pure fractional differential operators.The frequency-domain characteristics of fractional differential operators with different orders are compared as shown in Figure3 (assuming that the approximation frequency segment is 1,100 rad/s, the sampling frequency point is 40, the order of approximation expression is 6, and the order is 0.1, 0.3, 0.5, 0.7, and 0.9, respectively).Figure 3
Approximation effect of fractional differential operators with different characteristic orders.Based on the frequency-domain quadratic programming method, the fitting effect of fractional order operator is very good in the middle and high frequency band of approaching frequency band, and the fractional order differential operator of each characteristic order has no large amplitude fluctuation in the frequency band. When the order is smaller (close to 0), the effect of approximation is better, and the amplitude and phase characteristics almost have no oscillation.


### 7.2. Pure Fractional Integral Operator

The expression of fractional integral operator isG=s−γ (0<γ<1). Assuming that the approximation frequency segment is 1,100 rad/s and the sampling frequency is 40, the order of approximation expression is 6 and the order is −0.5. The amplitude phase characteristic curve is obtained (as shown in Figure 4) by using the Oustaloup method and quadratic programming method based on the frequency-domain error proposed in this paper.Figure 4
Comparison of the effect of different methods to approximate fractional integral operator.Both the Oustaloup method and the method proposed in this paper can approach the fractional integral operator well. The Oustaloup method cannot fit the amplitude characteristics of fractional integral operators well in the high frequency band. For the fitting of phase frequency characteristics, refer to the data comparison in Table5 (the ideal amplitude and phase characteristic curve of integral operator s−0.5 is selected as the reference value).From the data in Table5, we can see that the phase frequency characteristics obtained by using this method to approximate the fractional order integral operator are consistent with the ideal phase frequency characteristics in the range from 34 rad/s to100 rad/s. The mean and standard deviation of the proposed method are smaller, the data are more stable, and the mean value is closer to the ideal value, so the approximation method proposed in this paper can more accurately approximate the ideal value, and the approximation data are more stable. It can be seen from Figure 4 that the phase frequency characteristic curves obtained by this method coincide with the ideal phase frequency characteristic curves in the middle and high frequency range. A better approximation effect is obtained.The frequency-domain characteristics of different order fractional integral operators are compared as shown in Figure5 (assuming that the approximation frequency band is 1,100 rad/s, the sampling frequency point is 40, and the order of approximation expression is 6, the integral characteristic order is −0.1, −0.3, −0.5, −0.7, −0.9).Figure 5
Approximation effect of fractional integral operators with different characteristic orders.It can be seen from Figure5 that the new method proposed in this paper can well approximate the characteristic fractional integral operator without oscillation in the approximation frequency range. The smaller the absolute value of integral order is, the better the approximation effect is.


### 7.3. Fractional Order System Model

Some common fractional order system models are simulated. For example, the fractional order model of PMSM is given as [19](40)Gs=9098.5417s1.7781+397316.2297s0.8371+1393070.1646s3.7915+127.38s2.8834+9995.678s1.8834.Original system cutoff frequencyωc = 40.8 rad/s. The constant parameters in this simulation are the approximation frequency (1,100 rad/s), the sampling frequency point (40), and the order of approximation expression (6). The results are shown in Figure 6.Figure 6
Approximation of ideal fractional order model by different methods.It can be seen from the amplitude frequency characteristic curve in Figure6 that the amplitude frequency characteristics of the system obtained by the two approximation methods are close. For the phase frequency characteristics obtained by the two approximation methods, in the second half of the approximation frequency band, the quadratic programming method based on the frequency-domain error can better fit the ideal phase frequency characteristics of the system than the Oustaloup method.


### 7.4. Comparison of New Method and Direct Discretization Method

The Muir recursive method is a direct discretization method commonly used in engineering. This method was first applied to geophysical scientific data processing and oil exploitation [10]. Without losing generality, we assume that the order of calculus operator is γ∈−1,1. For simplicity, we only give the recursive formula when γ is positive:(41)ωz−1r=2Tr1−z−11+z−1r=2Trlimn⟶∞Anz−1,rAnz−1,−r,where(42)A0z−1,r=1,(43)Anz−1,r=An−1z−1,r−cnznAn−1z,r,(44)cn=rn,nis odd,0,nis even.MATLAB toolbox was used, and the result is obtained as(45)sr≈2TrAnz−1,rAnz−1,−r.When the sampling frequencyT is 0.001, the number of iterations is 3. The fractional differential operator s0.5 can be directly discretized by the Muir recursive method in the following form [10]:(46)Gz=44.72z3−22.36z2+3.727z−7.454z3+0.5z2+0.8333z+0.1667.Compared with the method proposed in this paper (the approximation frequency segment is0.01,1 rad/s, the sampling frequency point is 40, and the order of approximation expression is taken as 6), the Bode diagram as shown in Figure 7 can be obtained.Figure 7
Comparison of the Muir method with the present method to approximates0.5.As can be seen from the Bode diagram, the amplitude and phase characteristic curves of the proposed method almost coincide with the ideal fractional differential operator in the approximation frequency range. The amplitude and phase characteristics of the approximate expression obtained by the Muir recursive method are quite different from the real values, and the high frequency amplitude frequency characteristics are seriously distorted. This method is superior to the Muir recursive method.


## 7.1. Pure Fractional Differential Operators

The expression of fractional differential operator isGs=sγ0<γ<1. In engineering, the Oustaloup filter is often used to approximate the fractional calculus operator in a given frequency band, assuming that the frequency range of interest is ωb,ωh; then, a set of broken lines can be used to approximate the linear characteristics of fractional calculus. Based on this idea, French scholar Professor Oustaloup [12] put forward the design method of CRONE filter. All these broken lines are generated alternately by zeros and poles of integer order so that the slope of the amplitude frequency characteristic asymptote changes alternately between 0 dB/dec and −20 dB/dec, and the frequency-domain response itself will be very close to an oblique line.The standard form of Oustaloup filter is shown in formula (29):(32)Gs=K∏k=1Ns+ωk′s+ωk,where k = 1,2, …, N, and zeros ωk′, poles ωk, and gain K can be calculated as follows:(33)ωk′=ωbωu2k−1−γ/N,(34)ωk=ωbωu2k−1+γ/N,(35)K=ωhγ,(36)ωu=ωhωb,where ωu is called unit gain frequency or transition frequency. Assuming that the frequency band of approximation is 10,100 rad/s, the value of γ is 0.5, and the number of iterations is 5; then, using the Oustaloup filter method [12], the approximation function can be obtained as follows:(37)GOustaloups=10s5+1726s4+107800s3+3039000s2+38650000s+177800000s5+217.4s4+17090s3+606300s2+9709000s+56230000.Using the quadratic programming method based on frequency-domain error, the frequency sample point is taken as 40, and the new approximation function is obtained as follows:(38)Gnews=0.05053s6+20.08s5+1668s4+33730s3+100000s2+0.001s+0.0010.001s6+1.271s5+217.7s4+9057s3+80130s2+35630s+100000.The amplitude phase characteristic curve is obtained by MATLAB simulation.It can be seen from Figure2 that the quadratic programming method based on the frequency-domain error can approach the pure fractional differential operator well in the specified frequency band. For the approximation of phase angle, it can be seen from Table 1 that the approximation effect of the quadratic programming method based on the frequency-domain error is better than that of the Oustaloup method, and it is closer to the theoretical value.Figure 2
Ideal Bode graph ofs0.5 and Bode graph obtained by two approximation methods.According to the sample standard deviation, we can measure the dispersion of the sample data. The formula is shown in (36):(39)s=1N−1∑i=1Nxi−x¯2,where N represents the sample population, xi represents the sample, x¯ represents the sample mean. According to the results of sample standard deviation in Table 4, the proposed method has smaller standard deviation, more concentrated data, and better approximation effect. The mean value is also closer to the ideal value.Table 5
Phase angle of characteristic frequency in approximation frequency band of different rational implementation methods.
2 rad/s34 rad/s67 rad/s100 rad/sx¯ (sample mean)S (standard deviation)Ideal−45.0°−45.0°−45.0°−45.0°−45.0°0°Oustaloup−32.4°−35.0°−27.8°−22.2°−29.35°5.6199°This paper−45.1°−45.0°−45.0°−45.0°−45.025°0.05°With the increase of the number of sample points, the quadratic programming method based on frequency error can obtain the amplitude frequency characteristic curve which basically coincides with the theoretical curve. In general, the quadratic programming method based on frequency error is better to approximate pure fractional differential operators.The frequency-domain characteristics of fractional differential operators with different orders are compared as shown in Figure3 (assuming that the approximation frequency segment is 1,100 rad/s, the sampling frequency point is 40, the order of approximation expression is 6, and the order is 0.1, 0.3, 0.5, 0.7, and 0.9, respectively).Figure 3
Approximation effect of fractional differential operators with different characteristic orders.Based on the frequency-domain quadratic programming method, the fitting effect of fractional order operator is very good in the middle and high frequency band of approaching frequency band, and the fractional order differential operator of each characteristic order has no large amplitude fluctuation in the frequency band. When the order is smaller (close to 0), the effect of approximation is better, and the amplitude and phase characteristics almost have no oscillation.


## 7.2. Pure Fractional Integral Operator

The expression of fractional integral operator isG=s−γ (0<γ<1). Assuming that the approximation frequency segment is 1,100 rad/s and the sampling frequency is 40, the order of approximation expression is 6 and the order is −0.5. The amplitude phase characteristic curve is obtained (as shown in Figure 4) by using the Oustaloup method and quadratic programming method based on the frequency-domain error proposed in this paper.Figure 4
Comparison of the effect of different methods to approximate fractional integral operator.Both the Oustaloup method and the method proposed in this paper can approach the fractional integral operator well. The Oustaloup method cannot fit the amplitude characteristics of fractional integral operators well in the high frequency band. For the fitting of phase frequency characteristics, refer to the data comparison in Table5 (the ideal amplitude and phase characteristic curve of integral operator s−0.5 is selected as the reference value).From the data in Table5, we can see that the phase frequency characteristics obtained by using this method to approximate the fractional order integral operator are consistent with the ideal phase frequency characteristics in the range from 34 rad/s to100 rad/s. The mean and standard deviation of the proposed method are smaller, the data are more stable, and the mean value is closer to the ideal value, so the approximation method proposed in this paper can more accurately approximate the ideal value, and the approximation data are more stable. It can be seen from Figure 4 that the phase frequency characteristic curves obtained by this method coincide with the ideal phase frequency characteristic curves in the middle and high frequency range. A better approximation effect is obtained.The frequency-domain characteristics of different order fractional integral operators are compared as shown in Figure5 (assuming that the approximation frequency band is 1,100 rad/s, the sampling frequency point is 40, and the order of approximation expression is 6, the integral characteristic order is −0.1, −0.3, −0.5, −0.7, −0.9).Figure 5
Approximation effect of fractional integral operators with different characteristic orders.It can be seen from Figure5 that the new method proposed in this paper can well approximate the characteristic fractional integral operator without oscillation in the approximation frequency range. The smaller the absolute value of integral order is, the better the approximation effect is.


## 7.3. Fractional Order System Model

Some common fractional order system models are simulated. For example, the fractional order model of PMSM is given as [19](40)Gs=9098.5417s1.7781+397316.2297s0.8371+1393070.1646s3.7915+127.38s2.8834+9995.678s1.8834.Original system cutoff frequencyωc = 40.8 rad/s. The constant parameters in this simulation are the approximation frequency (1,100 rad/s), the sampling frequency point (40), and the order of approximation expression (6). The results are shown in Figure 6.Figure 6
Approximation of ideal fractional order model by different methods.It can be seen from the amplitude frequency characteristic curve in Figure6 that the amplitude frequency characteristics of the system obtained by the two approximation methods are close. For the phase frequency characteristics obtained by the two approximation methods, in the second half of the approximation frequency band, the quadratic programming method based on the frequency-domain error can better fit the ideal phase frequency characteristics of the system than the Oustaloup method.


## 7.4. Comparison of New Method and Direct Discretization Method

The Muir recursive method is a direct discretization method commonly used in engineering. This method was first applied to geophysical scientific data processing and oil exploitation [10]. Without losing generality, we assume that the order of calculus operator is γ∈−1,1. For simplicity, we only give the recursive formula when γ is positive:(41)ωz−1r=2Tr1−z−11+z−1r=2Trlimn⟶∞Anz−1,rAnz−1,−r,where(42)A0z−1,r=1,(43)Anz−1,r=An−1z−1,r−cnznAn−1z,r,(44)cn=rn,nis odd,0,nis even.MATLAB toolbox was used, and the result is obtained as(45)sr≈2TrAnz−1,rAnz−1,−r.When the sampling frequencyT is 0.001, the number of iterations is 3. The fractional differential operator s0.5 can be directly discretized by the Muir recursive method in the following form [10]:(46)Gz=44.72z3−22.36z2+3.727z−7.454z3+0.5z2+0.8333z+0.1667.Compared with the method proposed in this paper (the approximation frequency segment is0.01,1 rad/s, the sampling frequency point is 40, and the order of approximation expression is taken as 6), the Bode diagram as shown in Figure 7 can be obtained.Figure 7
Comparison of the Muir method with the present method to approximates0.5.As can be seen from the Bode diagram, the amplitude and phase characteristic curves of the proposed method almost coincide with the ideal fractional differential operator in the approximation frequency range. The amplitude and phase characteristics of the approximate expression obtained by the Muir recursive method are quite different from the real values, and the high frequency amplitude frequency characteristics are seriously distorted. This method is superior to the Muir recursive method.


## 8. Conclusion

In this paper, a new rational approximation method of fractional calculus operator is proposed. The main idea of this method is to make the difference between the frequency characteristic of the approximation transfer function and the ideal fractional calculus operator in the frequency domain and then make quadratic programming for the difference expression to obtain the optimal rational approximation transfer function coefficients. Compared with the previous methods, the method proposed in this paper has obvious improvement in approximation accuracy, precision, and algorithm encapsulation and can arbitrarily select the approximation frequency band, which makes the algorithm highly flexible and has great advantages. When the parameters of the approximation expression fluctuate by 10%, the phase angle change of the system is less than 3.473°, which indicates that the system has strong robustness. By using different methods to simulate and compare the fractional calculus operator and model, it is concluded that this method can effectively fit the frequency-domain characteristics of fractional calculus operator and fractional order model in any specified frequency range. However, for the fractional calculus operator whose fractional order is close to 1, there is a large error with the ideal value in the initial frequency range, so the initial constraint condition can be added to reduce the error. When the order is close to 0, the amplitude phase characteristic curve almost coincides with the amplitude phase characteristic curve of ideal fractional calculus operator, and a good approximation effect is obtained. This method can only approximate the fractional order system with pure fractional calculus operator, but it is difficult to deal with the approximate irrational fractional order system, such as single-fractal system Hs=1/1+s/pTm. But it is not impossible to deal with it. We can expand the irrational term into a system with only pure fractional calculus operator, and then we can use the method in this paper to approximate it.



---

*Source: 6646718-2021-02-12.xml*