© Wuhan University 2025

0 Introduction

Industrial robots are integral to modern automation technologies, extensively utilized in manufacturing for tasks such as assembly, welding, material handling, and precision machining. The drive system, as the core of an industrial robot, ensures precise movement and operational efficiency. However, these systems are susceptible to faults in their drive control mechanisms, which can degrade performance, reduce productivity, or even lead to catastrophic failures^[1-2]. As the adoption of industrial robots grows, fault detection and diagnosis in drive systems have become critical for ensuring reliability and minimizing operational downtime^[3].

One of the earliest approaches to improving system reliability, hardware redundancy, involves duplicating critical components to maintain operation during faults^[4-5]. Although effective, this approach increases system cost, weight, and complexity, making it less desirable for modern applications. In contrast, model-based analytical methods leverage mathematical models to represent nominal system behavior and detect anomalies by comparing predicted outputs with actual observations^[6]. These methods offer distinct advantages over hardware redundancy by enabling the detection of subtle faults and providing deeper insights into system dynamics. Among model-based methods, the Kalman filter (KF), a recursive algorithm, is widely employed to estimate the state of dynamic systems from noisy measurements. The KF has been extensively studied and applied across various fields^[7]. For instance, Lee et al^[8] developed a fault detection and diagnosis (FDD) algorithm for an open-cycle liquid propellant rocket engine using a conventional KF. Similarly, a residual generator based on the KF was proposed in Ref. [9] to diagnose faults and achieve fault-tolerant control in linear drive systems. However, despite its success, the KF's application is limited in systems exhibiting nonlinear dynamics, which are prevalent in physical and industrial systems.

To address the challenges posed by nonlinearity, advanced variants of the KF, such as the extended Kalman filter (EKF), unscented Kalman filter (UKF), and particle filter (PF), have been developed. The EKF, which linearizes nonlinear systems via Taylor series expansion, has been employed in various applications. For instance, Gao et al^[10] proposed an EKF-based fault diagnosis strategy for unmanned aerial vehicles, integrating multiple-model adaptive estimation to diagnose actuator faults efficiently. Similarly, an EKF-based fault detection model for the papermaking wastewater treatment process was introduced in Ref. [11]. However, the EKF's reliance on linearization introduces approximation errors, motivating the development of the UKF, which uses sigma-point sampling to achieve more accurate approximations of nonlinear dynamics. For example, Tian et al^[12] proposed a UKF-based online estimation algorithm for insulation resistance in battery energy storage systems. Furthermore, a distributed UKF framework for sensor fault detection, isolation, and accommodation was presented in Ref. [13]. While the particle filter offers another alternative by using particle sampling, its effectiveness can be compromised by state estimation inaccuracies^[14]. To mitigate this, Yin et al^[15] enhanced the PF with genetic algorithms to detect faults in a three-tank system. Despite these advancements, traditional KF-based FDD methods often neglect robustness and stability in state estimation, particularly in handling model uncertainties, external disturbances, and noise.

Industrial robot drive systems operate in complex environments subject to such disturbances, necessitating more robust fault detection frameworks. In this paper, we propose an enhanced UKF-based fault detection method to address these challenges. A dynamic time-varying boundary layer is incorporated into the UKF framework to balance estimation accuracy and robustness. Robust state estimates are used to generate residuals by comparing them with actual system outputs. To further improve fault detection reliability, z-score analysis, and signal-to-noise ratio (SNR) analysis are integrated into the detection process. This combined approach ensures efficient and accurate fault detection in industrial robot drive systems.

The remaining sections are structured as follows: Section 1 reviews the mathematical formulation of the drive system and the convention UKF algorithm. In Section 2, an enhanced UKF algorithm is introduced. Subsequently, Section 3 proposes the fault detection scheme. Section 4 presents the numerical simulation cases for validating the proposed method. The conclusion is drawn in Section 5.

1 Preliminary

This section presents the mathematical model of the industrial robot drive system to accurately describe its dynamic characteristics and illustrates the algorithm for state estimation using the traditional UKF.

1.1 Modeling of the Drive System

The joints of industrial robot arms require motor drives with high precision, high reliability, and high torque characteristics. The brushless DC motor, which meets these requirements, is widely used in industrial robots. A typical brushless DC motor drive system includes the motor body, sinusoidal pulse width modulation (SPWM) inverter, and current feedback loop, as shown in the structure of Fig. 1. Based on the effect of current feedback, the relationships between the system's current, voltage, and control vector, after applying the Park transformation, can be expressed as follows:

