Systematic Control Design for Adaptive Rejection of Narrow-Band Disturbances Based on a Novel Resonant Generalized Extended State Observer

Abstract

In this paper, the adaptive rejection of multiple narrow-band disturbances with time-varying frequencies is investigated. A systematic control approach based on a resonant generalized extended state observer (RGESO) is proposed. The frequency response of the RGESO is shaped into a desired form to reject the total disturbance, including the wide-band uncertainties and the narrow-band disturbances. The state-space realization of the RGESO is also provided to estimate both the states and the total disturbance. To handle the time-varying frequencies of the narrow-band disturbances, a direct frequency estimation approach based on the Levenberg–Marquardt method and an online supervision strategy is presented to improve convergence and reliability, and the estimations are then used to adapt the RGESO parameters. Parameter tuning and stability robustness analysis methods are also developed for the proposed scheme to facilitate its application. Finally, extensive simulations are conducted to demonstrate the effectiveness of our method.

1. Introduction

Disturbance rejection is a fundamental problem in control theory and practice. There are narrow-band disturbances in a surprisingly large number of applications, such as the hard disk drive [1,2], suspension system [3,4], and power active filter [5,6]. In some applications, the frequencies of the disturbances are unknown and time-varying. Thus, an adaptive narrow-band disturbance rejection design, which can accommodate to the frequency change and achieve high performance, is urgently needed.

The adaptive narrow-band disturbance rejection has been extensively studied, and many control methods have been proposed. The phase-locked loop (PLL)-based adaptive feedforward cancellation (AFC) [7,8], the internal model principle (IMP)-based adaptive regulation [3,9], and the disturbance observer-based control (DOBC) [1,4,5,10] are the main techniques. In the PLL method, both the magnitude and frequency can be estimated in real time. However, the frequency responses of the plant at all the concerned frequencies need to be known in advance for parameter adaptation, which limits its application to some complicated plants whose frequency responses are not available. The IMP is also a powerful design approach to reject narrow-band disturbances. According to the IMP, the controller should include the disturbance model and have resonant poles determined by the disturbance frequencies. Many adaptive regulation methods based on the IMP have been proposed under the assumption that the plants can be precisely modeled. However, the effects of model uncertainties cannot be ignored for some plants and must be taken into account together with the narrow-band disturbances. Although many researchers have considered using the adaptive control theory to solve the unknown disturbances as well as the model uncertainties, the system may be destabilized in the presence of unmodeled dynamics [11]. The DOBC is a promising alternative since both the narrow-band disturbances, and the model uncertainties can be rejected in a unified and simple manner. Thus, we focus on the DOBC for the adaptive narrow-band disturbance rejection problem in this paper.

Among the various types of disturbance observers, the resonant observer (RO) is a fundamental technique for narrow-band disturbance estimation [5,12], whereas it cannot deal with model uncertainties. Different from the RO design, which is performed in the state space, the design of the Q-filter based disturbance observer (QBDO) is conducted in the frequency domain; although, it is model-dependent. In most of the literature on the QBDO design, only the narrow-band disturbance can be estimated [1,4], and the core idea is to modify the frequency response of the traditional low-pass Q-filter based on the disturbance model. Jia et al. [1] designed the Q-filter as an adaptive band-pass filter with the frequency estimation. Since the designs of disturbance rejection and frequency estimation algorithms are performed separately, this scheme is categorized as an indirect method. In contrast, the direct method involves adapting the Q-filter directly to minimize the effect of the narrow-band disturbance [4]. In addition to rejecting the narrow-band disturbance, the capability of attenuating the low-frequency wide-band uncertainty is retained [13,14]. A period time-delay element, which can be updated in real time, is cascaded to a standard low-pass Q-filter [13]. However, additional phase delays are introduced at the disturbance frequencies by the low-pass filter, which is undesired, especially for the high-frequency disturbances. To overcome this deficiency, a novel Q-filter was proposed [14]. The frequency response of this Q-filter is unity at the frequency of the narrow-band disturbance, and it retains the characteristic of a low-pass filter at other frequencies. However, it cannot provide the estimations of states.

In recent years, the extended state observer (ESO), which is the core of the active disturbance rejection control (ADRC) [15,16], has become a popular approach for disturbance rejection and has been applied successfully in many fields. Unlike the QBDO design, which requires a nominal model of the plant, cascaded integrators are selected as the nominal model and a total disturbance (including both the internal uncertainties and the external disturbances) concept is used in this scheme. Therefore, the ESO can be applied to complicated plants whose models are difficult to obtain [17]. Various enhanced forms have been proposed to address specific challenges in disturbance rejection. Considering that the standard ESO [16] is limited in performance when high-frequency measurement noise is fed into it, cascade-parallel ESOs [18,19] are proposed to maintain the good noise suppression of the cascaded ESO and overcome its weakness of poor disturbance rejection. To handle the unknown input gains, adaptive ESOs are developed to estimate the input gains, unmeasured states, and total disturbance simultaneously [20,21]. To further promote the anti-disturbance property, different nonlinear mechanisms are introduced to the ESO design. To improve the convergence speed and estimation accuracy, the generalized super-twisting technique is introduced to construct the finite-time-convergent ESOs [22,23], which ensure that the estimation error can converge to zero in finite time. Although the convergence time becomes finite, it is closely related to the initial condition and grows unboundedly with the increase in the initial error. To address this problem, fixed-time-convergent ESOs are designed [24,25]. Additionally, with regards to the peaking phenomenon of the ESO, fractional-order ESOs [26] and switching ESOs [27] are developed to improve the transient performance. Many researchers [28,29,30] have investigated the performance of the generalized extended state observer (GESO) on the narrow-band disturbance rejection, but it is assumed that the mth-order derivative of the disturbance is zero, which is obviously not reasonable for the sinusoidal disturbance. Instead of increasing the ESO order, Zheng et al. [31] cascaded a phase compensator, which can provide a phase lead around the frequency of the narrow-band disturbance, to a low-pass ESO to compensate for the undesired phase delay at the high frequency. In these methods [28,29,30,31], the ESO is still in the form of a low-pass filter. Therefore, a high observer bandwidth is required to guarantee the rejection performance of the high-frequency narrow-band disturbance, which can affect the robustness to unmodelled dynamics and measurement noises.

In view of the above-mentioned facts, we propose a systematic control design method for adaptive narrow-band disturbance rejection based on a resonant generalized extended state observer (RGESO) and a robust frequency estimation method. The main contributions are

(i) The frequency response of the RGESO is shaped so that both the low-frequency wide-band uncertainties and the high-frequency narrow-band disturbances can be simultaneously rejected. Therefore, the high observer bandwidth necessary for an ordinary ESO is avoided, which can greatly improve the system robustness.

(ii) The bandwidths for the wide-band uncertainty and the narrow-band disturbance in the RGESO are represented by various parameters whose physical meanings are explicit, which can make the parameter tuning process quite straightforward. Furthermore, both parameter optimization and stability robustness analysis methods are provided to facilitate its application.

(iii) A robust adaptive notch filter (ANF) is proposed to estimate the unknown and time-varying disturbance frequencies and update the parameters of the RGESO. A direct algorithm is designed to avoid the computationally complicated factorization. Additionally, a recursive version of the Levenberg–Marquardt (LM) method is used to enhance the convergence.

The remainder of the paper is organized as follows. The design process of the RGESO is presented in Section 2. In Section 3, the algorithm for direct frequency estimation is proposed. Section 4 discusses parameter tuning and robustness analysis methods. Extensive simulation results on a voice coil motor are provided in Section 5. Finally, concluding remarks are given in Section 6.

2. Ordinary GESO Design

Before introducing the design process of the RGESO, a brief description of the GESO is necessary to clearly illustrate the differences and advantages of the proposed method. Consider a general single-input–single-output plant which is uncertain and nonlinear as follows:

$$y^{\left(n\right)}=f\left(y^{\left(n-1\right)},\dots ,y,d\right)+b_{0}u$$

(1)

where y is the controlled output, u is the control input, d represents the external disturbance, and $b_0$ denotes the nominal value of the input gain. In the framework of ADRC, all the internal uncertainties and external disturbances are lumped into a total signal f. Here, f can be further divided into a low-frequency wide-band uncertainty $f_l$ and a high-frequency narrow-band disturbance $f_n$.

For the considered plant, it is assumed that the relative order and an approximate estimate of $b_0$ are known. They are prerequisites essential for the design of ADRC and RGESO-based control discussed in this paper. For a wide range of systems, these parameters can be derived either through dynamic modeling or by utilizing identification techniques based on experimental data. Additionally, it is assumed that the number of the narrow-band disturbances in $f_n$ is known, which determines the number of the resonant terms in the RGESO designed subsequently. Finally, to augment the high-order derivatives of f as states, it is assumed that f is differentiable and let $f^{\left(i\right)}$ represent its ith-order derivative.

To reject the disturbance with high-order dynamics, the original plant can be augmented with m additional states, including the total disturbance and its high-order derivatives. Then, (1) can be reformulated into a state space form as follows:

