Research Article | Open Access

# A Robust Direct Parameter Identification of Exponentially Damped Low-Frequency Oscillation in Power Systems

**Academic Editor:**Pawel Malinowski

#### Abstract

Low-frequency oscillations in power systems can be modeled as an exponentially damped sinusoid (EDS) signal. Its frequency, damping factor, and amplitude are identified by the robust algorithm proposed in this paper. Under the condition of no noise, the exponentially convergent property of the proposed identification is proved by the use of time scale change, variable transformation, slow integral manifold, averaging method, and Lyapunov stability theorem in sequence. Under the condition of bounded additive noise, the antinoise performance of the identification of each parameter is investigated by the perturbed system theorem and error synthesis principle. The robustness of the proposed method is embodied in the following aspects: the exponential convergence for EDS signal with a wide range of frequency, especially with a rather low frequency; the boundary values of identification errors resulting from high-frequency sinusoidal noise of both frequency and damping factor can be adjusted by tuning the design parameters; and the different effects of the four design parameters on tracking performance and antinoise performance of each parameter identification. Simulation results demonstrate the performance of the algorithm and validate the conclusions.

#### 1. Introduction

Widespread used renewable and sustainable energy sources bring modern power systems electromechanical oscillations, synchronous, and subsynchronous oscillations [1–6]. The low-frequency oscillations (LFO) may result in unstable and unsecure operation of power systems, so the real-time identification of frequency, amplitude, and damping factor are still essential in recent years [7–10].

In the parameter identification filed, there are mainly two challenges, model selection, and parameter estimation [11, 12]. When number of sinusoids, i.e., model order, is unknown, some techniques, such as the Bayesian method, have been developed to identify the model order [13].

In a transient state of power systems, the electromechanical oscillations usually have amplitudes decaying as time and low oscillating frequency in range of 0.05-2.0 Hz [9, 10]. Because the oscillating frequency is much lower than the fundamental frequency of power systems, 50 Hz or 60 Hz, the LFOs in power systems are easily extracted from the measured voltage signals or current signals by low-pass filters. Furthermore, the most energy of oscillation usually concentrates on one frequency. So, the LFO can be modelled as low frequency exponentially damped sinusoid (EDS) signal. The key technique to identify the parameters of LFOs is now the parameter identification of one exponentially damped low-frequency sinusoid.

When number of sinusoids, i.e., model order, is known, numerous studies have been developed in the signal processing field. Singular value decomposition (SVD) has been used to estimate the parameters of exponentially damped sinusoidal signals in noise [14]. In SVD method, many sampled data in a rather long period of time must be required, and heavy computation burden is implemented with the number of sinusoid components more than one.

A type of methods to estimate the parameters of EDS signals are based on integration, such as a method combining Hankel singular value decomposition (HSVD) with the extended complex Kalman filter (ECKF) [2], an algorithm employing the DFT coefficients around the peaks of the components [3], a moving window discrete Fourier transform filter (MWDFT) in cooperation with a frequency-locked loop [8], and two novel interpolation iterative estimators based on DFT coefficients [10]. In these integral methods, the large computation restrains the online identification in a more wider engineering.

Another large class of techniques has been presented based on the differential model instead of the integration model. These differential methods employ the concept of frequency estimator.

In [6], the existing concept of the adaptive notch filter and second-order generalized integrator was modified for the estimation of the frequency, the magnitude, and the damping of low-frequency electromechanical oscillations with a direct component in an interconnected power system.

Two types of frequency estimators on the basis of adaptive notch filter (ANF), the normalized estimator and the nonnormalized estimator, were investigated in [15]. The actual value of signal amplitude has little effect on the convergence speed in the normalized methods while it has a great influence on the convergence speed in the nonnormalized methods. The nonnormalized version was extended in [16] to estimate frequencies and amplitudes of multisinusoidal components. The normalized version was improved in [17] in order to gain exponential convergence property in a wide range of frequency, especially in a rather low frequency. It was then applied in [18] to cancel sinusoidal disturbance with unknown low frequency of linear time-invariant systems.

Based on the adaptive identifier proposed in [19], an indirect algorithm was proposed in [20] for identifying the parameters of the EDS signal. The variables in this type of indirect algorithms converge to the value combining of the square of frequency and square of damping factor, which is inconvenient to provide the instantaneous estimate value of EDS parameters.

