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 operator 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.
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 into n 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]where 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 function has n+1 derivative in the interval . For any real number (the lower limit of integration at the initial time of 0 is 0), the α-order calculus of function is defined aswhere represents the number of approximate recurrence terms, in whichis 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 number , , the RL fractional order is defined aswhere is a gamma function.
RL fractional integral is defined as
2.1.3. Caputo Definition of Fractional Calculus
Caputo fractional calculus is defined as
Among them, . 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, represents the best approximate rational expression of , and is the error between rational approximation function and ideal fractional calculus operator in frequency domain. means to use some optimization algorithm for error , such as genetic algorithm, linear programming, model reduction method, and so on. represents the optimized frequency-domain error function. The parameter of is obtained by calculating the partial derivative of each parameter of . Different from other methods, in this paper, (fractional transfer function model to be optimized) is changed into the frequency-domain model of fractional order calculus operator (), reducing the time and complexity of the program.
For the fractional calculus operator to be approximated, it can be expanded into the following form:
After optimization by some optimization algorithm, the rational approximation function is obtained. Its form is shown in the following formula:where are the coefficient parameters of approximation expression . 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 . Suppose that characteristic frequency points are selected in the approximation frequency range. Among them,and at these points, the transfer function is approximated:
For each frequency point :and at each characteristic frequency point, the frequency error is . Then,and .
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:and the frequency-domain error is
Let
In order to get the minimum value of the frequency-domain error, needs to find the partial derivatives of parameters and , respectively. Let the partial derivative be zero and solve the system of equations; we can get the minimum value of .
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:and then
Let , and , where
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,and (19)–(23) and (25) can be arranged into the following matrix form:
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 operator , 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 ; 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]:where Q is a symmetric matrix, is a linear constraint condition, is equality constraint, and is independent variable definition domain.
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 vector 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 is . Next, we analyze the sensitivity caused by the change of parameters in the vector.
Firstly, the value of nonsensitive quantity is determined: assuming that the frequency band of approximation is 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
At the initial approximation frequency band, the rational approximation function can be reduced to
The main parameters that affect the system (the quantity that determines the sensitivity) are .
When 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.
Investigate the sensitivity index when each parameter changes; the following can be obtained:
Therefore, we can derive the sensitivity when a single weight parameter changes, as shown in Table 3.
It can be seen from Table 3 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.
6. Algorithm Complexity Analysis
6.1. Time Complexity
The time complexity of the proposed algorithm is dominated by two loop statements, which are for i = 1:n A(:,:,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]; end and for i = 1:n C(:,:,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)]; end
The number of execution steps of the two loop statements is n 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 size n, 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 is . 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 ; 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):where k = 1,2, …, N, and zeros , poles , and gain can be calculated as follows:where is called unit gain frequency or transition frequency. Assuming that the frequency band of approximation is 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:
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:
The amplitude phase characteristic curve is obtained by MATLAB simulation.
It can be seen from Figure 2 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.
According to the sample standard deviation, we can measure the dispersion of the sample data. The formula is shown in (36):where N represents the sample population, represents the sample, 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.
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 Figure 3 (assuming that the approximation frequency segment is 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).
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 is (). Assuming that the approximation frequency segment is 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.
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 Table 5 (the ideal amplitude and phase characteristic curve of integral operator is selected as the reference value).
From the data in Table 5, 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 Figure 5 (assuming that the approximation frequency band is 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).
It can be seen from Figure 5 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]
Original system cutoff frequency = 40.8 rad/s. The constant parameters in this simulation are the approximation frequency ( rad/s), the sampling frequency point (40), and the order of approximation expression (6). The results are shown in Figure 6.
It can be seen from the amplitude frequency characteristic curve in Figure 6 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 . For simplicity, we only give the recursive formula when is positive:where
MATLAB toolbox was used, and the result is obtained as
When the sampling frequency T is 0.001, the number of iterations is 3. The fractional differential operator can be directly discretized by the Muir recursive method in the following form [10]:
Compared with the method proposed in this paper (the approximation frequency segment is 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.
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 . 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.
Data Availability
The data used to support the findings of this study are available from the first author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.