$$\left\{\begin{array}{cc}\hfill \dot{\mathit{X}}=& \mathit{A}\mathit{X}+{\mathit{B}}_{\mathbf{1}}u+{\mathit{B}}_{\mathbf{2}}{f}^{\left(m\right)}\hfill \\ \hfill y=& \mathit{C}\mathit{X}\hfill \end{array}\right.$$

(2)

where

$$\begin{array}{cc}\hfill \mathit{X}=& {\left[x_1{x}_2\cdots x_{n+1}\cdots x_{n+m}\right]}^T={[y\dot{y}\cdots f\cdots f^{(m-1)}]}^T,\hfill \end{array}$$

$\mathit{A}={\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& \ddots & \ddots & ⋮\\ ⋮& ⋮& \ddots & \ddots & 0\\ 0& 0& \cdots & 0& 1\\ 0& 0& \cdots & 0& 0\end{array}\right]}_{(n+m)\times (n+m)}$, ${\mathit{B}}_{\mathbf{1}}={\left[\begin{array}{ccccccc}0& \cdots & 0& b_0& 0& \cdots & 0\end{array}\right]}_{(n+m)\times 1}^T$, ${\mathit{B}}_{\mathbf{2}}={\left[\begin{array}{cccc}0& \cdots & 0& 1\end{array}\right]}_{(n+m)\times 1}^T$, and $\mathit{C}={\left[\begin{array}{cccc}1& 0& \cdots & 0\end{array}\right]}_{1\times (n+m)}$.

To estimate the total disturbance, a GESO [28], which is in the form of the conventional Luenberger observer, can be designed as

$$\dot{\widehat{\mathit{X}}}=\mathit{A}\widehat{\mathit{X}}+{\mathit{B}}_{\mathbf{1}}u+\mathit{L}e$$

(3)

where $\widehat{\mathit{X}}={\left[{\widehat{x}}_1{\widehat{x}}_2\cdots {\widehat{x}}_{n+1}\cdots {\widehat{x}}_{n+m}\right]}^T$ is the estimation of $\mathit{X}$, $e=y-{\widehat{x}}_1$ is the output estimation error, and $\mathit{L}={\left[l_1{l}_2\cdots l_{n+m}\right]}^T$ is the observer gain vector. To facilitate tuning, $\mathit{L}$ can be designed based on the parameterized bandwidth ${\omega }_o$ [16].

Combining (2) and (3), the dynamics of the estimation error can be obtained as

$$\dot{\tilde{\mathit{X}}}=\left(\mathit{A}-\mathit{LC}\right)\tilde{\mathit{X}}+{\mathit{B}}_{\mathbf{2}}{f}^{\left(m\right)}$$

(4)

where $\tilde{\mathit{X}}=\mathit{X}-\widehat{\mathit{X}}={\left[{\tilde{x}}_1{\tilde{x}}_2\cdots {\tilde{x}}_{n+1}\cdots {\widehat{x}}_{n+m}\right]}^T$. Taking the Laplace transformation of (4), the transfer functions from the total disturbance to the estimation errors can be obtained as

$$\begin{array}{c}{G}_f^{{\tilde{x}}_i}\left(s\right)=\frac{{\tilde{x}}_i\left(s\right)}{f\left(s\right)}=\frac{s^m(s^{i-1}+{\sum }_{j=1}^{i-1}{l}_j{s}^{i-1-j})}{s^{n+m}+{\sum }_{k=1}^{n+m}{l}_k{s}^{n+m-k}},i=1,2,\cdots ,n+1\hfill \end{array}$$

(5)

where $G_f^{{\tilde{x}}_i}\left(s\right)$ is a high-pass filter whose bandwidth is determined by ${\omega }_o$.

For the above GESO, ${\omega }_o$ should be much higher than the largest frequency of the disturbance. When dealing with disturbances that include high-frequency narrow-band components, the frequency response of the GESO is band-pass for all frequencies that are less than ${\omega }_o$. However, the high ${\omega }_o$ will degrade the robustness to the unmodeled dynamics and amplify the sensitivity to the measurement noise. In fact, such a high ${\omega }_o$ is unnecessary since only the band-pass around the frequency of the narrow-band disturbance is needed just as the resonant controller [6,32]. Considering this, a novel RGESO is proposed in this paper to effectively reject the high-frequency narrow-band disturbance with a moderate ${\omega }_o$.

3. RGESO-Based Controller Design

3.1. Frequency-Domain RGESO Design

The RGESO is designed by adding a resonant term $R\left(s\right)$ to the GESO. The transfer functions of the estimation errors are

$$\begin{array}{cc}\hfill & G_f^{{\tilde{x}}_i}=\frac{{\tilde{x}}_i\left(s\right)}{f\left(s\right)}=\frac{s^m(s^{i-1}+{\sum }_{j=1}^{i-1}{l}_j{s}^{i-1-j})}{s^{n+m}+{\sum }_{k=1}^{n+m}{l}_k{s}^{n+m-k}+R\left(s\right)},i=1,2,\cdots ,n+1\hfill \end{array}$$

(6)

where the resonant term is designed as

$$\begin{array}{c}\hfill R\left(s\right)=\frac{R_n\left(s\right)}{R_d\left(s\right)}=\prod _{i=1}^N{R}_i\left(s\right)=\prod _{i=1}^N\frac{s^2+{\zeta }_{in}s+{\omega }_i^2}{s^2+{\zeta }_{id}s+{\omega }_i^2}\end{array}$$

(7)

where N is the number of narrow-band disturbances, ${\omega }_i$ is the resonant frequency, ${\zeta }_{in}$ and ${\zeta }_{id}$ are the adjustable parameters which determine the depth and width of the peak, respectively, and they satisfy ${\zeta }_{in}\gg {\zeta }_{id}>0$. Such a form can improve the robustness to the finite word length of digital implementation and maintain a high resonant gain. The frequency response of $R_i\left(s\right)$ is