The concept of enhanced phase-locked loop (EPLL) proposed in [21], a nonnormalized frequency estimation, was normalized and extended [9] for the estimation of EDS parameters. In the EPLL-based methods, it is inconvenient to obtain the instantaneous value of input EDS signal.

Based on the adaptive internal model (AIM) controller, a normalized frequency estimator proposed in [22] was extended in [23] to identify the frequency and the damping factor of the EDS signal and then was applied in [24] to the power system for damping of electromechanical oscillations. In addition to the frequency and the damping factor, the AIM-based algorithms directly provide the instantaneous estimate values of both the in-phase component and the quadrature-phase component of the EDS signals.

Compared with integral methods, the differential methods need a smaller number of data in a shorter period time and less computation burden.

As been addressed in [17], the AIM methods may lost convergence property with its design parameters unchanged when the frequency value of EDS signal changes from some number to a rather smaller number. So, it is significant to bring some improvement on transient performance and antinoise performance of the parameter identification of LFOs modeled as EDS signals.

Motivated by [15, 17, 23], the robustness of the new method proposed here includes exponential convergence property in a wide range of frequency under no noise and adjustable boundary of steady-state identification errors under bounded additive noise.

The rest of the paper is arranged as follows: The framework of the new algorithm is described in Section 2. The identification characteristics are investigated in Section 3. The simulation study is carried out in Section 4. Section 5 concludes the whole paper.

#### 2. Materials and Methods

##### 2.1. Framework of New Algorithm

A low-pass filtered LFO in the power system can be represented as EDS signal in time : where denotes exponentially damped oscillation and represents the additive measure noise. The amplitude , the damping factor , and the frequency are all unknown constant real positive numbers, and the initial phase .

Motivated by the normalized updating law proposed in [15] and the robust frequency estimator proposed in [17], we change the AIM-based identification method proposed in [23] into the new parameter identification and write as the following close-loop dynamical system that consists of three differential equations and three algebraic equations

In the equations abovementioned, , , , and are four positive real design parameters to adjust the performance of identification, , , , and are four variables, , , and represent the identification value of EDS signal , and the signal tracking error and the identification value of amplitude , respectively. The two variables, and , are just right the identification values of the frequency and the damping factor , respectively.

We may assume that the four variables have initial values expressed as

Generally, all the parameters of the EDS signal have finite values. So, the four differential variables are bounded for and can be assumed in a compact set defined as where all the boundary values, , , , and , can be experimentally preset.

#### 3. Results and Discussion

##### 3.1. Characteristics of the Identification

In this section, we characterize the new identification from two aspects: transient convergence property and steady-state antinoise performance.

###### 3.1.1. Convergence Property Analysis

The convergence property of the proposed identification algorithm can be summarized as the following theorem.

Theorem 1. *When the measure noise , there exist positive numbers , , , and , such that if , , , and , then all the identification values exponentially converge to their actual values, i.e.,
*

*Proof. *We use time-scale change, variable transformation, slow integral manifolds, and averaging method in sequence to prove the asymptotic convergence property of the proposed identification algorithm.

Firstly, we define the time scale change as
and rewrite the EDS signal in time as
where
Put (6) into (2), (3), (4), let then we rewrite the state-space representation in time as
where the time in the state variables is omitted for convenience of writing, the same below.

Then, we define the variable transformation as
Let . After substituting (17) into (14) (15) (16), we obtain a nonlinear dynamic system consisting of two coupled 2-dimension differential equations:
where
The definite domain of the transformed variables, corresponding to the compact set , is expressed as
Now, the truth of (10) is equivalent to that the nonlinear system consisting of (18) and (19) has the asymptotical convergence property on the compact set expressed as

Secondly, we prove (22) true by the schedule similar to that done in [16, 17].

The dynamic system including (18) and (19) just accords with standard form of the system addressed in [25]. So, we use slow integral manifold proposed in [25] to decouple the updating law (19) from the state equation (18).

Considering frozen parameter method, let , the state equation (18) becomes a linear time-invariant (LTI) system which constant coefficient matrix has two eigenvalues expressed as

The two eigenvalues have both negative real parts if the following condition in is always satisfied.

So return to time domain of it must be satisfied that

This declares the existence of and such that and should be ensured.