Fig. 1 Schematic diagram of brushless DC motor drive systems

$\{\begin{array}{l}{u}_q=K_p{K}_s{U}_{\mathrm{i}\mathrm{n}q}-K_p{K}_s{K}_f{i}_q,\\ u_d=K_p{K}_s{U}_{\mathrm{i}\mathrm{n}d}-K_p{K}_s{K}_f{i}_d,\end{array}$(1)

where $u_d$ and $u_q$ are the voltages applied to the $d$-axis and $q$-axis, respectively. $i_d$ and $i_q$ are the currents flowing through the $d$-axis and $q$-axis. $K_p$, $K_s$, and $K_f$ represent the current controller gain, the equivalent gain of the inverter, and the current feedback parameter. $U_{\mathrm{i}\mathrm{n}d}$ and $U_{\mathrm{i}\mathrm{n}q}$ are the reference signals for the $d$-axis and $q$-axis, respectively.

To achieve a higher torque coefficient and a bidirectional symmetric linear current-torque characteristic, the rotor position sensor is typically installed such that the reference phasor aligns with the $q$-axis, i.e., $U_{\mathrm{i}\mathrm{n}d}=\mathrm{0}$, and $U_{\mathrm{i}\mathrm{n}q}$ is numerically equal to the control input $U_{\mathrm{i}\mathrm{n}}$. Additionally, for a motor body, its state can be expressed as follows.

$\{\begin{array}{l}\frac{\mathrm{d}{i}_q}{\mathrm{d}t}=(-r{i}_q-P{\omega }_{r}L{i}_d-{\varphi }_{f}P{\omega }_r+u_q)/L,\\ \frac{\mathrm{d}{i}_d}{\mathrm{d}t}=(-r{i}_d+P{\omega }_{r}L{i}_q)/L,\\ \frac{\mathrm{d}{\omega }_r}{\mathrm{d}t}=\frac{\mathrm{3}}{\mathrm{2}J}P{\varphi }_f{i}_q/L,\end{array}$(2)

where $r$ and $L$ are the winding resistance and self-induction coefficient. $J$ is the rotor inertia. $P$ stands for the number of poles. ${\varphi }_f$ is the flux linkage generated by permanent magnetic poles and the stator winding. While ${\omega }_r$ indicates the motor speed.

Considering the accuracy requirements of the mathematical model, the effect of the no-load torque should not be neglected. Combining the previous descriptions with equations (1) and (2), the final mathematical model of the brushless DC motor in the $d$-$q$ coordinate frame is expressed as follows:

$\{\begin{array}{l}\begin{array}{l}\dot{\mathit{x}}=\left[\begin{array}{ccc}-\frac{r+K_p{K}_s{K}_f}{L}& \mathrm{0}& -\frac{P{\varphi }_f}{L}\\ \mathrm{0}& -\frac{r+K_p{K}_s{K}_f}{L}& \mathrm{0}\\ \frac{\mathrm{3}P{\varphi }_f}{\mathrm{2}J}& \mathrm{0}& -\frac{B}{J}\end{array}\right]\mathit{x}\\ +\left[\begin{array}{cc}\frac{K_p{K}_s}{L}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\frac{\mathrm{1}}{J}\end{array}\right]\mathit{u}+\left[\begin{array}{c}-P{\omega }_r{i}_d\\ P{\omega }_r{i}_q\\ \mathrm{0}\end{array}\right],\end{array}\\ \mathit{y}={\mathit{I}}_{\mathrm{3}\times \mathrm{3}}\mathit{x}+\mathit{\nu },\end{array}$(3)

where $\mathit{x}=[i_q,i_d,{\omega }_r{]}^{\mathrm{T}}$ is defined as the state vector. $\mathit{u}=[U_{\mathrm{i}\mathrm{n}},T_l{]}^{\mathrm{T}}$ represents the control vector consisting of the input reference voltage and load torque. $B$ is the rotor damping coefficient. ${\mathit{I}}_{\mathrm{3}\times \mathrm{3}}$ is the identity matrix. $\mathit{\nu }$ denotes the measurement noise. Clearly, equation (3) presents a fully measurable nonlinear system. In the UKF, the unscented transform can effectively approximate the characteristics of the target system.

1.2 Estimation with the Conventional UKF

By using sigma points to propagate the state distribution through nonlinear transformations, the conventional UKF can handle complex nonlinearities without the need for linearization, making it particularly suitable for applications in industry robotics, where accurate state estimation is critical. The UKF consists of three main phases: sigma point generation, prediction step, and update step, which is iteratively applied to estimate the state and its covariance.

a) Sigma point generation

At each time step, the UKF begins by generating a set of sigma points ${\mathit{\chi }}_{k-\mathrm{1}}^{(i)}$ based on the current state estimate ${\widehat{\mathit{x}}}_{k-\mathrm{1}}$ and the associated covariance matrix ${\mathit{P}}_{k-\mathrm{1}}$.

$\{\begin{array}{l}{\mathit{\chi }}_{k-\mathrm{1}}^{(\mathrm{0})}={\widehat{\mathit{x}}}_{k-\mathrm{1}},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i=\mathrm{0},\\ {\mathit{\chi }}_{k-\mathrm{1}}^{(i)}={\widehat{\mathit{x}}}_{k-\mathrm{1}}+\sqrt[]{(n+\lambda ){\mathit{P}}_{k-\mathrm{1}}},\hspace{1em}i=\mathrm{1},\cdots ,n,\\ {\mathit{\chi }}_{k-\mathrm{1}}^{(i)}={\widehat{\mathit{x}}}_{k-\mathrm{1}}-\sqrt[]{(n+\lambda ){\mathit{P}}_{k-\mathrm{1}}},\hspace{1em}i=n+\mathrm{1},\cdots ,\mathrm{2}n,\end{array}$(4)

where $n$ is the dimension of $\mathit{x}$. ${\mathit{P}}_{k-\mathrm{1}}$ denotes the state covariance matrix. $\lambda$ is a scaling parameter that controls the spread of the sigma points around the mean state estimate.

b) Prediction step

In the prediction phase, the $\mathrm{2}n$ sigma points are propagated through the nonlinear process model, as presented in equation (3). For ease of reading, the target system (3) is represented in a compact form $f({\mathit{x}}_k,{\mathit{u}}_k)$. Thus, one has

${\mathit{\chi }}_{k|k-\mathrm{1}}^{(i)}=f({\mathit{\chi }}_{k-\mathrm{1}}^{(i)},{\mathit{u}}_k)$(5)

where ${\mathit{\chi }}_{k|k-\mathrm{1}}^{(i)}$ denotes the predicted sigma points at time $k$. Subsequently, the predicted state estimate and covariance can be obtained as the weighted average of the sigma points, as shown in equation (6).

$\{\begin{array}{c}\begin{array}{l}{\widehat{\mathit{x}}}_{k|k-\mathrm{1}}={\displaystyle \sum _{i=\mathrm{0}}^{\mathrm{2}n}}{W}_i^m{\mathit{\chi }}_{k|k-\mathrm{1}}^{(i)},\\ {\mathit{P}}_{k|k-\mathrm{1}}={\displaystyle \sum _{i=\mathrm{0}}^{\mathrm{2}n}}{W}_i^c({\mathit{\chi }}_{k|k-\mathrm{1}}^{(i)}-{\widehat{\mathit{x}}}_{k|k-\mathrm{1}})({\mathit{\chi }}_{k|k-\mathrm{1}}^{(i)}-{\widehat{\mathit{x}}}_{k|k-\mathrm{1}}{)}^{\mathrm{T}}+{\mathit{Q}}_k,\end{array}\end{array}$(6)

where $W_i^m$ and $W_i^c$ are the weights for the mean and covariance, respectively. ${\mathit{Q}}_k$ denotes the process noise covariance.

c) Update step

In the update phase, the predicted sigma points are passed through the nonlinear measurement model (not including the measurement noise), denoted as $h({\mathit{x}}_k)$, resulting in the measurement and prediction similarly:

$\{\begin{array}{l}{\xi }_{k|k-\mathrm{1}}^{(i)}=h({\mathit{\chi }}_{k|k-\mathrm{1}}^{(i)}),{\widehat{\mathit{y}}}_{k|k-\mathrm{1}}={\displaystyle \sum _{i=\mathrm{0}}^{\mathrm{2}n}}{W}_i^m{\xi }_{k|k-\mathrm{1}}^{(i)},\\ {\mathit{S}}_{k|k-\mathrm{1}}={\displaystyle \sum _{i=\mathrm{0}}^{\mathrm{2}n}}{W}_i^c({\xi }_{k|k-\mathrm{1}}^{(i)}-{\widehat{\mathit{y}}}_{k|k-\mathrm{1}})({\xi }_{k|k-\mathrm{1}}^{(i)}-{\widehat{\mathit{y}}}_{k|k-\mathrm{1}}{)}^{\mathrm{T}}+{\mathit{R}}_k,\\ {\mathit{P}}_{xy}={\displaystyle \sum _{i=\mathrm{0}}^{\mathrm{2}n}}{W}_i^c({\mathit{\chi }}_{k|k-\mathrm{1}}^{(i)}-{\widehat{\mathit{x}}}_{k|k-\mathrm{1}})({\xi }_{k|k-\mathrm{1}}^{(i)}-{\widehat{\mathit{y}}}_{k|k-\mathrm{1}}{)}^{\mathrm{T}},\end{array}$(7)

where ${\mathit{R}}_k$ denotes the measurement noise covariance. ${\mathit{S}}_{k|k-\mathrm{1}}$ denotes the predicted measurement covariance. ${\mathit{P}}_{xy}$ is the cross-covariance matrix. Based on (7), the innovation gain ${\mathit{K}}_k$ is readily calculated as:

${\mathit{K}}_k={\mathit{P}}_{xy}{\mathit{S}}_{k|k-\mathrm{1}}^{-\mathrm{1}}.$(8)

Finally, the state estimate ${\widehat{\mathit{x}}}_k$ and its covariance ${\mathit{P}}_k$ are updated as follows:

${\widehat{\mathit{x}}}_k={\widehat{\mathit{x}}}_{k|k-\mathrm{1}}+{\mathit{K}}_k({\mathit{y}}_k-{\widehat{\mathit{y}}}_{k|k-\mathrm{1}}),\hspace{1em}{\mathit{P}}_k={\mathit{P}}_{k|k-\mathrm{1}}-{\mathit{K}}_k{\mathit{S}}_{k|k-\mathrm{1}}{\mathit{K}}_k^{\mathrm{T}},$(9)

where ${\mathit{y}}_k$ denotes the actual measurement at time $k$.

2 Enhanced UKF for Robust Estimation

As seen from equations (8) and (9) in the UKF algorithm, the state estimation is achieved by using the innovation gain ${\mathit{K}}_k$. This approach typically yields high estimation accuracy. Nonetheless, ensuring reliable state estimation despite system uncertainties or external noise continues to be a significant challenge.

In addressing the trade-off between accurate estimation and robustness against disturbances, we are inspired by the variable structure filters presented in Refs. [16-17]. We incorporate a sliding boundary layer into the UKF gain update process to enhance the robustness of the estimation. The enhanced innovation gain is defined as follows:

${\mathit{K}}_k={\mathit{H}}^{†}\mathrm{s}\mathrm{a}\mathrm{t}\left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\frac{\left|{\tilde{\mathit{y}}}_k\right|}{\mathit{\epsilon }})\right),\mathrm{f}\mathrm{o}\mathrm{r}{\tilde{\mathit{y}}}_k={\mathit{y}}_k-{\widehat{\mathit{y}}}_{k|k-\mathrm{1}},$(10)

where $\mathit{\epsilon }$ is the sliding boundary layer to be designed. ${\mathit{H}}^{†}$ denotes the pseudoinverse of the measurement matrix $h(\cdot )$. $\left|\cdot \right|$, $\mathrm{s}\mathrm{a}\mathrm{t}(\cdot )$, and $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\cdot )$ represent the absolute, saturation, and diagonal operations, respectively. By replacing the conventional innovation gain ${\mathit{K}}_k$ in (8) with the enhanced one in (10), the state estimation in the UKF is forced to converge within a boundary layer $\mathit{\epsilon }$.