$$\left\{\begin{array}{cc}\hfill |R_i\left(jw\right)|=& {\zeta }_{in}/{\zeta }_{id},w={\omega }_i\hfill \\ \hfill |R_i\left(jw\right)|\approx & 1,w\gg {\omega }_i\mathrm{or}w\ll {\omega }_i\hfill \end{array}\right.$$

(8)

The maximal magnitude is $k_{ir}={\zeta }_{in}/{\zeta }_{id}$, which should be designed to be sufficiently large to effectively reject the high-frequency narrow-band disturbance. In addition, the width of the resonant peak can be adjusted by changing the value of ${\zeta }_{id}$. For a given $k_{ir}$, the larger ${\zeta }_{id}$ is, the wider the peak range is. Based on (6), we can have the transfer function from the total disturbance f to its estimate $\widehat{f}$ as

$$\begin{array}{c}\hfill G_f^{\widehat{f}}\left(s\right)\frac{{\sum }_{j=n+1}^{n+m}{l}_j{s}^{n+m-j}+R\left(s\right)}{s^{n+m}+{\sum }_{k=1}^{n+m}{l}_k{s}^{n+m-k}+R\left(s\right)}\end{array}$$

(9)

By employing (8), (9) can be approximated by

$$G_f^{\widehat{f}}\left(s\right)\approx \left\{\begin{array}{c}1,w={\omega }_i\hfill \\ \frac{{\sum }_{j=n+1}^{n+m}{l}_j{s}^{n+m-j}+1}{s^{n+m}+{\sum }_{k=1}^{n+m}{l}_k{s}^{n+m-k}+1},w\gg {\omega }_i\mathrm{or}w\ll {\omega }_i\hfill \end{array}\right.$$

(10)

According to (10), the addition of $R\left(s\right)$ can hardly affect the attenuation performance of the low-frequency wide-band uncertainties, and the perfect estimation of the narrow-band disturbances can be approximately achieved. Furthermore, we can use ${\omega }_o$ and ${\zeta }_{ik}(k=n,d)$ to represent the bandwidths for the low-frequency wide-band uncertainties and the high-frequency narrow-band disturbances, respectively. Such a manipulation can greatly facilitate the subsequent tuning process. To be specific, only the frequency band of the wide-band disturbance needs to be considered for the selection of ${\omega }_o$, while ${\zeta }_{ik}(k=n,d)$ can be chosen according to the specified requirements for the magnitude and width of the resonant peak. First, $k_{ir}$ should be large enough to guarantee unity gain and zero phase at the disturbance frequency. In addition, the width of the resonant peak should ensure the robustness to the estimation error of the disturbance frequency. According to the Bode sensitivity integral constraint [33], spillovers in the high frequencies can be caused by wide resonant peaks. Thus, the selection of the resonant peak width should make a trade-off between the robustness to the estimation error and the spillover effect.

3.2. Time-Domain RGESO Design

After the frequency-domain expression of the RGESO is designed, the next key step is the time-domain realization of (6) to obtain the estimations of the states and total disturbance. This process is outlined in the theorem that follows.

Theorem 1.

The state-space realization of (6) in the observable canonical form can be derived as

$$\left\{\begin{array}{cc}\hfill \dot{\mathit{Z}}=& {\mathit{A}}_{\mathit{z}}\mathit{Z}+{\mathit{B}}_{\mathbf{1}\mathit{z}}y+{\mathit{B}}_{\mathbf{2}\mathit{z}}u\hfill \\ \hfill \widehat{f}=& {\mathit{C}}_{\mathit{z}}\mathit{Z}\hfill \end{array}\right.$$

(11)

where

$$\begin{array}{c}{\mathit{A}}_{\mathit{z}}={\left[\begin{array}{ccccc}0& 0& \cdots & 0& -{\beta }_{n+m+2N}\\ 1& 0& \cdots & 0& ⋮\\ 0& \ddots & \ddots & ⋮& ⋮\\ ⋮& \ddots & \ddots & 0& -{\beta }_2\\ 0& \cdots & 0& 1& -{\beta }_1\end{array}\right]}_{(n+m+2N)\times (n+m+2N)},{\mathit{B}}_{\mathbf{1}\mathit{z}}={\left[\begin{array}{c}0\\ ⋮\\ 0\\ {\alpha }_{m+2N-1}\\ ⋮\\ {\alpha }_1\\ {\alpha }_0\end{array}\right]}_{(n+m+2N)\times 1},\hfill \\ {\mathit{B}}_{\mathbf{2}\mathit{z}}=-b_0{\left[\begin{array}{c}{\alpha }_{m+2N-1}\\ ⋮\\ {\alpha }_1\\ {\alpha }_0\\ 0\\ ⋮\\ 0\end{array}\right]}_{(n+m+2N)\times 1},and{\mathit{C}}_{\mathit{z}}={\left[\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right]}_{(n+m+2N)\times 1}.\hfill \end{array}$$

By using the disturbance estimation, the state estimations can be obtained as

$${\dot{\widehat{\mathit{X}}}}_{\mathit{M}}={\mathit{A}}_{\mathit{M}}{\widehat{\mathit{X}}}_{\mathit{M}}+{\mathit{B}}_{\mathbf{1}\mathit{M}}u+{\mathit{B}}_{\mathbf{2}\mathit{M}}\widehat{f}+{\mathit{L}}_{\mathit{M}}e$$

(12)

where ${\widehat{\mathit{X}}}_{\mathit{M}}={\left[{\widehat{x}}_1{\widehat{x}}_2\cdots {\widehat{x}}_n\right]}^T$ is the vector of state estimations,

$${\mathit{A}}_{\mathit{M}}={\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& \ddots & \ddots & ⋮\\ ⋮& ⋮& \ddots & \ddots & 0\\ 0& 0& \cdots & 0& 1\\ 0& 0& \cdots & 0& 0\end{array}\right]}_{n\times n},{\mathit{B}}_{\mathbf{1}\mathit{M}}={\left[\begin{array}{c}0\\ ⋮\\ 0\\ b_0\end{array}\right]}_{n\times 1},{\mathit{B}}_{\mathbf{2}\mathit{M}}={\left[\begin{array}{c}0\\ ⋮\\ 0\\ 1\end{array}\right]}_{n\times 1},and{\mathit{L}}_{\mathit{M}}={\left[\begin{array}{c}{l}_1\\ l_2\\ ⋮\\ l_n\end{array}\right]}_{n\times 1}.$$

Proof. 

According to (1), the total disturbance can be expressed by

$$f\left(s\right)=y\left(s\right)s^n-b_{0}u\left(s\right)$$

(13)

Substituting (13) into (9) can obtain the frequency-domain expression of the disturbance estimation as

$$\begin{array}{c}\hfill \widehat{f}\left(s\right)=G_y^{n+1}\left(s\right)y\left(s\right)+G_u^{n+1}\left(s\right)u\left(s\right)\end{array}$$

(14)

where $G_y^{n+1}\left(s\right)$ and $G_u^{n+1}\left(s\right)$ are in the following forms:

$$\begin{array}{cc}\hfill G_y^{n+1}\left(s\right)=& \frac{s^n({\sum }_{j=n+1}^{n+m}{l}_j{s}^{n+m-j}+R\left(s\right))}{s^{n+m}+{\sum }_{k=1}^{n+m}{l}_k{s}^{n+m-k}+R\left(s\right)},G_u^{n+1}\left(s\right)=-\frac{b_0({\sum }_{j=n+1}^{n+m}{l}_j{s}^{n+m-j}+R\left(s\right))}{s^{n+m}+{\sum }_{k=1}^{n+m}{l}_k{s}^{n+m-k}+R\left(s\right)}\hfill \end{array}$$

(15)

which are both proper. It can be seen that the disturbance estimation can be calculated by using the output measurement and the control signal. Then, (15) can be further rewritten as

$$\begin{array}{cc}\hfill \widehat{f}\left(s\right)=& \frac{{\sum }_{j=0}^{m+2N-1}{s}^{m+2N-1-j}({\alpha }_j{s}^{n}y\left(s\right)-b_0{\alpha }_{j}u\left(s\right))}{s^{n+m+2N}+{\sum }_{k=1}^{n+m+2N}{\beta }_k{s}^{n+m+2N-k}}\hfill \end{array}$$

(16)

where both ${\alpha }_j$ and ${\beta }_k$ are the functions of ${\omega }_o$, ${\omega }_i$, ${\zeta }_{in}$, and ${\zeta }_{id}$.

Next, we validate the correctness of (12). Based on (12) and the first n differential equations in (2), the transfer functions from the disturbance estimation error to the state estimation errors can be obtained as

$$\begin{array}{cc}\hfill G_{{\tilde{x}}_{n+1}}^{{\tilde{x}}_i}\left(s\right)=& \frac{G_f^{{\tilde{x}}_i}\left(s\right)}{G_f^{{\tilde{x}}_{n+1}}\left(s\right)}=\frac{s^{i-1}+{\sum }_{j=1}^{i-1}{l}_j{s}^{i-1-j}}{s^n+{\sum }_{j=1}^n{l}_j{s}^{n-j}},i=1,2,\cdots ,n\hfill \end{array}$$

(17)

Combining (17) with $G_f^{{\tilde{x}}_{n+1}}\left(s\right)$, the transfer functions from the total disturbance to the state estimation errors can be obtained as

$$\begin{array}{cc}\hfill G_f^{{\tilde{x}}_i}\left(s\right)=& G_f^{{\tilde{x}}_{n+1}}\left(s\right)G_{{\tilde{x}}_{n+1}}^{{\tilde{x}}_i}\left(s\right)=\frac{s^m(s^{i-1}+{\sum }_{j=1}^{i-1}{l}_j{s}^{i-1-j})}{s^{n+m}+{\sum }_{k=1}^{n+m}{l}_k{s}^{n+m-k}+R\left(s\right)},i=1,2,\cdots ,n\hfill \end{array}$$

(18)

which is exactly the form that is required in (6). □

3.3. Control Law Design

After deriving the disturbance and state estimations, the control law can be designed as

$$\begin{array}{cc}\hfill u=& [k_1\left(r-{\widehat{x}}_1\right)+k_2\left(\dot{r}-{\widehat{x}}_2\right)+\cdots +k_n\left(r^{\left(n-1\right)}-{\widehat{x}}_n\right)-{\widehat{x}}_{n+1}]/b_0\hfill \end{array}$$

(19)

where $k_i$ $(i=1,2,\cdots ,n)$ is the control gain. Substituting the control law into (1), the closed-loop system can be derived as

$$\begin{array}{cc}\hfill y^{\left(n\right)}=& k_1\left(r-{\widehat{x}}_1\right)+k_2\left(\dot{r}-{\widehat{x}}_2\right)+\cdots +k_n\left(r^{\left(n-1\right)}-{\widehat{x}}_n\right)+\left(f-{\widehat{x}}_{n+1}\right)\hfill \end{array}$$

(20)

Although the RGESO-based control method is linear, it can be applied to nonlinear, time-varying, and uncertain processes with little model information. Since the total disturbance is treated as one signal, there is no difference between the linear and nonlinear dynamics from a signal point of view. In the ordinary ESO design, only the order and the input gain of the plant need to be known in advance. For the RGESO, the approximate frequency range of the disturbance is an extra requirement for parameter tuning.

4. ANF for Disturbance Frequency Estimation

In this section, the ANF is designed to provide real-time estimation of the instantaneous frequency (IF) of the narrow-band disturbance.The definition of the IF for a real sinusoidal signal is given as

$${\omega }_{IF}=\frac{d{p}_{\varphi }\left(t\right)}{dt}$$

(21)

where $p_{\varphi }\left(t\right)$ is the phase of the signal. The ideal frequency response of an ANF satisfies

$$E\left(e^{jw}\right)=\left\{\begin{array}{cc}0,\hfill & w={\omega }_i\hfill \\ 1,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(22)

There are many forms of finite or infinite impulse response approximations of (22). The minimal parameter ANF with zeros distributed on the unit circle [34] is the most popular one, which can be expressed as

$$H\left(z^{-1}\right)=\frac{F\left(z^{-1}\right)}{F\left(\gamma z^{-1}\right)}=\prod _{i=1}^N\frac{1-2cos{\theta }_i{z}^{-1}+z^{-2}}{1-2cos{\theta }_i\gamma z^{-1}+{\gamma }^2{z}^{-2}}$$

(23)

where $z^{-1}$ is the unit delay operator, N is the number of sinusoidal disturbances, ${\theta }_i={\omega }_i\Delta T$ (${\omega }_i$ and $\Delta T$ are the notch frequency and the sampling time, respectively), and $\gamma$ ($0\le \gamma <1$) is the design parameter that determines the notch width and depth. When $\gamma <1$, $H\left(z^{-1}\right)$ is a stable filter. As $\gamma$ decreases, the notch width widens, resulting in a faster convergence speed but lower estimation accuracy. Due to the constraint that the zeros are on the unit circle, the coefficients of $F\left(z^{-1}\right)$ have a mirror symmetric form which can be written as

$$\begin{array}{c}\hfill F\left(z^{-1}\right)=\sum _{i=0}^{N-1}{a}_i(z^{-i}+z^{-2N+i})+a_N{z}^{-N}\end{array}$$

(24)

where $a_0=1$. When the ANF input is $v\left(t\right)$, the output prediction $e\left(t\right)$ can be expressed as

$$\begin{array}{cc}\hfill e\left(t\right)=& v\left(t\right)+v(t-2N)-{\gamma }^{2N}e(t-2N)-a_N(t-1)\left[{\gamma }^{N}e(t-N)-v(t-N)\right]-\hfill \\ & \sum _{i=1}^{N-1}{a}_i(t-1)\left[{\gamma }^{i}e(t-i)+{\gamma }^{2N-i}e(t-2N+i)-v(t-i)-v(t-2N+i)\right]\hfill \\ \hfill =& v\left(t\right)+v(t-2N)-{\gamma }^{2N}e(t-2N)-{\mathit{\varphi }}^T\left(t\right)\mathit{a}(t-1)\hfill \end{array}$$

(25)

where

$$\begin{array}{ccc}\hfill & \mathit{a}={\left[a_1{a}_2\cdots a_N\right]}^T,\hfill & \hfill \mathit{\varphi }={\left[{\varphi }_1{\varphi }_2\cdots {\varphi }_N\right]}^T\end{array}$$

(26)

and

$${\varphi }_i\left(t\right)=\left\{\begin{array}{c}{\gamma }^{i}e(t-i)+{\gamma }^{2N-i}e(t-2N+i)-v(t-i)-v(t-2N+i),1\le i\le N-1\hfill \\ {\gamma }^{N}e(t-N)-v(t-N),i=N\hfill \end{array}\right.$$

(27)

Each component of $\mathit{a}$ is a polynomial of $cos{\theta }_i$ $(i=1,2,\cdots ,N)$.

In the traditional ANF design, $\mathit{a}$ is a directly adapted parameter vector, and the complicated factorization is subsequently required to obtain the frequency estimation. However, in this paper, the parameter vector $\mathit{\theta }={\left[{\theta }_1{\theta }_2\cdots {\theta }_N\right]}^T$ is adapted directly. The adaptation of $\mathit{\theta }$ aims to minimize the following weighed quadratic prediction-error criterion:

$$V_t\left(\mathit{\theta }\right)=\frac{1}{2}\sum _{\tau =1}^t{\lambda }^{t-\tau }e{(\tau ,\mathit{\theta })}^2$$

(28)

where $\lambda$ ($0<\lambda \le 1$) is the forgetting factor. A small value of $\lambda$ implies that only recent data are included in the criterion. This allows a fast-varying $\mathit{\theta }$ to be tracked properly, which results in a low bias error in the estimate. The drawback is that the estimate will be more influenced by the noise, which gives a high variance error in the estimate. Therefore, the selection of $\lambda$ is a tradeoff between how fast the filter should track the changes in and the energy of the noise. From (28), we can see that the identification priority is given to the frequency component with larger power. Since there are both narrow-band and wide-band signals in the total disturbance, a band-pass filter is added before the disturbance estimation enters the ANF to eliminate the effects of the wide-band components. The identification of $\mathit{\theta }$ is a highly nonlinear least square problem. To develop the LM algorithm, the error gradient should be determined first. The error gradient with respect to the parameter vector is

$$\begin{array}{c}\hfill {{\mathit{\psi }}_{\mathit{\theta }}}^T\left(t\right)=-\frac{\partial e\left(t\right)}{\partial {\mathit{\theta }}^T}=-\frac{\partial e\left(t\right)}{\partial {\mathit{a}}^T}\frac{\partial \mathit{a}}{\partial {\mathit{\theta }}^T}=-\left[\frac{\partial e\left(t\right)}{\partial {\theta }_1}\frac{\partial e\left(t\right)}{\partial {\theta }_2}\cdots \frac{\partial e\left(t\right)}{\partial {\theta }_N}\right]\end{array}$$

(29)

Differentiating both sides of (25) can obtain

$${{\mathit{\psi }}_{\mathit{a}}}^T\left(t\right)=-\frac{\partial e\left(t\right)}{\partial {\mathit{a}}^T}=\frac{{\mathit{\varphi }}^T\left(t\right)}{F\left(\gamma z^{-1}\right)}$$

(30)

This result implies that the gradient vector is effectively the regression ${\mathit{\varphi }}^T\left(t\right)$ filtered by $F\left(\gamma z^{-1}\right)$. Such filtering can improve the convergence of the algorithm. To reduce the computation complexity, the following filtered variables:

$$v_F\left(t\right)=\frac{v\left(t\right)}{F\left(\gamma z^{-1}\right)},e_F\left(t\right)=\frac{e\left(t\right)}{F\left(\gamma z^{-1}\right)}$$

(31)

are introduced to calculate ${{\mathit{\psi }}_{\mathit{a}}}^T\left(t\right)$. The term $\partial \mathit{a}/\partial {\mathit{\theta }}^T$ in (29) is a Jacobian matrix:

$$\frac{\partial \mathit{a}}{\partial {\mathit{\theta }}^T}=\left[\begin{array}{ccc}{\displaystyle \frac{\partial a_1}{\partial {\theta }_1}}& \cdots & {\displaystyle \frac{\partial a_1}{\partial {\theta }_N}}\\ ⋮& \ddots & ⋮\\ {\displaystyle \frac{\partial a_N}{\partial {\theta }_1}}& \cdots & {\displaystyle \frac{\partial a_N}{\partial {\theta }_N}}\end{array}\right]$$

(32)

After some mathematical manipulations, the recursive expression can be obtained as

$$\frac{\partial a_i}{\partial {\theta }_j}=2cos{\theta }_j\frac{\partial a_{i-1}}{\partial {\theta }_j}-\frac{\partial a_{i-2}}{\partial {\theta }_j}+2a_{i-1}sin{\theta }_j$$

(33)

for $2\le i\le N$ and $1\le j\le N$. The first two terms can be easily derived as

$$\frac{\partial a_0}{\partial {\theta }_j}=0,\frac{\partial a_1}{\partial {\theta }_j}=2sin{\theta }_j$$

(34)

Thus, the partial derivatives can be calculated efficiently based on the recursion. To improve the convergence speed, the posteriori error $\overline{e}\left(t\right)$ is used for the error prediction and the gradient calculation in (25) and (30), since it is expected to yield a better estimation than $e\left(t\right)$. $\overline{e}\left(t\right)$ is calculated using the latest parameters $\mathit{a}\left(t\right)$ as

$$\begin{array}{c}\hfill \overline{e}\left(t\right)=v\left(t\right)+v(t-2N)-{\gamma }^{2N}\overline{e}(t-2N)-{\overline{\mathit{\varphi }}}^T\left(t\right)\mathit{a}\left(t\right)\end{array}$$

(35)

where

$$\overline{\mathit{\varphi }}={\left[{\overline{\varphi }}_1{\overline{\varphi }}_2\cdots {\overline{\varphi }}_N\right]}^T$$

(36)

and

$$\overline{{\varphi }_i}\left(t\right)=\left\{\begin{array}{c}{\gamma }^i\overline{e}(t-i)+{\gamma }^{2N-i}\overline{e}(t-2N+i)-v(t-i)-v(t-2N+i),1\le i\le N-1\hfill \\ {\gamma }^N\overline{e}(t-N)-v(t-N),i=N\hfill \end{array}\right.$$

(37)

After obtaining the error gradient, the next step is to design the recursive LM (RLM) algorithm to estimate the frequencies. The RLM method can be summarized as

$$\left\{\begin{array}{cc}\hfill \mathit{R}\left(t\right)=& \lambda \mathit{R}(t-1)+\overline{\lambda }[{\mathit{\psi }}_{\mathit{\theta }}\left(t\right){\mathit{\psi }}_{\mathit{\theta }}{\left(t\right)}^T+\delta \mathit{I}]\hfill \\ \hfill \mathit{\theta }\left(t\right)=& \mathit{\theta }(t-1)+\overline{\lambda }{\mathit{R}}^{-1}\left(t\right){\mathit{\psi }}_{\mathit{\theta }}\left(t\right)e\left(t\right)\hfill \end{array}\right.$$

(38)

where $\overline{\lambda }=1-\lambda$, $\mathit{R}\left(t\right)=\overline{\lambda }\mathit{H}\left(\mathit{t}\right)+\delta \mathit{I}$, $\mathit{H}\left(\mathit{t}\right)$ is the approximation of the Hessian matrix, and $\mathit{I}$ is a identity matrix. When $\delta \to \infty$, the diagonal elements of $\mathit{R}\left(t\right)$ will be dominant, and thus, the parameters will adapt in the steepest descent direction with a short step length. The direction is acceptable, but the convergence speed may be very low. When $\delta =0$, the algorithm just becomes the SGN method, which may not achieve the convergence. The LM method can make a trade-off between the acceptance of the direction and the convergence speed through adding the term $\delta \mathit{I}$. To avoid the matrix inversion in (38), the approximation method is used as follows:

$$\left\{\begin{array}{cc}\hfill \mathit{S}\left(t\right)=& {{\mathit{\psi }}_{\mathit{\theta }}}^{*}{\left(t\right)}^T\mathit{P}(t-1){{\mathit{\psi }}_{\mathit{\theta }}}^{*}\left(t\right)+\lambda \Lambda (t-1)\hfill \\ \hfill \mathit{P}\left(t\right)=& \frac{1}{\lambda }[\mathit{P}(t-1)-\mathit{P}(t-1){{\mathit{\psi }}_{\mathit{\theta }}}^{*}\left(t\right){\mathit{S}}^{-1}\left(t\right){{\mathit{\psi }}_{\mathit{\theta }}}^{*}{\left(t\right)}^T\mathit{P}(t-1)]\hfill \\ \hfill \mathit{\theta }\left(t\right)=& \mathit{\theta }(t-1)+\mathit{P}\left(t\right){\mathit{\psi }}_{\mathit{\theta }}\left(t\right)e\left(t\right)\hfill \end{array}\right.$$

(39)

where ${{\mathit{\psi }}_{\mathit{\theta }}}^{*}$ is a $d_m\times 2$ matrix ($d_m$ is the dimension of the parameter vector) with the first column being ${\mathit{\psi }}_{\mathit{\theta }}$ and the second column being a $d_m\times 1$ zero vector except for its ${[\left(t\mathrm{mod}{d}_m\right)+1]}^{th}$ element being one, which can be represented as

$${{\mathit{\psi }}_{\mathit{\theta }}}^{*}\left(t\right)={\left[\begin{array}{ccccc}& & {\mathit{\psi }}_{\mathit{\theta }}{\left(t\right)}^T& & \\ 0& \cdots & 1& \cdots & 0\end{array}\right]}^T$$

(40)

and $\Lambda (t-1)$ satisfies

$${\Lambda }^{-1}(t-1)=\left[\begin{array}{cc}1& 0\\ 0& d_m\delta (t-1)\end{array}\right]$$

(41)

$\delta \left(t\right)$ should be adapted on-line such that the cost function value in (28) decreases using the RLM method. The predicted reduction of (28) at time t is

$$r_p\left(t\right)={\Xi }^T\mathit{P}\left(t\right)\Xi$$

(42)

where $\Xi ={\mathit{\psi }}_{\mathit{\theta }}\left(t\right)e\left(t\right)$. The actual reduction $r_a\left(t\right)$ is defined as

$$r_a\left(t\right)=\frac{1}{2d_m}\left[\sum _{\tau =t-2d_m+1}^{t-d_m}{\overline{e}}^2\left(\tau \right)-\sum _{\tau =t-d_m+1}^t{\overline{e}}^2\left(\tau \right)\right]$$

(43)

Then, the value of $\delta \left(t\right)$ can be adjusted by the following rules:

$$\delta \left(t\right)=\left\{\begin{array}{c}\kappa \delta (t-1),\mathrm{if}{r}_a\left(t\right)<ϵr_p\left(t\right)\hfill \\ \delta (t-1)/\kappa ,\mathrm{if}{r}_a\left(t\right)>(1-ϵ)r_p\left(t\right)\hfill \\ \delta (t-1),\mathrm{otherwise}\hfill \end{array}\right.$$

(44)

where $ϵ$ ($0<ϵ<0.5$) and $\kappa$ ($\kappa >1$) determine the time and magnitude of the adjustment, respectively.

At the start of the data processing, it is advisable to apply the algorithm with wider notches (i.e., a smaller value of $\gamma$), thus increasing filter sensitivity to the presence of input sine waves. After convergence, it is recommended that a larger $\gamma$ is used, which improves the asymptotic performance. Since the estimates of the variables in the gradient are initially quite poor, they should therefore be assigned a lower weight in the criterion compared to later measurements. Therefore, transient processes are designed for $\gamma$ and $\lambda$ to improve the convergence speed at the beginning as

$$\begin{array}{cc}\hfill \lambda (t+1)=& {\lambda }_r\lambda \left(t\right)+(1-{\lambda }_r)\lambda (\infty ),\gamma (t+1)={\gamma }_r\gamma \left(t\right)+(1-{\gamma }_r)\gamma (\infty )\hfill \end{array}$$

(45)

where $\lambda (\infty )$ and $\gamma (\infty )$ are the desired steady values, and ${\lambda }_r$ and ${\gamma }_r$ are the time constants of the processes.

Although the RLM is designed to improve the convergence performance, convergence cannot always be ensured since both the signal and the algorithm are stochastic in the presence of the measurement noise. Once the estimation converges to a false value well beyond our design consideration, the closed-loop performance will be degraded, and the stability may even be violated. Considering that the reliability is the most important aspect in practice, monitoring the estimation results based on the prior knowledge about the frequency range is necessary. Thus, the following strategy is proposed to further enhance the convergence performance:

$${\theta }_i=\left\{\begin{array}{c}{\theta }_i^{+},{\theta }_i>{\theta }_i^{+}\hfill \\ {\theta }_i^{-},{\theta }_i<{\theta }_i^{-}\hfill \\ {\theta }_i,\mathrm{otherwise}\hfill \end{array}\right.$$

(46)

where ${\theta }_i^{+}$ and ${\theta }_i^{-}$ are the upper and lower bounds for ${\theta }_i$, respectively.

Until now, the complete form of the proposed control scheme has been presented.

5. Parameter Tuning and Stability Robustness Analysis

In this section, parameter tuning and stability robustness analysis are presented for the proposed method.

5.1. Parameter Tuning

The DOBC has a two-degree-of-freedom (2-DOF) scheme in nature [35,36], under which the tracking and disturbance rejection performance can be designed separately. By taking advantage of this property, a two-step tuning process is used in this paper.

First, the error feedback gains in (19) are determined based on the desired tracking performance including the settling time and the overshoot. For a second-order plant, the error feedback gains $k_1$ and $k_2$ can be selected based on the desired closed-loop transfer function:

$$G_{cl}=\frac{k_1}{s^2+k_{2}s+k_1}$$

(47)

where $k_1={\omega }_c^2$, $k_2=2\varsigma {\omega }_c$, and ${\omega }_c$ and $\varsigma$ are the natural frequency and the damping ratio, respectively. Then, the RGESO parameters are derived by using the optimization method. Based on (6) and (13), the state and disturbance estimations can be expressed as

$$\begin{array}{cc}\hfill {\widehat{x}}_i\left(s\right)=& y\left(s\right)s^{(i-1)}-G_f^{{\tilde{x}}_i}\left(s\right)[y\left(s\right)s^n-b_{0}u\left(s\right)]=G_y^i\left(s\right)y\left(s\right)+G_u^i\left(s\right)u\left(s\right),i=1,2,\cdots ,n\hfill \end{array}$$

(48)

where $G_y^i\left(s\right)=s^{(i-1)}-G_f^{{\tilde{x}}_i}\left(s\right)s^n$ and $G_u^i\left(s\right)=b_0{G}_f^{{\tilde{x}}_i}\left(s\right)$. Substituting (48) into the control signal in (19) has

$$\begin{array}{c}\hfill u\left(s\right)=C_R\left(s\right)r\left(s\right)-C_F\left(s\right)y\left(s\right)\end{array}$$

(49)

By defining

$$\begin{array}{c}\hfill C_R\left(s\right)=\frac{{\sum }_{i=1}^n{k}_i{s}^{(i-1)}}{b_0+{\sum }_{i=1}^n{k}_i{G}_u^i\left(s\right)+G_u^{n+1}\left(s\right)},C_F\left(s\right)=\frac{{\sum }_{i=1}^n{k}_i{G}_y^i\left(s\right)y\left(s\right)+G_y^{n+1}\left(s\right)}{b_0+{\sum }_{i=1}^n{k}_i{G}_u^i\left(s\right)+G_u^{n+1}\left(s\right)}\end{array}$$

(50)

the 2-DOF control structure can be represented by the block diagram in Figure 1, where $G_p\left(s\right)$ denotes the plant dynamics. Therefore, the closed-loop transfer function from the disturbance d to the output y can be derived as

$$G_d^y\left(s\right)=\frac{G_p\left(s\right)}{1+C_F\left(s\right)G_p\left(s\right)}$$

(51)

Then, the RGESO parameters can be obtained by solving the following constrained optimization problem:

$$\left\{\begin{array}{c}\underset{{\omega }_o,{\zeta }_{in,}{\zeta }_{id}}{min}{\omega }_o\hfill \\ \mathrm{subject}\mathrm{to}\hfill \\ {\omega }_o>{\omega }_l,{\zeta }_{in}>0,{\zeta }_{id}>0,i=1,2,\cdots ,N\hfill \\ |G_d^y\left(j{\omega }_i\right)|<S_i\mathrm{dB}\hfill \\ |G_d^y\left(j{\omega }_i(1\pm {\gamma }_i)\right)|<P_i\mathrm{dB}\hfill \\ {∥G_d^y\left(s\right)∥}_{\infty }<M\mathrm{dB}\hfill \\ \mathrm{Re}\left(s_k\right)<0,k=1,2,\cdots ,R\hfill \end{array}\right.$$

(52)

To be specific, the optimization objective is to minimize ${\omega }_o$ and the constraint ${\omega }_o>{\omega }_l$ mean to reduce the value of ${\omega }_o$ as much as possible while ensuring the attenuation performance of the low-frequency wide-band uncertainty. A smaller ${\omega }_o$ can enhance the robustness to the unmodeled dynamics (such as the time delay) [37,38], decrease the sensitivity to the measurement noise, and reduce the requirement for the sampling rate [39]. For a DOBC system, the observer bandwidth is generally much larger than that of the closed-loop bandwidth, and it was suggested that the bandwidth of the ESO is three to five times the control bandwidth [16]. Therefore, ${\omega }_l=3{\omega }_c$ is used in this paper. $S_i$ can be determined by the expected magnitude of the steady error caused by the narrow-band disturbance. To guarantee the robustness to the estimation errors of the ANF, the bandwidths of the resonant peaks are adjusted by the parameters ${\gamma }_i$ and $P_i$. To avoid the spillovers at other frequencies, the maximal singular value of $G_d^y$ is constrained to be less than M dB. Since a higher resonant peak is required to achieve the desired attenuation magnitude $S_i$ for a higher-frequency disturbance, the desired attenuation magnitudes can generally be guaranteed across all possible disturbance frequencies if they can be met for the maximal disturbance frequencies. Thus, the optimization is performed with the disturbance frequencies being their possible maxima. The last constraint is to guarantee that the closed-loop system is stable, and the real part of the pole $s_k$ is restricted to be negative. Many mature methods such as the interior point method, sequential quadratic programming, and intelligent algorithms can be used to solve the above optimization problem.

5.2. Stability Robustness Analysis

The above tuning process is performed with the disturbance frequencies being their possible maxima. However, the disturbance frequencies are time-varying and the RGESO is adapted based on the frequency estimations provided by the ANF. Thus, the closed-loop stability for all the possible disturbance frequencies should be verified. Define $\mathit{\omega }=[{\omega }_1{\omega }_2\cdots {\omega }_N]$, and suppose ${\omega }_i\in \left[{\omega }_i^{-},{\omega }_i^{+}\right]$ $(i=1,2,\cdots ,N)$ where ${\omega }_i^{-}$ and ${\omega }_i^{+}$ are the lower and upper bounds for ${\omega }_i$, respectively. The set $\mathit{Q}$ is defined as

$$\mathit{Q}=\left\{\mathit{\omega }\in {\mathit{R}}^N:{\omega }_i^{-}\le {\omega }_i\le {\omega }_i^{+}\mathrm{for}i=1,2,\cdots ,N\right\}$$

(53)

The closed-loop characteristic polynomial can be obtained as

$$\begin{array}{cc}\hfill G_c(s,\mathit{\omega })=& b_0+\sum _{i=1}^n{k}_i(G_u^i\left(s\right)+G_p\left(s\right)G_y^i\left(s\right))+G_u^{n+1}\left(s\right)+G_p\left(s\right)G_y^{n+1}\left(s\right)\hfill \end{array}$$

(54)

Define a family of polynomials ${\mathit{P}}_{\mathit{f}}=\left\{G_c(s,\mathit{\omega }):\mathit{\omega }\in \mathit{Q}\right\}$. The stability of ${\mathit{P}}_{\mathit{f}}$ for $\mathit{\omega }\in \mathit{Q}$ can be verified by calculating the eigenvalues of the polynomials for all frequency combinations, whereas the computations are rather time-consuming. Since the polynomial coefficients in (54) are not independent with each other, the traditional Kharitonov theorem cannot be applied. We find that the polynomial coefficients are all multilinear functions of ${\omega }_i^2$. Based on this characteristic, an effective and intuitive graphical method that combines the zero exclusion condition [40] and the mapping theorem [41] is proposed in this paper. To facilitate understanding, these two theorems are presented below:

Theorem 2

([40]). For a family of polynomials ${\mathit{P}}_{\mathit{f}}=\left\{G_c(s,\mathit{\omega }):\mathit{\omega }\in \mathit{Q}\right\}$, the value set at frequency $w\in \mathit{R}$ is defined as

$$G_c(jw,\mathit{Q})=\left\{G_c(jw,\mathit{\omega }):\mathit{\omega }\in \mathit{Q}\right\}$$

(55)

If the polynomials in the family have an invariant degree, and the polynomial coefficients are all continuous functions of $\mathit{\omega }$, and furthermore, there exists at least one stable member $G_c(s,{\mathit{\omega }}^0)$ in the family, then the family ${\mathit{P}}_{\mathit{f}}$ is robustly stable if the complex plane origin is excluded from the value set $G_c(jw,\mathit{Q})$ at all the frequencies $w\ge 0$, i.e., ${\mathit{P}}_{\mathit{f}}$ is robustly stable if

$$0\notin G_c(jw,\mathit{Q}),\forall w\ge 0$$

(56)

Theorem 3

([41]). Suppose $G_c(jw,\mathit{\omega })$ is a multi-affine function of $\mathit{\omega }$ ($\mathit{\omega }\in \mathit{Q}$). Let ${\mathit{Q}}_{\mathit{v}}$ denote the corners of $\mathit{Q}$, and conv${\mathit{A}}_{\mathit{v}}$ denote the convex hull of a set ${\mathit{A}}_{\mathit{v}}\subset {\mathit{R}}^2$. Then,

$$\mathrm{conv}{G}_c(jw,\mathit{Q})=\mathrm{conv}{G}_c(jw,{\mathit{Q}}_{\mathit{v}})$$

(57)

According to Theorems 1 and 2, at each frequency $w\ge 0$, we only need to calculate the values of $G_c(jw,\mathit{Q})$ on the corner points. Then, the convex hull can just be found by using the Graham scanning method, which can enhance the computation efficiency greatly. Finally, the stability robustness can just be verified by plotting all the convex hulls for each sampled frequency $w_s\ge 0$.

6. Numerical Simulations and Discussion

In this section, numerical simulations are performed on the voice coil motor (VCM) model [42] to find whether the proposed method can achieve the high-precision position control under the multi-band disturbances. The nominal model of the VCM is

$$P\left(s\right)=\frac{Y\left(s\right)}{U\left(s\right)}=\frac{b_0}{s^2+p_{1}s+p_2}$$

(58)

where $b_0=2.239\times {10}^4$, $p_1=131.5$, and $p_2=84.9$. The external disturbance is added in the input channel as

$$d\left(t\right)=\sum _{i=1}^3{m}_{i}sin\left({\omega }_{i}t\right)+n_g\left(t\right)$$

(59)

where $m_i$ and ${\omega }_i$ are the magnitude and frequency of each sinusoidal signal, respectively, and $n_g\left(t\right)$ represents a white Gaussian noise. In the following simulations, let $m_i=1$, ${\omega }_i=200i$ rad/s, and set the standard derivation of $n_g\left(t\right)$ as 0.3544. According to the control objective, the parameters in (52) are selected as

$$\begin{array}{c}\hfill {\omega }_l=300,S_1=S_2=S_3=-100,P_1=P_2=P_3=-80,{\gamma }_1={\gamma }_2={\gamma }_3=0.02,M=-40,\end{array}$$

where ${\omega }_l=3{\omega }_c$, and ${\omega }_c=100$ is chosen to make the regulation time be less than 0.03 s. Additionally, the closed-loop system is required to have stability robustness against the variations of narrow-band frequencies within $[-0.1,0.1]$. To make it clearer, the block diagram of the overall control structure for the VCM, including the plant, the controller, the observer, and the adding disturbances, is presented in Figure 2. MATLAB R2023a is used for conducting the simulations, and the simulations are performed on a computer equipped with an Intel i7 processor and memory of 16 GB.

6.1. Performance Test for the ANF

In this subsection, comparisons between the proposed ANF and the SGN-based ANF are performed to demonstrate the advantages of the RLM algorithm in terms of convergence. The periodic signal includes three sinusoidal components with $m_i=1$ and ${\omega }_i=200i$ rad/s. The parameters for the two algorithms are presented in

Table 1. Parameters of two types of ANFs.

RLM Method	SGN Method
$\lambda \left(0\right)=0.7$, ${\lambda }_r=0.99$, $\lambda (\infty )=0.992$,	$\lambda \left(0\right)=0.7$, ${\lambda }_r=0.99$, $\lambda (\infty )=0.992$
$\gamma \left(0\right)=0.8$, ${\gamma }_r=0.99$, $\gamma (\infty )=0.95$,	$\gamma \left(0\right)=0.8$, ${\gamma }_r=0.99$, $\gamma (\infty )=0.95$,
$\mathit{P}\left(0\right)=1000\mathit{I}$, $\mathit{\theta }\left(0\right)=0.01{\mathit{\theta }}^{*}$,	$\mathit{P}\left(0\right)=1000\mathit{I}$, $\mathit{\theta }\left(0\right)=0.01{\mathit{\theta }}^{*}$
$\delta \left(0\right)=1.0\times {10}^{-4}$, $ϵ=0.25$, $\kappa =2$

, where ${\mathit{\theta }}^{*}$ represents the vector of true values. To ensure a fair comparison, the common parameters have been selected identically for both algorithms. A large $\mathit{P}\left(0\right)$ can make $\mathit{\theta }\left(t\right)$ have a large update initially to reflect the little confidence in $\mathit{\theta }\left(0\right)$. Therefore, we let $\mathit{P}\left(0\right)=1000\mathit{I}$. Here, we let $\mathit{\theta }\left(0\right)=0.01{\mathit{\theta }}^{*}$ to make the simulation scenario more challenging. The values of the remaining parameters are determined based on the value ranges of each parameter, the parameter tuning criteria, and extensive simulation results. The number of data samples is 10,000, and the convergence is defined as when the final estimation errors of all three frequencies are less than 5%. Without the online supervision strategy, 100 simulations are conducted for different SNRs and the numbers of convergence for the two algorithms are shown in

Table 2. Number of convergence for 100 simulations under different SNRs.

SNR	RLM Method	SGN Method
3 dB	69 (100)	52 (100)
0 dB	63 (100)	46 (100)
−10 dB	39 (100)	31 (100)

. It is clear that the convergence cannot be achieved for all the scenarios, whereas the RLM-based ANF is more advantageous especially when the SNR is relatively high. Furthermore, it is found that the convergence is guaranteed for all the simulations by using the supervision strategy.

Next, the simulations of tracking time-varying frequencies under different SNRs are performed. The standard deviation (SD) ${\tilde{\omega }}_i$ of the estimation bias for 100 simulations is shown in

Table 3. Standard deviation of estimation bias under different SNRs.

SNR	${\tilde{\mathit{\omega }}}_1$ (rad/s)	${\tilde{\mathit{\omega }}}_2$ (rad/s)	${\tilde{\mathit{\omega }}}_3$ (rad/s)	${\tilde{\mathit{\omega }}}_1$ (rad/s)	${\tilde{\mathit{\omega }}}_2$ (rad/s)	${\tilde{\mathit{\omega }}}_3$ (rad/s)
3 dB	1.6088	1.5558	1.5886	(0.9337)	(0.9784)	(0.9148)
0 dB	1.8664	1.8611	1.8324	(1.2929)	(1.3841)	(1.2931)
−10 dB	2.8658	2.8936	2.7436	(2.5265)	(2.6377)	(2.3819)

. The SD with and without the bracket correspond to the frequency variation rate of 2 $\mathrm{rad}/{\mathrm{s}}^2$ and 5 $\mathrm{rad}/{\mathrm{s}}^2$, respectively. It is obvious that the estimation precision degrades with decreasing the SNR and increasing the variation rate of the frequency. When there is no parameter supervision, some estimations may be false due to the noise. However, the estimations can track the frequencies accurately in all scenarios with the addition of the supervision strategy.

6.2. Parameter Tuning and Stability Robustness Analysis Results

Suppose that the required closed-loop bandwidth is 100 rad/s, and then, $k_1$ and $k_2$ can be selected as $k_1=20b_0$ and $k_2=0.2b_0$, respectively. Since the velocity signal can be measured directly, the RGESO is designed based on the first-order ESO as follows:

$$\begin{array}{c}\hfill G_f^{\widehat{f}}=\frac{{\omega }_o+{\prod }_{i=1}^3{R}_i\left(s\right)}{(s+{\omega }_o)+{\prod }_{i=1}^3{R}_i\left(s\right)}\end{array}$$

(60)

The problem in (52) is solved using the interior point method, and the optimized control parameters are presented in

Table 4. Parameters of the RGESO- and ${\mathrm{GESO}}_i$-based control methods.

RGESO	${\mathbf{GESO}}_1$	${\mathbf{GESO}}_2$	${\mathbf{GESO}}_3$	${\mathbf{GESO}}_4$
$k_1=20b_0$, $k_2=0.2b_0$, ${\omega }_o=712.19,$	$k_1=20b_0$,	$k_1=20b_0$,	$k_1=20b_0$,	$k_1=20b_0$,
${\zeta }_{1n}=1.00\times {10}^6,{\zeta }_{2n}=33.85,{\zeta }_{3n}=95.61,$	$k_2=0.2b_0$,	$k_2=0.2b_0$,	$k_2=0.2b_0$,	$k_2=0.2b_0$,
${\zeta }_{1d}=2.21,{\zeta }_{2d}=0.23,{\zeta }_{3d}=0.36$	${\omega }_o=5.00\times {10}^5$	${\omega }_o=1.00\times {10}^6$	${\omega }_o=3.15\times {10}^4$	${\omega }_o=9.52\times {10}^3$

. The closed-loop performance transfer functions $G_r^y$ and $G_d^y$ are shown in Figure 3. It can be seen that the required tracking bandwidth of 100 rad/s is realized, and the disturbance rejection specifications in (52) are all met. Due to the decoupling characteristic of the 2-DOF scheme, $G_r^y$ is almost the same with the desired closed-loop dynamics in (47), even though the error feedback gains are selected first without considering the dynamics of the RGESO. Figure 4 and

Table 5. Stability margin and crossover frequency for different values of ${\omega }_o$.

${\mathit{\omega }}_{\mathit{o}}$ (rad/s)	2000	1000	900	800	712.19
$G_m$ (Gain Margin, dB)	−21.28	−14.07	−13.05	−11.93	−10.85
$P_m$ (Phase Margin, deg)	77.12	79.79	80.24	80.73	81.21
$t_m$ (Time-Delay Margin, ms)	0.2057	0.2572	0.2643	0.2718	0.2790
$W_c$ (Crossover Frequency, rad/s)	6544	5415	5299	5183	5080

show the effect of ${\omega }_o$ on the closed-loop performance. According to Figure 4, the tracking performance for different values of ${\omega }_o$ is almost the same when ${\omega }_o$ is large enough. For the disturbance rejection performance, the frequency responses around the narrow-band frequencies are very close, and they are mainly determined by the parameters of the resonant term, while in other frequencies, a higher ${\omega }_o$ can achieve better disturbance attenuation performance. The expense of a higher ${\omega }_o$ is a higher crossover frequency and thus a smaller phase margin (as well as the time-delay tolerance), as can be seen in Table 5. Since the optimized ${\omega }_o$ can effectively reject the wide-band uncertainty and the narrow-band disturbance, a higher ${\omega }_o$ is unnecessary and may degrade the robustness.

Next, the verification of the stability robustness within the variation range $\left[-10\%,10\%\right]$ of the frequency is performed based on the above-mentioned graphical method. The closed-loop characteristic polynomial can be obtained as

$$\begin{array}{c}\hfill G_c(s,\mathit{\omega })=R_d\left(s\right)C_1\left(s\right)+R_n{C}_2\left(s\right)\end{array}$$

(61)

where

$$\begin{array}{cc}\hfill R_q\left(s\right)=& s^6+\sum _{i=0}^5{c}_{iq}{s}^i,c_{0q}=\prod _{i=1}^3{\omega }_i^2,c_{1q}=\sum _{i=1}^3\left({\zeta }_{iq}\prod _{j=1,j\ne i}^3{\omega }_j^2\right),\hfill \\ \hfill c_{2q}=& \sum _{i=1}^3\left({\omega }_i^2\prod _{j=1,j\ne i}^3{\zeta }_{jq}\right)+\sum _{i=1}^3\prod _{j=1,j\ne i}^3{\omega }_i^2,c_{3q}=\sum _{i=1}^3\left({\omega }_i^2\sum _{j=1,j\ne i}^3{\zeta }_{jq}\right)+\prod _{i=1}^3{\zeta }_{iq},\hfill \\ \hfill c_{4q}=& \sum _{i=1}^3\left({\omega }_i^2+\prod _{j=1,j\ne i}^3{\zeta }_{jq}\right),c_{5q}=\sum _{i=1}^3{\zeta }_{iq},\hfill \\ \hfill C_1\left(s\right)=& s^3+({\omega }_o+k_2+p_1)s^2+({\omega }_o{k}_2+k_1+p_0)s+{\omega }_o{k}_1,C_2\left(s\right)=s^2+k_{2}s+k_1,q=n,d.\hfill \end{array}$$

The value sets of the polynomial family are depicted in Figure 5. It can be observed that the origin is not included in the value sets, and thus, the system is robustly stable.

6.3. Comparative Analysis

To illustrate the advantages of the proposed scheme, comparisons with ${\mathrm{GESO}}_i$, whose order is i $(i=1,2,3,4)$, are performed. The error feedback gains are selected to be the same as those designed above for a fair comparison. For ${\mathrm{GESO}}_i$, the transfer function from the practical disturbance to the disturbance estimation is

$$\begin{array}{cc}\hfill G_{{\mathrm{GESO}}_1}=& \frac{{\omega }_o}{s+{\omega }_o},G_{{\mathrm{GESO}}_2}=\frac{{\omega }_o^2}{{(s+{\omega }_o)}^2},\hfill \\ \hfill G_{{\mathrm{GESO}}_3}=& \frac{3{\omega }_o^{2}s+{\omega }_o^3}{{(s+{\omega }_o)}^3},G_{{\mathrm{GESO}}_4}=\frac{6{\omega }_o^2{s}^2+4{\omega }_o^{3}s+{\omega }_o^4}{{(s+{\omega }_o)}^4}.\hfill \end{array}$$

(62)

To perform a fair comparison between the ${\mathrm{GESO}}_i$ and the proposed RGESO, the bandwidths of the ${\mathrm{GESO}}_i$ are determined to achieve specified disturbance attenuation performance (as defined in (52)) as

$$\begin{array}{c}\hfill |{\left.F_i\left(s\right)\right|}_{j{\omega }_k}|<S_k\mathrm{dB},i=1,2,3,4,k=1,2,3\end{array}$$

(63)

where $F_i\left(s\right)$ is the transfer function from the disturbance to the output. For the closed-loop systems designed based on the ${\mathrm{GESO}}_i$, $F_i\left(s\right)$ can be obtained as

$$\begin{array}{cc}\hfill & F_1\left(s\right)=\frac{b_{0}s}{s^3+\left(p_1+k_2+{\omega }_o\right)s^2+\left(p_2+k_1+k_2{\omega }_o\right)s+k_1{\omega }_o}\hfill \\ \hfill & F_2\left(s\right)=\frac{b_0{s}^2+2b_0{\omega }_{o}s}{\left[\begin{array}{c}{s}^4+(p_1+k_2+2{\omega }_o)s^3+(p_2+k_1+2p_1{\omega }_o+2k_2{\omega }_o+{{\omega }_o}^2)s^2+\hfill \\ (2p_2{\omega }_o+2k_1{\omega }_o+k_2{{\omega }_o}^2)s+k_1{{\omega }_o}^2\hfill \end{array}\right]}\hfill \\ \hfill & F_3\left(s\right)=\frac{b_0{s}^3+3b_0{\omega }_o{s}^2}{\left[\begin{array}{c}{s}^5+(p_1+k_2+3{\omega }_o)s^4+(p_2+k_1+3p_1{\omega }_o+3k_2{\omega }_o+3{{\omega }_o}^2)s^3+\hfill \\ (3p_2{\omega }_o+3k_1{\omega }_o+3k_2{{\omega }_o}^2+{{\omega }_o}^3)s^2+(k_2{{\omega }_o}^3+3k_1{{\omega }_o}^2)s+k_1{{\omega }_o}^3\hfill \end{array}\right]}\hfill \\ \hfill & F_4\left(s\right)=\frac{(b_0{s}^4+4b_0{\omega }_o{s}^3)}{\left[\begin{array}{c}(s^6+(p_1+k_2+4{\omega }_o)s^5+(p_2+k_1+4p_1{\omega }_o+4k_2{\omega }_o+6{{\omega }_o}^2)s^4+\hfill \\ (4p_2{\omega }_o+4k_1{\omega }_o+6k_2{{\omega }_o}^2+4{{\omega }_o}^3)s^3+\hfill \\ ({{\omega }_o}^4+4k_2{{\omega }_o}^3+6k_1{{\omega }_o}^2)s^2+(k_2{{\omega }_o}^4+4k_1{{\omega }_o}^3)s+k_1{{\omega }_o}^4\hfill \end{array}\right]}\hfill \end{array}$$

(64)

Because $F_i\left(s\right)$ is in the form of a high-pass filter, only the following condition for the highest narrow-band disturbance is needed to be satisfied:

$$|{\left.F_i\left(s\right)\right|}_{j{\omega }_3}|<S_3\mathrm{dB},i=1,2,3,4$$

(65)

To reduce ${\omega }_o$ as much as possible, we let

$$|{\left.F_i\left(s\right)\right|}_{j{\omega }_3}|=S_3\mathrm{dB},i=1,2,3,4$$

(66)

Based on (66), we can obtain ${\omega }_o$ for ${\mathrm{GESO}}_i$ by solving the high-order equation in terms of ${\omega }_o$. After some calculations, the required bandwidths for ${\mathrm{GESO}}_i$ $(i=1,2,3,4)$ are derived and presented in Table 4. The required bandwidth of ${\mathrm{GESO}}_4$ is the smallest, which is consistent with the previous results [28,29,30]. The frequency responses of closed-loop transfer functions from the disturbance to the output are presented in Figure 6. We can see that the disturbance attenuation requirement defined in (52) is satisfied for the proposed and the comparative methods.

For RGESO- and ${\mathrm{GESO}}_i$-based control systems, the transfer functions from the reference signal to the output are depicted in Figure 7. Due to the sufficiently high bandwidths of RESO and ${\mathrm{GESO}}_i$, the decoupling characteristic of the 2-DOF scheme is achieved, resulting in $G_r^y$ of all methods closely resembling the desired closed-loop dynamics.

To evaluate the sensitivity to the measurement noise, $G_n^u$ (the transfer function from the measurement noise to the control input) is analyzed. The frequency-domain responses of $G_n^u$ for various methods are shown in Figure 8. As evident in Figure 8, the significantly higher bandwidths of ${\mathrm{GESO}}_i$ result in a substantial amplification of the noise in the control signal, potentially accelerating actuator wear and tear. By contrast, the RESO causes the least amount of amplification on noise, indicating its superior robustness to measurement noise.

To achieve the desirable attenuation magnitude at the narrow-band frequencies, the bandwidths have to be designed to be sufficiently large. With such high bandwidths, the robustness to the unmodeled dynamics (such as the time-delay) will be degraded. In addition, a high sampling rate is required, which is expensive and often unimplementable in practice. The gain margin ($GM$), phase margin ($PM$), crossover frequency ($W_c$), and time delay tolerance $t_m$($t_m=P_m/\phantom{P_m{W}_c}{W}_c$), as well as the required Shannon sampling frequency $f_s$ are shown

Table 6. Stability margin and crossover frequency for different types of observers.

	${\mathbf{GESO}}_1$	${\mathbf{GESO}}_2$	${\mathbf{GESO}}_3$	${\mathbf{GESO}}_4$	RGESO
$G_m$ (dB)	Inf	Inf	−20.07	−10.01	−10.85
$P_m$ (deg)	89.51	75.85	51.03	43.32	81.21
$t_m$ (ms)	0.0031	0.0027	0.0258	0.0458	0.2790
$W_c$ (rad/s)	$5.04\times {10}^5$	$4.90\times {10}^5$	$3.45\times {10}^4$	$1.65\times {10}^4$	5080
$f_s$ (Hz)	$1.60\times {10}^5$	$1.56\times {10}^5$	$1.10\times {10}^4$	$5.25\times {10}^3$	$1.62\times {10}^3$

. It can be seen that increasing the order of GESO can help enhance the robustness to time delay and lower the required sampling time. However, the performance of the ${\mathrm{GESO}}_4$ is still poor compared with the RGESO.

To validate the robustness of various methods, we consider parameter perturbations and time delay. First, we performed 1000 Monte Carlo simulations, with the parameters ($p_1$, $p_2$, and $b_0$) randomly perturbed between −50% and 50% in each simulation, including the cases of extreme perturbations. The simulation results indicate that all methods are capable of maintaining closed-loop stability. Next, to simulate the effects of time delay, the VCM model is modeled as $P\left(s\right)e^{-\tau s}$, where $\tau$ represents the time delay. Our analysis revealed that the closed-loop systems become unstable when $\tau >t_m$. Notably, the proposed method exhibits the strongest robustness to a time delay.

6.4. Numerical Simulations

Numerical simulations are performed with the optimized controller. The sampling time is 0.0001 s. The disturbance frequencies are set as nominal values at the beginning and then have step changes to the maximal and minimal values at $t=6$ s and $t=9$ s, respectively. Only the PD (proportional derivative) control is employed in the first three seconds, and the disturbance estimation is added to the control signal thereafter. The simulation results are shown in Figure 9, Figure 10 and Figure 11. As can be seen in Figure 9, the disturbance rejection performance is greatly improved after using the disturbance estimation in the control signal calculation. Furthermore, Figure 9 depicts the outputs with and without the adaptation to the frequency variation, and it can be found that the adaptation is necessary to achieve better disturbance attenuation. The tracking error increases slightly with the increase in the disturbance frequency since the attenuation magnitude becomes smaller at a higher frequency. The disturbance and frequency estimations are presented in Figure 10 and Figure 11, respectively. According to Figure 10, the disturbance estimation can track the true value accurately with very small phase delay. In addition, the high-frequency noise can be filtered out effectively due to a moderate selection of ${\omega }_o$. As shown in Figure 11, the frequency estimations can track the true values fast and accurately, which provides an important precondition for effective rejection of time-varying narrow-band disturbances.

7. Conclusions

In this paper, a novel control scheme for simultaneously rejecting wide-band uncertainties and multiple narrow-band disturbances with time-varying frequencies was proposed. This control scheme is based on an RGESO. In the RGESO, the bandwidths for the wide-band uncertainty and the narrow-band disturbance are represented by separate parameters explicitly, which can be easily optimized. To estimate the frequencies of the narrow-band disturbances, a direct frequency estimation method based on the RLM algorithm and an on-line supervision strategy was provided. It was verified that the proposed frequency estimation method can achieve high convergent performance and estimation precision even when the SNR is very low. In the proposed tuning procedure, the error feedback gains were determined according to the desired tracking performance, and then, the parameters of the RGESO were optimized considering multiple constraints. Moreover, a precise yet simple graphical method was employed to evaluate the stability robustness of the control system, which has narrow-band disturbances with time-varying frequencies. Finally, extensive simulations were conducted to demonstrate the performance of the proposed control scheme.

In the future, we will establish a VCM experiment platform and thoroughly test our algorithm to ensure its efficiency.