For the LTI system (18) when , the negative real parts of also reveal that its transient response to the initial value of vanishes and its steady-state response forced by input signal sin t can be calculated by use of frequency characteristics.

Let . Under the condition of initial relaxation, the frequency characteristics of (18) from the input signal, i.e., , to state variable vector can be written as

We represent the complex denominator in the form of amplitude-frequency characteristics and phase-frequency characteristics expressed as where.

Now, the steady-state response of LTI system (18) is a 2*π*- period orbit represented by

Then, the steady-state identification error is

It can be verified that all the following items are continuous, bounded, and Lipschitzian with respect to both and in . This ensues both Assumption 2.2 and Assumption 3.1 of [25]. From (25), Assumption 2.1 of [25] has been ensured by the negative real parts of .

So, according to Theorem 3.1 of [25], there exists and such that if and then system combining (18) and (19) has -family slow integral manifolds.

On these manifolds, the state variable vector can be expressed as

Then, the identification error is represented by

Now, the parameter updating law (19) can be decoupled from the state equation (18) and expressed as two coupled almost periodical dynamic systems with respect to .

Thirdly, we use the averaging method to prove the convergence property of (34) and (35).

From (29), (30), (32), and (33), we let and get where

Similarly, we also get

The integration of the following periodical function is zero.

According to the integral formula, we get where

According to Section 10.4 of [26], we get

Thus, the averaged system of updating law (18) is

Considering , the point is the unique equilibrium point of the averaged system (43). Let Lyapunov candidate function be where is a sufficiently large positive real number. Its derivative along the trajectory of the system (43) is

When , then and obviously. When , it can be discussed as follows: (1)If , then and (2)If , then we can let be enough large such that resulting from

So, the derivative function is negatively defined and the point is the unique globally stable equilibrium point of the averaged system (43).

Invoking Theorem 10.4 of [20], for sufficiently small positive real number , the two variables and in (34) and (35) are asymptotically convergent to their equilibrium point,

respectively, i.e., .

In steady state, put and into (29), we get and then prove (22) true.

Now, recalling the time scale change (11), the variable transformation (17), and the condition (25), the theorem and the conclusion of (10) have been proved true.

###### 3.1.2. Noise Characteristics Analysis

Let the high-frequency measure noise in the signal (1) be bounded and represented by where the amplitude , the frequency , and the phase angle are constants.

The noise characteristics of the new identification can be summarized as

Theorem 2. *Under high-frequency noise of (46), if the design parameters, , , , and satisfy the conditions for Theorem 1, then the steady-state identification errors will be bounded by
where and are positive real constants in .*

*Proof. *By use of the time scale change (11), the noise (46) is rewritten in time as
Firstly, consider , let denote the state variable vector of the LTI system (14) forced only by noise instead of EDS signal , we get
The frequency characteristics of the system (51) from to state variable vector can be written as
So the steady-state response of the LTI system (51) is
where
In this instance, the identification error resulting from is rewritten as
where has some appropriate value and
Let denote the steady-state response of the LTI system (14) forced only by the EDS signal . According to the superposition principle of LTI systems, the steady-state response of the LTI system (14) resulting from both the EDS signal and the noise is represented by
Secondly, the decoupled updating laws (15) and (16) can be written as the following perturbed system:
Theorem 1 abovementioned shows that the nominal system
is an exponentially stable system.

The perturbation items are bounded by
According to Lemma 9.2 of [26], there exit boundary values of and such that (47) and (48) hold true.

Thirdly, from (7), by error synthesis principle, we know
where

After the trajectory goes into the small neighborhood of the equilibrium point, we get and . If and , then , .

Now, the Theorem 2 has been consequently proved true.

###### 3.1.3. Performance Evaluation.

For the measured signal (1), the identification method proposed in [17] with two positive design parameters, and can be written in time as

With time scale change (11) applied, it is written in time as

After time scale change (11), the identification algorithm proposed in [6] can be rewritten in time as where , , and are four positive design parameters.

*Remark 3. *When has a rather large value, both the algorithm (64) and the new algorithm exemplify exponential convergence property. When has a rather small value, the explicit reciprocal in (65) may bring the algorithm (64) unconvergence property. Contrarily, the reciprocal in (15) and (16) is compensated by the frequency estimate value , which brings the new method better robustness for low-frequency signals.