Clearly, the selection of the boundary layer width is crucial. The simplest approach is to set a fixed width based on the system's maximum uncertainty. However, this method lacks flexibility. Therefore, we next focus on designing an adaptive time-varying boundary layer to adapt to the dynamic estimation process.

Based on (10), and utilizing a first-order Taylor approximation, one has the posterior error as follows:

$\begin{array}{l}{\mathit{e}}_k={\widehat{\mathit{x}}}_k-{\mathit{x}}_k={\widehat{\mathit{x}}}_{k|k-\mathrm{1}}-{\mathit{x}}_k+{\mathit{K}}_k({\mathit{y}}_k-{\widehat{\mathit{y}}}_{k|k-\mathrm{1}})\\ ={\widehat{\mathit{e}}}_k+{\mathit{K}}_k[h({\mathit{x}}_k)-h({\widehat{\mathit{x}}}_{k|k-\mathrm{1}})]+{\mathit{K}}_k{\mathit{\nu }}_k\\ \approx {\widehat{\mathit{e}}}_k+{\mathit{K}}_k\frac{\partial h}{\partial \mathit{x}}({\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k-\mathrm{1}})+{\mathit{K}}_k{\mathit{\nu }}_k\\ ={\widehat{\mathit{e}}}_k-{\mathit{K}}_k\mathit{H}{\widehat{\mathit{e}}}_k+{\mathit{K}}_k{\mathit{\nu }}_k\\ =(\mathit{I}-{\mathit{K}}_k\mathit{H}){\widehat{\mathit{e}}}_k+{\mathit{K}}_k{\mathit{\nu }}_k,\end{array}$(11)