When being distorted by the noise (50), the perturbation terms of the algorithm (65) similar to (60) and (61) are

In a small neighborhood of the equilibrium point, , so the boundary values of identification errors of frequency and damping factor in (65) can be represented as

*Remark 4. *When the signal is distorted by high-frequency sinusoidal noise with constant amplitude, the identification errors of frequency and damping factor in the algorithm (65) increase exponentially over time, while that in the new method have both bounded values governed by design parameter .

*Remark 5. *For the identification of both frequency and damping factor in the new method, a greater value of the design parameter brings smaller steady-state identification errors under noise and slower transient convergence speeds simultaneously.

*Remark 6. *The state variable in (66) denotes the estimated value of the direct current component. Furthermore, compare (65) with (66), we know that the same conclusions as Remark 3 and Remark 4 can be drawn for the algorithm in [6].

#### 4. Simulation Results

In this section, the performance of the proposed algorithm is exemplified by simulation results employing Matlab/Simulink.

##### 4.1. Performance Evaluation

Firstly, we let the EDS signal be .

At , the frequency changes from down to , and the damping factor changes from 0.03 down to 0.01. All the time, the design parameters of our new identification are set to , , , and , the design parameters of the Browns algorithm proposed in [23] are and . The design parameters of the Mojiri’s algorithm proposed in [6] are , , , and . The simulation results of the three algorithms are illustrated in Figure 1 for comparison.

**(a)**

**(b)**

In Figure 1, the two black dashed lines of the algorithm in [6] and the two red solid lines of the algorithm in [23] show that both the frequency and the damping factor are accurately estimated when and neither the frequency nor the damping factor are correctly identified when with the same values of their design parameters. Contrarily, the two blue solid curves, results from the new algorithm, exemplify that both frequency and damping factor are exactly identified no matter what value is. This validates the conclusion of Remark 3: the convergence robustness of the new identification algorithm is better than that of the algorithm in [23] and the algorithm in [6].

Secondly, let the EDS signal be distorted by sinusoidal noise . We set , , , and in the new algorithm, and let and in algorithm of [23], and let , , , and in algorithm in [16]. The simulation results of the three algorithms are illustrated in Figure 2.

In Figure 2, the two red curves from the algorithm in [23] and the two black curves from the algorithm in [6] show that the identification errors of both frequency and damping factor, respectively, increase exponentially over time. In the two blue curves from the new algorithm, the boundary value of frequency identification error is ±0.0004 and the boundary value of identification error of the damping factor is ±0.0005. This fact verifies the conclusion of Remark 4: the antinoise performance of the new algorithm is far better than that of the algorithm in [23].

##### 4.2. Amplitude Identification

Thirdly, the amplitude identification proposed in this paper is exemplified in Figure 3. We let the EDS signal be , let the sinusoidal noise be , and let the design parameters be , , , and .

**(a)**

**(b)**

In Figure 3, the subplot (a) shows that the proposed algorithm can exactly identify the amplitude of the EDS signal under no noise while the subplot (b) declares that the amplitude identification error resulting from high-frequency sinusoidal noise grows exponentially over time. These simulation results verify the conclusion of (49).

##### 4.3. Effects of Design Parameters

Finally, we exemplify the effect of the design parameter on tracking performance and antinoise performance by simulation examples with different value of . We let the EDS signal be and let the three design parameters be , , and .

In each subplot of Figure 4, from the same initial values of , the convergence speed of the red dashed curves corresponding to is slower than that of the blue solid curves corresponding to , while faster than that of the black dashed curves corresponding to . This validates one conclusion of Remark 5: the larger value of the design parameter brings slower convergence speed to the identification of frequency, damping factor, and amplitude.

In Figure 5, let the EDS signal in Figure 4 be distorted by the noise . In the three frequency identification errors illustrated in the left three subplots, the boundary values of the blue curve, the red curve, and the black curve, corresponding to , , and are ±0.0041, ±0.0025, and ±0.0009, respectively.

In the three damping factor identification errors illustrated in the right three subplots, the boundary values of the blue curve, the red curve, and the black curve, corresponding to ,, and , are ±0.0022, ±0.0012, and ±0.0005, respectively. So, a larger value of the design parameter brings smaller steady-state identification error resulting from sinusoidal noise to the identification of frequency and damping factor. This verifies the other conclusion of Remark 5.