where ${\widehat{\mathit{e}}}_k={\widehat{\mathit{x}}}_{k|k-\mathrm{1}}-{\mathit{x}}_k$ denotes the prior error. $\mathit{H}=\partial h/\partial \mathit{x}$ is the Jacobian matrix with respect to $h(\cdot )$. Consequently, the state error covariance matrix ${\mathit{P}}_k$ is readily calculated as (12), which indicates the characteristics of the estimated state.

${\mathit{P}}_k=E({\mathit{e}}_k{\mathit{e}}_k^{\mathrm{T}})=(\mathit{I}-{\mathit{K}}_k\mathit{H}){\widehat{\mathit{P}}}_k(\mathit{I}-{\mathit{K}}_k{\mathit{H})}^{\mathrm{T}}+{\mathit{K}}_k{\mathit{R}}_k{\mathit{K}}_k^{\mathrm{T}}$(12)

${\widehat{\mathit{P}}}_k={\widehat{\mathit{e}}}_k{\widehat{\mathit{e}}}_k^{\mathrm{T}}$ is the prior state error covariance. ${\mathit{R}}_k={\mathit{\nu }}_k{\mathit{\nu }}_k^{\mathrm{T}}$ represents the measurement noise covariance as described before. It is evident that ${\mathit{P}}_k$ is a function of the boundary layer parameter $\mathit{\epsilon }$, and thus, the optimal value of $\mathit{\epsilon }$ can be determined by minimizing the trace of ${\mathit{P}}_k$, i.e.,