The design parameters should be tuned according to the compromise between transient performance and steady-state performance. From (43)–(49) as well as Figures 4 and 5, big values of the three parameters, , , and, bring faster-tracking performance and worse antinoise performance while big value of results in slower tracking speed and better antinoise performance.

#### 5. Conclusion

A robust algorithm is proposed for asymptotically identifying the frequency, the damping factor, and the amplitude of a LFO modeled as an exponentially damped sinusoidal signal. The transient convergence property under no noise and the boundary of the steady-state identification errors under sinusoidal noise are exhibited. The new algorithm exemplifies good convergence property for a rather wide range of frequency signal, especially for low-frequency signal, as well as bounded identification errors being adjusted by values of design parameters under bounded noise. The steady-state identification errors of frequency and damping factors are bounded with boundary values adjustable by the value of design parameters. The design parameters have different effects on transient performance and steady-state performance of the identification of frequency, damping factor, and amplitude, respectively. The proposed amplitude identification will be improved in antinoise performance in expectation.

#### Conflicts of Interest

The author declares that there are no conflicts of interest.

#### Acknowledgments

This work was supported in part by the Open Fund of State Key Laboratory of Operation and Control of Renewable Energy & Storage Systems (China Electric Power Research Institute) (No. EPRI4124-190856).

#### References