$\mathit{\epsilon }=\underset{\mathit{\epsilon }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}[\mathrm{t}\mathrm{r}({\mathit{P}}_k)]$(13)

For ease of calculation, the enhanced innovation gain in (10) is denoted in a compact form ${\mathit{K}}_k={\mathit{H}}^{†}\mathit{Y}{\mathit{\epsilon }}^{-\mathrm{1}}$, where $\mathit{Y}$ and ${\mathit{\epsilon }}^{-\mathrm{1}}$ represent the saturation of associated diagonal elements in (10), respectively. Therefore, by combining (12), the optimal solution to (13) is calculated as follows:

$\begin{array}{l}\partial [\mathrm{t}\mathrm{r}({\mathit{P}}_k)]/\partial \mathit{\epsilon }\\ =\partial \{\mathrm{t}\mathrm{r}[{\widehat{\mathit{P}}}_k-{\mathit{H}}^{†}\mathit{Y}{\mathit{\epsilon }}^{-\mathrm{1}}\mathit{H}{\widehat{\mathit{P}}}_k-{\widehat{\mathit{P}}}_k{\mathit{H}}^{\mathrm{T}}{\mathit{\epsilon }}^{-\mathrm{T}}{\mathit{Y}}^{\mathrm{T}}{\mathit{H}}^{†\mathrm{T}}\\ +{\mathit{H}}^{†}\mathit{Y}{\mathit{\epsilon }}^{-\mathrm{1}}{\mathit{\Gamma }}_k{\mathit{\epsilon }}^{-\mathrm{T}}{\mathit{Y}}^{\mathrm{T}}{\mathit{H}}^{†\mathrm{T}}]\}/\partial \mathit{\epsilon }\\ =\partial [\mathrm{t}\mathrm{r}({\widehat{\mathit{P}}}_k)-\mathrm{t}\mathrm{r}({\mathit{H}}^{†}\mathit{Y}{\mathit{\epsilon }}^{-\mathrm{1}}\mathit{H}{\widehat{\mathit{P}}}_k)-\mathrm{t}\mathrm{r}({\widehat{\mathit{P}}}_k{\mathit{H}}^{\mathrm{T}}{\mathit{\epsilon }}^{-\mathrm{T}}{\mathit{Y}}^{\mathrm{T}}{\mathit{H}}^{†\mathrm{T}})\\ +\mathrm{t}\mathrm{r}({\mathit{H}}^{†}\mathit{Y}{\mathit{\epsilon }}^{-\mathrm{1}}{\mathit{\Gamma }}_k{\mathit{\epsilon }}^{-\mathrm{T}}{\mathit{Y}}^{\mathrm{T}}{\mathit{H}}^{†\mathrm{T}})]/\partial \mathit{\epsilon }\\ =\mathrm{0},\end{array}$(14)