- H. Ye, Y. Liu, and X. Niu, “Low frequency oscillation analysis and damping based on Prony method and sparse eigenvalue technique,” in
*2006 IEEE International Conference on Networking, Sensing and Control*, pp. 1006–1010, Ft. Lauderdale, FL, USA, 2006. View at: Publisher Site | Google Scholar - J. Zhang, A. K. Swain, and S. K. Nguang, “Parameter estimation of exponentially damped sinusoids using HSVD based extended complex Kalman filter,” in
*TENCON 2008 - 2008 IEEE Region 10 Conference*, pp. 1–6, Hyderabad, India, 2008. View at: Publisher Site | Google Scholar - D. Agrez, “A frequency domain procedure for estimation of the exponentially damped sinusoids,” in
*2009 IEEE Intrumentation and Measurement Technology Conference*, pp. 1321–1326, Singapore, 2009. View at: Publisher Site | Google Scholar - J. Rueda, C. Juarez, and I. Erlich, “Wavelet-based analysis of power system low-frequency electromechanical oscillations,”
*IEEE Transactions on Power Systems*, vol. 26, no. 3, pp. 1733–1743, 2011. View at: Publisher Site | Google Scholar - D. Lauria and C. Pisani, “On Hilbert transform methods for low frequency oscillations detection,”
*IET Generation, Transmission & Distribution*, vol. 8, no. 6, pp. 1061–1074, 2014. View at: Publisher Site | Google Scholar - M. Mansouri, M. Mojiri, M. A. Ghadiri-Modarres, and M. Karimi Ghartemani, “Estimation of electromechanical oscillations from phasor measurements using second-order generalized integrator,”
*IEEE Transactions on Instrumentation and Measurement*, vol. 64, no. 4, pp. 943–950, 2015. View at: Publisher Site | Google Scholar - A. Q. Zhang, L. L. Zhang, M. S. Li, and Q. H. Wu, “Identification of dominant low frequency oscillation modes based on blind source separation,”
*IEEE Transactions on Power Systems*, vol. 32, no. 6, pp. 4774–4782, 2017. View at: Publisher Site | Google Scholar - S. Tomar and P. Sumathi, “Amplitude and frequency estimation of exponentially decaying sinusoids,”
*IEEE Transactions on Instrumentation and Measurement*, vol. 67, no. 1, pp. 229–237, 2018. View at: Publisher Site | Google Scholar - H. Zamani, M. Karimi-Ghartemani, and M. Mojiri, “Analysis of power system oscillations from PMU data using an EPLL-based approach,”
*IEEE Transactions on Instrumentation and Measurement*, vol. 67, no. 2, pp. 307–316, 2018. View at: Publisher Site | Google Scholar - A. Serbes, “Fast and efficient sinusoidal frequency estimation by using the DFT coefficients,”
*IEEE Transactions on Communications*, vol. 67, no. 3, pp. 2333–2342, 2019. View at: Publisher Site | Google Scholar - G. L. Bretthorst, “Bayesian analysis. III. Applications to NMR signal detection, model selection, and parameter estimation,”
*Journal of Magnetic Resonance (1969)*, vol. 88, no. 3, pp. 571–595, 1990. View at: Publisher Site | Google Scholar - C. Andrieu, P. M. Djuri, and A. Doucet, “Model selection by MCMC computation,”
*Signal Processing*, vol. 81, no. 1, pp. 19–37, 2001. View at: Publisher Site | Google Scholar - K. Xu, G. Marrelec, S. Bernard, and Q. Grimal, “Lorentzian-model-based Bayesian analysis for automated estimation of attenuated resonance Spectrum,”
*IEEE Transactions on Signal Processing*, vol. 67, no. 1, pp. 4–16, 2019. View at: Publisher Site | Google Scholar - R. Kumaresan and D. Tufts, “Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise,”
*IEEE Transactions on Acoustics, Speech, and Signal Processing*, vol. 30, no. 6, pp. 833–840, 1982. View at: Publisher Site | Google Scholar - L. Hsu, R. Ortega, and G. Damm, “A globally convergent frequency estimator,”
*IEEE Transactions on Automatic Control*, vol. 44, no. 4, pp. 698–713, 1999. View at: Publisher Site | Google Scholar - Z. Chu, M. Ding, S. Du, and X. Feng, “Exponential stability, semistability, and boundedness of a multi-ANF system,”
*IEEE Transactions on Circuits and Systems I: Regular Papers*, vol. 58, no. 2, pp. 326–335, 2011. View at: Publisher Site | Google Scholar - Z. Chu, C. Sheng, M. Zhu, B. Chen, and H. Li, “A robust adaptive identification of sinusoidal signal with unknown frequency,”
*IEEE Transactions on Circuits and Systems II: Express Briefs*, vol. 66, no. 9, pp. 1562–1566, 2019. View at: Publisher Site | Google Scholar - Z. Chu, W. Zhang, M. Zhu, X. Dong, and B. Chen, “A robust adaptive cancellation of unknown sinusoidal disturbance,”
*IEEE Transactions on Circuits and Systems II: Express Briefs*, vol. 67, no. 4, pp. 705–709, 2020. View at: Publisher Site | Google Scholar - X. Xia, “Global frequency estimation using adaptive identifiers,”
*IEEE Transactions on Automatic Control*, vol. 47, no. 7, pp. 1188–1193, 2002. View at: Publisher Site | Google Scholar - M. A. Ghadiri-Modarres and M. Mojiri, “Control of oscillatory damped disturbances-an indirect approach,” in
*The 2nd International Conference on Control, Instrumentation and Automation*, pp. 1154–1159, Shiraz, Iran, 2011. View at: Publisher Site | Google Scholar - M. Karimi-Ghartemani and A. K. Ziarani, “Performance characterization of a non-linear system as both an adaptive notch filter and a phase-locked loop,”
*International Journal of Adaptive Control Signal Process*, vol. 18, no. 1, pp. 23–53, 2004. View at: Publisher Site | Google Scholar - Q. Zhang and L. J. Brown, “Noise analysis of an algorithm for uncertain frequency identification,”
*IEEE Transactions on Automatic Control*, vol. 51, no. 1, pp. 103–110, 2006. View at: Publisher Site | Google Scholar - J. Lu and L. J. Brown, “Internal model principle-based control of exponentially damped sinusoids,”
*International Journal of Adaptive Control Signal Process*, vol. 24, no. 3, pp. 219–232, 2010. View at: Publisher Site | Google Scholar - M. Mansouri, A. M. Knight, Z. Moradi-shahrbabak, and M. Mojiri, “Damping of electromechanical oscillations based on the internal model principle,” in
*2017 IEEE International Electric Machines and Drives Conference (IEMDC)*, pp. 1–6, Miami, FL, USA, 2017. View at: Publisher Site | Google Scholar - B. Riedle and P. Kokotovic, “Integral manifolds of slow adaptation,”
*IEEE Transactions on Automatic Control*, vol. 31, no. 4, pp. 316–324, 1986. View at: Publisher Site | Google Scholar - H. K. Khalil,
*Nonlinear System*, Macmillan, New York, 3rd edition, 2007.

#### Copyright

Copyright © 2020 Zhaobi Chu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.