where ${\mathit{\Gamma }}_k=\mathit{H}{\widehat{\mathit{P}}}_k{\mathit{H}}^{\mathrm{T}}+{\mathit{R}}_k$. By applying the matrix trace operation rules to the terms in (14), the equation is simplified as:

${\mathit{\epsilon }}^{-\mathrm{T}}({\mathit{Y}}^{\mathrm{T}}{\mathit{H}}^{†\mathrm{T}}{\widehat{\mathit{P}}}_k{\mathit{H}}^{\mathrm{T}}-{\mathit{Y}}^{\mathrm{T}}{\mathit{H}}^{†\mathrm{T}}{\mathit{H}}^{†}\mathit{Y}{\mathit{\epsilon }}^{-\mathrm{1}}{\mathit{\Gamma }}_k){\mathit{\epsilon }}^{-\mathrm{T}}=\mathrm{0}$(15)

For (15), the solution can be easily obtained as:

$\mathit{\epsilon }={\mathit{\Gamma }}_k(\mathit{H}{\widehat{\mathit{P}}}_k{\mathit{H}}^{\mathrm{T}}{)}^{-\mathrm{1}}\mathit{Y}$(16)

Combining equations (9), (10), and (16), one obtains the new state estimation with the enhanced innovation gain, as expressed in equation (17).

${\widehat{\mathit{x}}}_k={\widehat{\mathit{x}}}_{k|k-\mathrm{1}}+{\mathit{H}}^{†}\mathrm{s}\mathrm{a}\mathrm{t}\left(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\frac{\left|{\tilde{\mathit{y}}}_k\right|}{\mathit{\epsilon }})\right)({\mathit{y}}_k-{\widehat{\mathit{y}}}_{k|k-\mathrm{1}})$(17)

Obviously, by introducing a switching term in the state estimation, the estimated state is driven to converge within the boundary layer. Unlike traditional smooth variable structure filters (SVSF, refer to Refs. [16-17]), the time-varying boundary layer proposed in this paper employs a saturation function, denoted as $\mathrm{s}\mathrm{a}\mathrm{t}(\cdot )$, instead of a sign function, denoted as $\mathrm{s}\mathrm{g}\mathrm{n}(\cdot )$, which mitigates the chattering phenomenon, as illustrated in Fig. 2. Note that $\mathit{\epsilon }$ incorporates both the state covariance and measurement error (see equation (16)), which presents the system uncertainties and measurement noise during the estimation process. This is precisely why the proposed method enhances the robustness of the estimation.

Fig. 2 State estimation comparison with different methods

3 Fault Detection Strategy

As mentioned earlier, fault detection is crucial when anomalies occur in the industrial robot drive system. The generation and analysis of residuals serve as the fundamental basis for fault detection. This section will focus on the fault detection framework and residual analysis methods based on the enhanced UKF.

3.1 Fault Detection Framework

The fault detection process based on the enhanced UKF is initiated by estimating the system state using the enhanced UKF algorithm detailed in Section 2. This state estimation provides a reference for subsequent analysis. Subsequently, residuals are generated by comparing the actual system outputs with the estimates from the enhanced UKF model. These residuals encapsulate the discrepancies between the actual and expected behaviors of the system, serving as primary indicators of potential faults. The principle framework for fault detection is shown in Fig. 3.

Fig. 3 Principle fault detection framework based on the enhanced UKF

Specifically, to ensure the reliability and effectiveness of the fault detection framework, a comprehensive monitoring system is established. As shown in Fig. 3, the fault detection system tracks performance metrics of the industrial robot drive system, such as motor speed ${\omega }_r$ and current fluctuations of $i_d$ and $i_q$, as illustrated in equation (3). By integrating these monitored data with the residual analysis, a more holistic understanding of the system's health status is achieved.

3.2 Residual Analysis

This paper first applies the improved z-score method for normalization. The z-score process is given by the following formula:

${\mathit{Z}}_k=\frac{\mathrm{e}\mathrm{r}{\mathrm{r}}_k-\mathrm{m}\mathrm{e}\mathrm{d}(\mathrm{e}\mathrm{r}{\mathrm{r}}_{(k-m+\mathrm{1})}^{(k)})}{\mathrm{M}\mathrm{A}\mathrm{D}},\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{M}\mathrm{A}\mathrm{D}=\mathrm{m}\mathrm{e}\mathrm{d}(\mathrm{e}\mathrm{r}{\mathrm{r}}_k-\mathrm{m}\mathrm{e}\mathrm{d}(\mathrm{e}\mathrm{r}{\mathrm{r}}_{(k-m+\mathrm{1})}^{(k)}))$(18)

where $\mathrm{m}\mathrm{e}\mathrm{d}(\cdot )$ represents the median operation. $\mathrm{e}\mathrm{r}{\mathrm{r}}_k$ is the residual vector at time $k$ (shown as "residual" in Fig.3). Additionally, $\mathrm{e}\mathrm{r}{\mathrm{r}}_{(k-m+\mathrm{1})}^{(k)}$ represents a residual set consisting of a window of residual samples with the length of $m$. Without using the insensitive mean and standard deviation in the traditional z-score, the improved z-score is based on the median and median absolute deviation, providing a better measure of central tendency and dispersion for skewed or non-normal data^[18]. Therefore, the improved z-score offers an effective way of indicating potential anomalies in the residual signals.

Considering the measurement noise or unknown external disturbances affecting the residuals, misleading results of the improved z-score, we treat the normalized results ${\mathit{Z}}_k$ as noise-contaminated signals and apply the SNR method for filtering to detect anomalies, calculated as follows.

$P_{\mathrm{S}\mathrm{N}\mathrm{R}}=\frac{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([{\mathit{Z}}_k-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}({\mathit{Z}}_{(k-m+\mathrm{1})}^{(k)})][{\mathit{Z}}_k-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}({\mathit{Z}}_{(k-m+\mathrm{1})}^{(k)}{)]}^{\mathrm{T}})}{{\sigma }^{\mathrm{2}}({\mathit{Z}}_{(k-m+\mathrm{1})}^{(k)})},$(19)

where $\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}(\cdot )$ and ${\sigma }^{\mathrm{2}}(\cdot )$ represent the mean operation and variance operation, respectively. ${\mathit{Z}}_{(k-m+\mathrm{1})}^{(k)}$ is a set of ${\mathit{Z}}_k$ within a window of $m$ samples.

The SNR metric in equation (19) can be understood as the ratio of the power of specific correlated noise to the average power of the entire sampling noise. When the system is operating normally, the SNR metric has a lower value. However, when a fault occurs, the PSNR metric will display higher amplitudes, as illustrated in the following rule:

$\{\begin{array}{cc}{P}_{\mathrm{S}\mathrm{N}\mathrm{R}}\ge {\tau }_{\mathrm{t}\mathrm{h}},& \mathrm{f}\mathrm{a}\mathrm{u}\mathrm{l}\mathrm{t},\\ P_{\mathrm{S}\mathrm{N}\mathrm{R}}<{\tau }_{\mathrm{t}\mathrm{h}},& \mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l},\end{array}$(20)

where ${\tau }_{\mathrm{t}\mathrm{h}}$ represents the threshold, selected as the mean of $P_{\mathrm{S}\mathrm{N}\mathrm{R}}$ over a window of samples from the initial fault-free moment.

4 Numeral Simulation

This section validates the effectiveness of the proposed method (denoted as EUKF) through two simulation cases and compares it with the traditional UKF method in Ref. [19] (denoted as UKF) and the adaptive robust UKF method proposed in Ref. [20] (denoted as ARUKF). To initiate the simulation, the parameters of the brushless DC motor are configured as follows. $K_p=K_s=K_f=\mathrm{1}$, $r=\mathrm{0.2}\mathrm{\Omega }$, $L=\mathrm{0.009}\mathrm{H}$, $B=\mathrm{0.008}\mathrm{N}\cdot \mathrm{m}\cdot \mathrm{s}$, ${\phi }_f=\mathrm{0.193}\mathrm{W}\mathrm{b}$, $J=\mathrm{0.001}\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{\mathrm{2}}$, and $T_l=\mathrm{0.05}\mathrm{N}\cdot \mathrm{m}$. The measurement noise is Gaussian white noise with a power level of $-\mathrm{30}\mathrm{d}\mathrm{B}$. Moreover, the sampling step is configured as $t_s=\mathrm{0.002}\mathrm{s}$.

4.1 Case 1: Robust Estimation Performance under Fault-Free Conditions

To validate the robustness of the proposed method, in this case, the system is subjected to Gaussian white noise interference without any faults occurring, and the initial value of the system state is set to ${\mathit{x}}_{\mathrm{0}}={[\mathrm{0.01,0.01,0}]}^{\mathrm{T}}$. To save space, we only plot the $I_q$ state output and the estimation curves of different methods in the left subplot of Fig. 4, where the actual output is labelled as $I_q$, and state estimates of the different methods are labelled with EUKF, ARUKF, and UKF. To clearly illustrate the estimation errors of different methods under noise interference, we present the estimation error curves (as labelled with Err) and box plots for the different methods on the right of Fig. 4. Furthermore, for ease of comparison, RMSE statistics of different estimation errors are provided in Table 1. The statistics show that, for all state estimations, the proposed method EUKF exhibits the smallest RMSE compared to the other methods, demonstrating the robustness of the proposed enhanced UKF.

Fig. 4 State estimations and estimation errors with different methods for ${\mathit{I}}_{\mathit{q}}$

4.2 Case 2: Fault Detection under Fault Conditions

Due to the complex working environment and continuous long-term operation of industrial robots, the drive system is prone to failures, such as increased operating temperature and lubrication failure, which can lead to sudden abnormalities in the motor rotor damping parameter $B$. In this case, it is assumed that $B$ changes from $\mathrm{0.008}\mathrm{N}\cdot \mathrm{m}\cdot \mathrm{s}$ to $\mathrm{0.01}\mathrm{N}\cdot \mathrm{m}\cdot \mathrm{s}$ between 2.6 s and 2.63 s. Residuals are generated using the proposed method to perform fault detection. Similarly, the state estimation and fault detection results for $I_q$ are illustrated in Fig. 5.

Fig. 5 State estimation and fault detection for ${\mathit{I}}_{\mathit{q}}$ under the fault condition

From the state estimation results in Fig. 5 (shown in the left subplot), it can be observed that when $B$ experiences an anomaly, the system state $I_q$ deviates from its normal value, and the estimation results of all methods exhibit fluctuations. The SNR results for each comparative method are shown in the three subplots on the right,

labelled as P_SNR-EUKF, P_SNR-ARUKF, and P_SNR-UKF. The corresponding thresholds are labelled as τ-EUKF, τ-ARUKF, and τ-UKF. When an anomaly occurs, the SNR statistics can quickly detect the system anomaly and exceed the anomaly threshold. To visually compare the fault detection results of different methods, we statistically analyze the false detection rate (FDR) and missed detection rate (MDR) for each state (i.e., $I_q$, $I_d$, and $\omega {}_r{}^{}$) with different methods, as shown in Table 2. The data show that EUKF, compared with the other methods, achieves the lowest FDR and MDR across all states. This indicates that the proposed state estimation and residual analysis strategy provides the best performance for fault detection in the system.

5 Conclusion

This paper presents an enhanced unscented Kalman filter (EUKF) approach for fault detection in industrial robot drive systems. By incorporating a dynamic, time-varying boundary layer into the UKF and using a combination of improved z-score and SNR analysis for residual evaluation, the proposed method achieves a better trade-off between estimation accuracy and robustness. The numerical simulation results under fault-free and faulty conditions suggest that the proposed approach can enhance the reliability and efficiency of fault detection in industrial robot drive systems.

Future work will focus on extending the proposed method to handle more complex fault scenarios and integrating it with fault-tolerant control strategies to improve the overall performance and safety of industrial robots.

References

All Tables

All Figures

	Fig. 1 Schematic diagram of brushless DC motor drive systems
In the text

	Fig. 2 State estimation comparison with different methods
In the text

	Fig. 3 Principle fault detection framework based on the enhanced UKF
In the text

	Fig. 4 State estimations and estimation errors with different methods for ${\mathit{I}}_{\mathit{q}}$
In the text

	Fig. 5 State estimation and fault detection for ${\mathit{I}}_{\mathit{q}}$ under the fault condition
In the text