2.1 Establishment of degradation model for quadrotor UAV

For a quadrotor UAV shown in Fig. 1, ignoring wind interference, the Newton-Euler formula is used to model it in a hybrid coordinate system as shown in Equation (1) [Reference Zhang, Chamseddine, Rabbath, Gordon, Su, Rakheja, Fulford, Apkarian and Gosselin27].

(1) \begin{align} & \ddot x = {{\left( {\cos \phi \sin \theta \cos \psi + \sin \phi \sin \psi } \right){U_z}} / m} \nonumber\\[4pt] & \ddot y = {{\left( {\cos \phi \sin \theta \sin \psi - \sin \phi \cos \psi } \right){U_z}} / m} \nonumber\\[4pt] & \ddot z = {{\left( {\cos \phi \cos \theta } \right){U_z}} / m} - g \nonumber\\[4pt] & \ddot \phi = {{\left( {\dot \theta \dot \psi \left( {{I_{yy}} - {I_{zz}}} \right) + {U_\phi }} \right)} / {{I_{xx}}}}\\[4pt] & \ddot \theta = {{\left( {\dot \phi \dot \psi \left( {{I_{zz}} - {I_{xx}}} \right) + {U_\theta }} \right)} / {{I_{yy}}}} \nonumber\\[4pt] & \ddot \psi = {{\left( {\dot \theta \dot \phi \left( {{I_{xx}} - {I_{yy}}} \right) + {U_\psi }} \right)} / } {{I_{zz}}} \nonumber \end{align}

where $x,y,z$ are the actual positions of the quadrotor UAV in the inertial frame with the original point at the GC, $\phi ,\theta ,\psi $ are the roll angle, pitch angle and yaw angle, respectively. The parameters $m,{\rm{g}}$ refer to the total airframe mass and gravitational acceleration of the UAV. ${I_{xx}},{I_{yy}},{I_{zz}}$ are the inertial moments about each axis. ${U_z},{U_\phi },{U_\theta },{U_\psi }$ are utilised as control inputs for position regulation. Equation (2) demonstrates the correlation between UAV and motor inputs ${u_i}$ , while ${U_z}$ is employed to regulate the position in $x,y,z$ directions. Additionally, ${U_\phi },{U_\theta },{U_\psi }$ are implemented to govern roll, pitch and yaw respectively.

(2) \begin{equation}\left[ {\begin{array}{*{20}{c}}{{U_z}}\\[4pt]{{U_\theta }}\\[4pt]{{U_\phi }}\\[4pt]{{U_\psi }}\end{array}} \right] = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{K_u}} & {{K_u}} & {{K_u}} & {{K_u}}\\[4pt]{{K_u}{L_d}} & { - {K_u}{L_d}} & 0 & 0\\[4pt]0 & 0 & {{K_u}{L_d}} & { - {K_u}{L_d}}\\[4pt]{{K_y}} & {{K_y}} & { - {K_y}} & { - {K_y}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{u_1}}\\[4pt]{{u_2}}\\[4pt]{{u_3}}\\[4pt]{{u_4}}\end{array}} \right]\end{equation}

The variables ${u_i}(i = 1,2,3,4)$ represent the pulse-width modulation (PWM) input of the ith motor; ${L_d}$ denotes the distance from the quadrotor UAV motor to the centre of mass; and ${K_u}$ is the thrust gain related to the force generated by the propeller, ${K_y}$ is the torque gain related to the torque generated by the propeller.

Assuming that the quadrotor UAV is hovering with no yaw ( $\psi = 0$ ), small roll and pitch angles, ignoring the drag term, gyroscopic and Coriolis effects, Equation (1) is simplified as:

(3) \begin{align}& \ddot x = \theta g \qquad\qquad\quad \ddot \phi = {{{U_\phi }} / {{I_{xx}}}} \nonumber\\[4pt] & \ddot y = - \phi g \qquad\qquad \ddot \theta = {{{U_\theta }} / {{I_{yy}}}}\nonumber\\[4pt] & \ddot z = {U_z} / m - g\qquad \ddot \psi = {{{U_\psi }} / {{I_{zz}}}} \end{align}

Let ${\ddot z_1} = \ddot z + g$ , define state vectors $x = {[x,y,{z_1},\phi ,\theta ,\psi ,\dot x,\dot y,{\dot z_1},\dot \phi ,\dot \theta ,\dot \psi ]^{\rm{T}}}$ , $y = {[x,y,{z_1},\phi ,\theta ,\psi ]^{\rm{T}}}$ , the system inputs $u = {[{u_1},{u_2},{u_3},{u_4}]^{\rm{T}}}$ . By utilising the state vectors, the linearised model of the quadrotor UAV is obtained as shown in Equation (4).

(4) \begin{equation}\left\{ \begin{array}{l}{\dot x(t) = Ax(t) + Bu(t)}\\[4pt]y(t) = Cx(t)\end{array} \right.\end{equation}

where

\begin{equation*}A = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\[4pt]0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\[4pt]0 & 0 & 0 & 0 & g & 0 & 0 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & { - g} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}} \right]\end{equation*}

\begin{equation*}B = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0\\[4pt]{\dfrac{1}{m}{K_u}} & {\dfrac{1}{m}{K_u}} & {\dfrac{1}{m}{K_u}} & {\dfrac{1}{m}{K_u}}\\[12pt]0 & 0 & {\dfrac{1}{{{I_{xx}}}}{L_d}{K_u}} & { - \dfrac{1}{{{I_{xx}}}}{L_d}{K_u}}\\[8pt]{\dfrac{1}{{{I_{yy}}}}{L_d}{K_u}} & { - \dfrac{1}{{{I_{yy}}}}{L_d}{K_u}} & 0 & 0\\[14pt]{\dfrac{1}{{{I_{zz}}}}{K_y}} & {\dfrac{1}{{{I_{zz}}}}{K_y}} & { - \dfrac{1}{{{I_{zz}}}}{K_y}} & { - \dfrac{1}{{{I_{zz}}}}{K_y}}\end{array}} \right]\end{equation*}

\begin{equation*}C = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\[4pt]0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\end{array}} \right]\end{equation*}

Whether the actuator can fully realise the ability of the controller to control the action $u(t)$ depends largely on its health level. Assuming that the initial execution capability of the ith actuator is ${c_{ia0}}$ and the degradation state is ${d_i}(t)$ , then the degradation rate of the actuator is ${\lambda _i}(t) = {d_i}(t)/{c_{ia0}}$ . The loss matrix of actuator control efficiency is $\Gamma = diag\{ {\lambda _1}(t),{\lambda _2}(t), \cdots ,{\lambda _m}(t)\} $ , $m$ represents the number of actuators, ${\lambda _i}(t) = 0$ indicating that the actuator has not degraded and ${\lambda _i}(t) = 1$ indicating that the actuator has completely failed. Thus, the actual output of the actuator at time t can be expressed as ${u_a}{\rm{(}}t) = {\rm{(}}I{\rm{ - }}\Gamma )u{\rm{(}}t)$ .

The linear model of a quadrotor UAV considering actuator degradation is mathematically represented as Equation (5):

(5) \begin{equation}\left\{ {\begin{array}{l}{\dot x(t) = Ax(t) + B(I - \Gamma )u(t)}\\[4pt]y(t) = Cx(t)\end{array}} \right.\end{equation}

2.2 Construction of autonomous PdM architecture

The autonomous PdM of a UAV refers to the real-time prediction of RUL based on the airframe condition monitoring data, and the autonomous maintenance through the appropriate adjustment of the control strategy as needed. Therefore, to make a quadrotor UAV with multi-actuator degradation have PdM function should first start with the construction of an autonomous PdM architecture,as shown in Fig. 2. The system includes modules such as the airframe, actuators, LF-MPC and PdM. The core functions of the PdM module are twofold: first, TSKF is used to estimate the UAV state and actuators degradation states in real-time when they are degraded, and then determine the UAV failure threshold, predict the RUL of the airframe; the second is to determine whether the UAV needs autonomous maintenance based on RUL and provide maintenance strategies. That is, when RUL does not meet the expected working time limit, LF-MPC is used to enhance the constraints imposed on the control action on the actuator, slow down actuator degradation and extend the effective service time of the quadrotor UAV to achieve autonomous maintenance.

Figure 2. Structure diagram of autonomous maintenance of quadrotor UAV.

Where ${y_d}(t)$ is the expected trajectory, $y(t)$ is the actual trajectory, $e(t)$ is the system deviation, and $u(t)$ is the control effect imposed on the actuator.

3.1 Description of the actuator degradation process

The brushless DC motor is the power source of the quadrotor UAV, its health status directly determines the service life of the UAV. During UAV flight, the execution force of the motor will decrease due to stator coil aging, excessive rotor shaft friction, and magnetic degradation. Considering that the above degradation processes all have nonlinear stationary degradation characteristics, the degradation of each actuator of the UAV can be modeled using the Wiener process, which is described as follows:

(6) \begin{equation}{d^{(m)}}(t){\rm{ = }}{d^{(m)}}(0){\rm{ + }}\int_0^{\rm{t}} {{\mu ^{(m)}}\left( {\tau ;\,{\theta ^{(m)}}} \right)} d\tau + \sigma _B^{(m)}{B^{(m)}}(t)\end{equation}

where $m = 1, \cdots ,4$ denotes the number of actuators. For the mth actuator, ${d^{(m)}}(0)$ represents the initial degradation of the actuator, without loss of generality, let ${d^{(m)}}(0) = 0$ . represents the standard Brownian motion, while ${\mu ^{(m)}}(t;\,{\theta ^{(m)}})$ and $\sigma _B^{(m)}$ represent the drift and diffusion coefficients of the degeneration process, respectively. ${\mu ^{(m)}}(t;\,{\theta ^{(m)}})$ is a nonlinear function of time, and this paper uses ${\mu ^{(m)}}(t;\,{\theta ^{(m)}}){\rm{ = }}{\alpha ^{(m)}}{\beta ^{(m)}}{t^{{\beta ^{(m)}} - 1}}$ to characterise the nonlinear characteristics of the model.

It is a prerequisite and foundation for autonomous PdM to obtain online observation data to estimate the hidden state and predict the RUL. Taking actuator degradation as a type of unknown bias in the system, TSKF can be used to estimate the system state and degradation rate while minimising the impact of measurement noise, thereby achieving more reliable estimation performance.

Discretise Equation (5) to obtain

(7) \begin{equation}\left\{ \begin{array}{l}{{x_{k + 1}} = {A_k}{x_k} + {B_k}u{}_k - {B_k}U{}_k{\gamma _k} + \omega _k^x}\\[4pt]{{y_k} = {C_k}{x_k} + {v_k}}\end{array} \right.\end{equation}

Here, the system noise $\omega _k^x$ and measurement noise ${v_k}$ are uncorrelated Gaussian white noise with zero mean, $Q_k^x$ and ${R_k}$ are the variance matrices of the two, respectively.

\begin{equation*}{\gamma _k} = \left[ {\begin{array}{*{20}{c}}{\lambda _k^{(1)}}\\[4pt]{\lambda _k^{(2)}}\\[4pt] \vdots \\[4pt]{\lambda _k^{(m)}}\end{array}} \right],0 \lt \lambda _k^{(m)} \lt 1;\,m = 1,2, \cdots ,4,\,{U_k} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{u_k^{(1)}} & 0 & \cdots & 0\\[4pt]0 & {u_k^{(2)}} & {} & \vdots \\[4pt] \vdots & {} & \ddots & 0\\[4pt]0 & \cdots & 0 & {u_k^{(m)}}\end{array}} \right],\,{\gamma _{k + 1}} = {\gamma _k} + \omega _k^\gamma \end{equation*}

The bias process propagation noise $\omega _k^\gamma $ is an uncorrelated Gaussian white noise with zero means, and its covariance matrix is represented by $Q_k^\gamma $ .

The augmented discrete linear time-varying state space model with bias is formulated as

(8) \begin{equation}\left\{ \begin{array}{l}{x_{k + 1}} = {A_k}{x_k} + {B_k}u{}_k - {B_k}U{}_k{\gamma _k} + \omega _k^x\\[4pt]{\gamma _{k + 1}} = {\gamma _k} + \omega _k^\gamma \\[4pt]{{y_{k + 1}} = {C_k}{x_{k + 1}} + {v_{k + 1}}}\end{array} \right.\end{equation}

The minimum variance solution for estimating actuator degradation rate and state truth using TSKF is as follows:

Bias estimation

(9) \begin{equation}{\hat \gamma _{k + 1|k}} = {\hat \gamma _{k|k}}\end{equation}

(10) \begin{equation}P_{k + 1|k}^\gamma = P_{k|k}^\gamma + Q_k^\gamma \end{equation}

(11) \begin{equation}{\hat \gamma _{k + 1|k + 1}} = {\hat \gamma _{k + 1|k}} + K_{k + 1}^\gamma \left({\tilde r_{k + 1}} - {H_{k + 1|k}}{\hat \gamma _{k|k}}\right)\end{equation}

(12) \begin{equation}K_{k + 1}^\gamma = P_{k + 1|k}^\gamma H_{k + 1|k}^T{\left({H_{k + 1|k}}P_{k + 1|k}^\gamma H_{k + 1|k}^T + {\tilde S_{k + 1}}\right)^{ - 1}}\end{equation}

(13) \begin{equation}P_{k + 1|k + 1}^\gamma = \left(I - K_{k + 1}^\gamma {H_{k + 1|k}}\right)P_{k + 1|k}^\gamma \end{equation}

Unbiased state estimation

(14) \begin{equation}{\tilde x_{k + 1|k}} = {A_k}{\tilde x_{k|k}} + {B_k}{u_k} + {W_k}{\hat \gamma _{k|k}} - {V_{k + 1|k}}{\hat \gamma _{k|k}}\end{equation}

(15) \begin{equation}\tilde P_{k + 1|k}^x = {A_k}\tilde P_{k|k}^xA_k^T + Q_k^x + {W_k}P_{k|k}^\gamma W_k^T - {V_{k + 1|k}}P_{k + 1|k}^\gamma V_{k + 1|k}^T\end{equation}

(16) \begin{equation}{\tilde x_{k + 1|k + 1}} = {\tilde x_{k + 1|k}} + \tilde K_{k + 1}^x{\tilde r_{k + 1}}\end{equation}

(17) \begin{equation}\tilde K_{k + 1}^x = \tilde P_{k + 1|k}^xC_{k + 1}^T{\tilde S_{k + 1}}^{ - 1}\end{equation}

(18) \begin{equation}\tilde P_{k + 1|k + 1}^x = \left(I - \tilde K_{k + 1}^x{C_{k + 1}}\right)\tilde P_{k + 1|k}^x\end{equation}

Filter residual and its covariance matrix

(19) \begin{equation}{\tilde r_{k + 1}} = {y_{k + 1}} - {C_{k + 1}}{\tilde x_{k + 1|k}}\end{equation}

(20) \begin{equation}{\tilde S_{k + 1}} = {C_{k + 1}}\tilde P_{k + 1|k}^xC_{k + 1}^T + Q_{k + 1}^v\end{equation}

Coupled equations

(21) \begin{equation}{W_k} = {A_k}{V_{k|k}} - {B_k}U{}_k\end{equation}

(22) \begin{equation}{V_{k + 1|k}} = {W_k}P_{k|k}^\gamma {\left(P_{k + 1|k}^\gamma \right)^{ - 1}}\end{equation}

(23) \begin{equation}{H_{k + 1|k}} = {C_{k + 1}}{V_{k + 1|k}}\end{equation}

(24) \begin{equation}{V_{k + 1|k + 1}} = {V_{k + 1|k}} - \tilde K_{k + 1}^x{H_{k + 1|k}}\end{equation}

The state optimal estimate and its estimation error covariance matrix are denoted by

(25) \begin{equation}{\hat x_{k + 1|k + 1}} = {\tilde x_{k + 1|k + 1}} + {V_{k + 1|k + 1}}{\hat \gamma _{k + 1|k + 1}}\end{equation}

(26) \begin{equation}{P_{k + 1|k + 1}} = {\tilde P_{k + 1|k + 1}} + {V_{k + 1|k + 1}}P_{k + 1|k + 1}^\gamma V_{k + 1|k + 1}^T\end{equation}

The real degradation information $d_k^{(m)}$ of each actuator can be accurately obtained based on the actuator degradation rate ${\hat \gamma _k}$ estimated by TSKF, which serves as a fundamental basis for the precise prediction of the RUL for a UAV. According to the definition of degradation rate, $d_k^{(m)}$ can be obtained as shown in Equation (27).

(27) \begin{equation}d_k^{(m)} = {\hat \gamma _k}*c_{a0}^{(m)}\end{equation}

According to the estimated hidden degradation of the actuator, a more accurate distribution of the RUL for a UAV can be obtained by real-time updating.

3.2 Analysis of the degradation process of quadrotor UAV under actuator hidden degradation

Usually, the motors with the same parameters are installed on each rotor of the quadrotor UAV. During flight, degradation occurs due to factors such as coil aging and magnetic degradation or other factors. At the same time, due to various uncertain factors such as air interference and foreign matter, there are individual differences in the degradation process and the coupled superposition effect of common degradation. There are significant differences between multi-actuator degradation and single-actuator degradation on the performance of quadrotor UAVs, resulting in the uncertainty of ${t_n}$ , ${t_{sf}}$ and ${t_{af}}$ . Where ${t_n}$ , ${t_{sf}}$ and ${t_{af}}$ represent the time when the airframe enters degradation from the normal state, the time when the airframe performance does not meet the requirements and the time when the actuator fails.

For the quadrotor UAV, failure means that the airframe performance no longer meets the requirements, and the actuator has not reached its failure threshold at this time. As an important performance index to measure the health status of equipment or system, a HD can reflect the health status of equipment or system more comprehensively. Therefore, the health state of the quadrotor UAV is quantified as a real-time HD to measure the health of the quadrotor UAV at each monitoring time, also the system failure threshold can be defined according to the HD.

Assuming that the desired state of the hovering quadrotor UAV is represented by the sample mean $\mu = {x_d} = {({x_{1d}},{x_{2d}}, \cdots ,{x_{12d}})^T}$ and the actual state estimated by the kth monitoring point is denoted as $\tilde x = {({\tilde x_1},{\tilde x_2}, \cdots ,{\tilde x_{12}})^T}$ , then the Mahalanobis distance can be calculated.

(28) \begin{equation}MD(\tilde x) = \sqrt {{{(\tilde x - \mu )}^T}{\Sigma ^{ - 1}}(\tilde x - \mu )} \end{equation}

where $\Sigma $ is the covariance matrix.

The Mahalanobis distance reflects the degree to which the current state deviates from the desired state. A negative function is used to establish the relationship between Mahalanobis distance $MD(\tilde x)$ and health degree $HD(\tilde x)$

(29) \begin{equation}HD({\tilde{\textbf{x}}}) = e{}^{ - b \cdot MD({\tilde{\textbf{x}}})}\end{equation}

where b represents the shape parameter, with b>0.

Let $H{D_n}$ , $H{D_{sf}}$ denote the health degree of the UAV from the normal working state to the degraded state and the HD when the airframe performance does not meet the requirements, respectively. Considering the full life cycle, it is assumed that multiple actuators have different degradation processes, and curves 1, 2, 3 and 4 in Fig. 3 give one of the degradation combinations. In Fig. 3, ${l_n}$ , ${l_{sf}}$ , and ${l_{af}}$ represent the critical threshold when the quadrotor UAV from normal work to degradation, the degradation threshold when the airframe performance does not meet the requirements, and the comprehensive degradation threshold when actuator failure, respectively. For each actuator, the degree of degradation is different, resulting in ${t_n}$ , ${t_{sf}}$ and ${t_{af}}$ are different. Therefore, to judge the influence of the degradation degree of each actuator on the airframe performance more conveniently, the comprehensive degradation quantity $ds(t)$ (shown in the red curve in Fig. 3) can be constructed by fusing the degradation state of each actuator. In the case of multi-actuator degradation, it becomes more direct to predict the RUL of the airframe with the $HD$ as a constraint.

Figure 3. Schematic diagram of the actuator degradation process based on HD.

The comprehensive degradation of actuators based on the definition of airframe HD can be specifically expressed as follows:

Stage 1, $HD \in (H{D_n},1]$ , the UAV meets the work performance requirements, the integrated degradation $ds(t) \lt {l_n}$ , and the degradation of each actuator has little impact on the overall performance of the UAV.

Stage 2, $HD \in [H{D_{sf}},H{D_n}]$ , the performance of the UAV decreases, and the comprehensive degradation $ds(t) \in [{l_{sf}},{l_n}]$ , indicating that the degradation degree of the actuator is gradually aggravated, which leads to the degradation of the UAV performance, and the RUL may not meet the working deadline requirements.

Stage 3, $HD \in [0,H{D_{{\rm{s}}f}})$ , the performance of the UAV does not meet the requirements, and the comprehensive degradation $ds(t) \gt {l_{sf}}$ , even if each actuator does not fail, but the HD of the quadrotor UAV is less than $H{D_{sf}}$ , and the system fails at this time.

3.3 RUL prediction of quadrotor UAV under hidden degradation of actuators

The comprehensive degradation variable, which combines the degradation state of each actuator, is undoubtedly a healthy indicator of the airframe performance. According to the expected airframe performance index, the system failure threshold of unsatisfactory performance is obtained, and then the RUL of a UAV with multi-actuator degradation can be predicted.

Since the degradation degree of each actuator has different effects on the airframe performance, each degradation quantity should not be treated the same when calculating the comprehensive degradation quantity, and different weights should be given according to the degradation process. The basic idea of the entropy weight method is to determine the objective weight according to the degradation degree of each degradation variable [Reference Xia, Yang, Lin, Xie and Zhang28]. In general, if the information entropy of a certain actuator degradation is smaller, it indicates that the degree of actuator degradation is greater, the more information it provides, the greater role it can play in the comprehensive evaluation, and the greater its weight. And vice versa, the steps for calculating the comprehensive degradation quantity using the entropy weight method are as follows.

1. (1) Data standardisation

   Let $dat{a_{kj}} = ({d^{(m)}})^{\prime},(k = 1, \cdots ,n;\,j = 1, \cdots ,m)$ denote m actuators and n samples, and the min-max normalisation method is used to normalise the degradation data of each actuator

   (30) \begin{equation}{\rho _{kj}} = \frac{{dat{a_{kj}} - {\rm{min}}(dat{a_k})}}{{{\rm{max}}(dat{a_k}) - {\rm{min}}(dat{a_k})}}\end{equation}

2. (2) Information entropy of degradation quantity

   (31) \begin{equation}{E_j} = - ln{(n)^{ - 1}}\sum\limits_{k = 1}^n {{P_{kj}}} ln{P_{kj}}\end{equation}
   where, ${P_{kj}} = {\rho _{kj}}/\sum\limits_{k = 1}^n {{\rho _{kj}}} $ , ${\rho _{kj}}$ represent the jth actuator degradation of the kth sample; if ${P_{kj}} = 0$ , then $\mathop {\lim }\limits_{{P_{kj}} \to \infty } {P_{kj}}{\rm{ln}}{P_{kj}} = 0$ is defined.

3. (3) Degradation quantity weight

   The weight of each degradation is calculated using the information entropy.

   (32) \begin{equation}{W_j} = {{1 - {E_j}} \over {m - \sum\limits_{j = 1}^m {{E_j}} }}{\rm{,}}\,\,\,\,{\rm{j}} = 1,2, \cdots {\rm{4}}\end{equation}

4. (4) Comprehensive degradation

   (33) \begin{equation}d{s_k} = \sum\limits_{j = 1}^m {{W_j}{\rm{dat}}{{\rm{a}}_{kj}}} \end{equation}

The original degradation of the BLDC motor is random, so the comprehensive degradation quantity obtained based on the entropy weight method can be modeled by the Wiener process in (6). The increment $\Delta d{s_k} = d{s_{k + 1}} - d{s_k}$ of the amount of degradation between two consecutive detection instants follows a normal distribution, $\Delta d{s_k}\sim N(\alpha (t_{k + 1}^\beta - t_k^\beta ),\sigma _B^2\Delta {t_k})$ . Let $DS({t_k}) = d{s_k}$ denote the integrated degradation quantity at the current time ${t_k}$ , and the integrated degradation state vector $D{S_{1:k}} = (d{s_1},d{s_2}, \cdots ,d{s_k})$ .

Thus, for the unknown parameter $\theta = (\alpha ,\beta ,{\sigma _B})$ of the model, the likelihood function of the unknown parameter $\theta $ in the degraded data $D{S_{1:k}}$ is given by.

(34) \begin{equation}L({\rm{\theta }}|\Delta d{s_i}) = \prod\nolimits_{i = 1}^k {\frac{1}{{\sqrt {2\pi (\sigma _B^2\Delta {t_i})} }}exp\left[ { - \frac{{{{(\Delta d{s_{\mathop{\rm i}\nolimits} } - \alpha (t_{i + 1}^\beta - t_i^\beta ))}^2}}}{{2(\sigma _B^2\Delta {t_i})}}} \right]} \end{equation}

Taking the logarithm of both sides of Equation (34) yields:

(35) \begin{equation}\ln L{\rm{(}}\theta |\Delta d{s_i}) = \sum\nolimits_{i = 1}^k [ ln\frac{1}{{\sqrt {2\pi } }} + ln\frac{1}{{\sqrt {\sigma _B^2\Delta {t_i}} }} - \frac{{{{(\Delta d{s_{\mathop{\rm i}\nolimits} } - \alpha (t_{i + 1}^\beta - t_i^\beta ))}^2}}}{{2(\sigma _B^2\Delta {t_i})}}{\rm{]}}\end{equation}

The logarithmic likelihood function (35) was maximised using the ‘fminsearch’ function in MATLAB simulation software, resulting in obtaining the maximum likelihood estimate $\hat \theta $ for the unknown parameter $\theta = (\alpha ,\beta ,{\sigma _B})$ .

Assuming that the comprehensive failure threshold of multi-actuator degradation is ${l_{af}}$ , and according to the definition of the first reach threshold, its corresponding life distribution is as follows [Reference Si, Hu, Zhang, He and Zhou29]:

(36) \begin{equation}{{\mathop{f}\nolimits} _T}{\rm{(}}t|\hat \theta ) \approx \frac{{{l_{af}} - \hat \alpha {t^{\hat \beta }}{\rm{(}}1 - \hat \beta )}}{{{{\hat \sigma }_B}\sqrt {2\pi {t^3}} }}\exp \left\{ {\frac{{ - {{{\rm{(}}{l_{af}} - \hat \alpha {t^{\hat \beta }})}^2}}}{{2\hat \sigma _B^2t}}} \right\}\end{equation}

$H{D_{sf}}$ is selected as the index to determine the failure threshold of the quadrotor UAV, and the failure threshold ${l_{sf}}$ is defined as the corresponding multi-actuator comprehensive degradation when the health degree $HD(\tilde x)$ is less than $H{D_{sf}}$ , then:

(37) \begin{equation}{l_{sf}} = \left\{ {\left. {ds(t)|HD(\tilde x) \lt H{D_{sf}}} \right\}} \right.\end{equation}

The life of the quadrotor UAV with hidden degradation of multi-actuators refers to the life of the airframe when it works normally to the minimum HD requirement, so the life of the UAV can be defined as

(38) \begin{equation}{t_{sf}} = inf {\rm{\{ }}t:ds{\rm{(}}t) \gt {l_{sf}}|ds{\rm{(}}0) \lt {l_{sf}}{\rm{\} }}\;\end{equation}

where inf is the infimum operator.

The RUL of the system is expressed as follows.

(39) \begin{equation}RUL{\rm{(}}t) = inf{\rm{ \{ }}\;{t_{sf}} > t,\,\,\,\,ds{\rm{(}}t) > {l_{sf}}\;{\rm{\} }}\; - t\end{equation}

Let the comprehensive degradation variable be $d{s_k}$ and ${l_k}$ be the RUL of the UAV predicted at the current time ${t_k}$ . Referring to the literature [Reference Si, Hu, Zhang, He and Zhou29] and combining with the comprehensive degradation threshold determined by Equation (37), the PDF of the RUL for a quadrotor UAV is:

(40) \begin{equation}{f_L}{\rm{(}}{t_{lk}}|\theta ) \approx \frac{1}{{\sqrt {2\pi \hat \sigma _B^2{t_{lk}}^3} }}{\rm{(}}d{s_s} - \hat \alpha {\rm{(}}\eta ({t_{lk}}) - \hat \beta {t_{lk}}{({t_{lk}} + {t_k})^{\hat \beta - 1}}{\rm{))}}exp\left\{ { - \frac{{{{{\rm{(}}d{s_s} - \hat \alpha \eta ({t_{lk}}))}^2}}}{{2\hat \sigma _B^2{t_{lk}}}}} \right\}\end{equation}

where $\eta ({t_{lk}}){\rm{ = }}{({t_k} + {t_{lk}})^{\hat \beta }} - t_k^{\hat \beta }$ , $d{s_s} = {l_{sf}} - d{s_k}$ .

According to the analytical expression, the real-time prediction of the RUL of a quadrotor UAV can be quickly obtained.

4.1 Life extension mechanism based on LF-MPC

To implement extended life control using MPC, the first step is to discretise model (4) based on the basic steps of MPC model prediction, feedback correction and rolling optimisation

(41) \begin{equation}\left\{ {\begin{array}{l}{{x_m}(k + 1) = {A_m}{x_m}(k) + {B_m}u(k)}\\[4pt]{y_m}(k) = {C_m}{x_m}(k)\end{array}} \right.\end{equation}

where $A_m = A^{\ast} \Delta t + I, B_m = B^{\ast} \Delta t$ .

The original state space model (41) is augmented to obtain $\Delta u$ as a manipulated variable for the future system output so that the feedback correction and rolling optimisation can be performed simultaneously.

Let $\Delta {x_m}(k) = {x_m}(k) - {x_m}(k - 1)$ , $\Delta u(k) = u(k) - u(k - 1)$ .

(42) \begin{equation}\left[ {\begin{array}{c}{\Delta {x_m}(k + 1)}\\[4pt]{{y_m}(k + 1)}\end{array}} \right] = \left[ {\begin{array}{c@{\quad}c}{{A_m}} & {{0_{12 \times 6}}}\\[4pt]{{C_m}{A_m}} & {{I_{6 \times 6}}}\end{array}} \right]\left[ {\begin{array}{c}{\Delta {x_m}(k)}\\[4pt]{{y_m}(k)}\end{array}} \right] + \left[ {\begin{array}{c}{{B_m}}\\[4pt]{{C_m}{B_m}}\end{array}} \right]\Delta u(k)\end{equation}

Make ${x_u}(k) = {[\Delta {x_m}{(k)^T},{y_m}{(k)^T}]^T}$ , Equation (41) be expressed as

(43) \begin{equation}\left\{ {\begin{array}{l}{{x_u}(k + 1) = {A_u}{x_u}(k) + {B_u}\Delta u(k)}\\[4pt]{{y_m}(k){\rm{ = }}{C_u}{x_u}(k)}\end{array}} \right.\end{equation}

where ${A_u} = \left[ {\begin{array}{c@{\quad}c}{{A_m}} & {{0_{12 \times 6}}}\\[4pt]{{C_m}{A_m}} & {{I_{6 \times 6}}}\end{array}} \right]$ , ${B_u} = \left[ {\begin{array}{*{20}{c}}{{B_m}}\\[4pt]{{C_m}{B_m}}\end{array}} \right]$ , ${C_u} = \left[{0_{6 \times 12}}{\begin{array}{*{20}{c}}{} I\end{array}_{6 \times 6}}\right]$ .

Assuming that the current time is $k$ , ${N_c}$ is the control step and ${N_p}$ is the prediction step, a recursive Equation (44) can be obtained.

(44) \begin{align}Y & = F{x_u}(k) + \Phi \Delta U\\[5pt]F & = \left[\begin{array}{c}C_u A_u\\[3pt]C_u A_u^2\\[3pt]\vdots\\[3pt]C_u A_u^{N_p}\end{array}\right], \Phi = \left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}C_u B_u & 0 & \cdots & 0\\[3pt]C_u A_u B_u & C_u B_u & \cdots & 0\\[1pt]\vdots & \vdots & \vdots &\vdots\\[3pt]C_u A_u ^{N_p -1}B_u & C_u A_u ^{N_p -2}B_u & \cdots & C_u A_u ^{N_p- N_c}B_u\end{array}\right]\nonumber\end{align}

Define the prediction output $Y = [{y_m}(k + 1),{y_m}(k + 2), \cdots ,{y_m}(k + {N_p})]$ , control increment $\Delta U(k) = {[\Delta u(k){\rm{,}}\Delta u(k + 1){\rm{ ,}} \cdots {\rm{ ,}}\Delta u(k + {N_c} - 1)]^T}$ .

Set ${Y_d} = [{y_d}(k + 1),{y_d}(k + 2), \cdots ,{y_d}(k + {N_p})]$ be the system reference input at time $k$ . To make the system output as close as possible to the input, for MPC control, the cost function is defined as follows.

(45) \begin{align}J & = ({Y_d} - Y{)^T}Q({Y_d} - Y) + \Delta {U^T}R\Delta U\nonumber\\[4pt]& = \sum\limits_{i = 1}^{{N_p}} {{{({y_d}(k + i) - {y_m}(k + i))}^T}{Q_i}({y_d}(k + i) - {y_m}(k + i))} {\rm{ + }}\sum\limits_{j = 0}^{{N_c} - 1} {{{(\Delta u(k + j))}^T}{R_j}(\Delta u(k + j))} \end{align}

In the above equation, the first term reflects the ability of the UAV to follow the reference trajectory, and the second term reflects the requirement for the smooth change of the control quantities. $Q$ and $R$ are the block diagonal matrices of ${Q_i}$ and ${R_i}$ , respectively; where, ${Q_i} = diag{\{ {q_0}, \cdots ,{q_0}\} _{n \times n}}$ and ${R_i} = diag{\{ {r_0}, \cdots ,{r_0}\} _{m \times m}}$ are the error weight and control quantity weight, ${q_0},{r_0}$ are the initial values of MPC. To adjust the parameters conveniently, the diagonal elements of matrices $Q$ and $R$ are selected as the same values here.

By seeking the extreme value of Equation (43), the optimal control increment at the current moment can be obtained.

(46) \begin{equation}\Delta U(k) = {\left({\Phi ^T}Q\Phi + R\right)^{ - 1}}\left({\Phi ^T}Q{Y_r} - {\Phi ^T}QFx(k)\right)\end{equation}

Based on the idea of rolling optimisation, taking the first element of the sequence as the current control, then

(47) \begin{equation}\Delta u(k) = [10 \cdots 0]\Delta U(k)\end{equation}

In MPC, there are many parameters in each optimisation prediction window of the control sequence, which reduces the solution efficiency. The multi-objective control algorithm based on LF reduces the calculation amount by approaching the finite time domain control trajectory with fewer optimisation parameters, which can save more time for the implementation of autonomous PdM on a UAV.

At time k, the elements within the control trajectory $\Delta U$ are represented in the following impulse response form using discrete $\delta $ -functions combined with $\Delta U$ .

(48) \begin{equation}\Delta u(k + i) = [\delta (i)\delta (i - 1) \cdots \delta (i - {N_c} + 1)]\Delta U(k)\end{equation}

where

\begin{equation*}\left\{\begin{array}{l}\delta(i) = diag([\delta_1(i)\, \delta_2(i) \cdots \delta_m(i)] = diag \underbrace{([1\,\,1\cdots 1])}_m), i = 0\\[6pt]\delta(i) = diag \underbrace{([1\,\,1\cdots 1])}_m), i \neq 0\end{array}\right.\end{equation*}

According to the definition of the LF, $\Delta u(k + i)$ in Equation (48) is denoted by ${L^T}(i)\eta $

(49) \begin{equation}\Delta u(k + i) = {L^T}(i)\eta \end{equation}

where $\eta = [{\eta _1}\quad {\eta _2}\quad \cdots\quad {\eta _m}]{}^T$ is the Laguerre coefficient and $L^T(i) = diag(L_1(i)^T,$ $L_2(i)^T,\quad \cdots ,\quad L_m(i)^T),\quad i \in [1,{N_c}]$ denotes the Laguerre function.

Let ${L_q}(i) = {[{l_{1q}}(i)\quad {l_{2q}}(i)\quad \cdots\quad {l_{Nq}}(i)]^T}(q \in [1,m])$ , where $Nq$ is the number of approximation factors of the qth actuator. The relationship of LF at adjacent times is obtained by the iterative operation as follows.

(50) \begin{equation}{L_q}(i + 1) = {A_{ql}}{L_q}(i)\end{equation}

where matrix ${A_{ql}}$ is a function of parameters ${\alpha _q}$ and ${\beta _q} = 1 - {\alpha _q}$ , and the initial condition is

(51) \begin{equation}{L_q}(0) = \sqrt {{\beta _q}} \left[1\quad {\rm{ - }}{\alpha _q}\quad {\alpha _q}^2\quad \cdots\quad {\rm{ - }}{\alpha _q}^{{N_q} - 1}\right]\end{equation}

(52) \begin{equation}{A_{ql}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}{{\alpha _q}}& 0& 0& \cdots& 0\\[4pt]{{\beta _q}}& {{\alpha _q}}& 0& \cdots& 0\\[4pt]{ - {\alpha _q}{\beta _q}}& {{\beta _q}} &{{\alpha _q}} & \cdots &0\\[4pt] \vdots & \vdots & \vdots & \ddots & \vdots \\[4pt]{{\rm{ - }}{\alpha _q}^{{N_q} - 2}{\beta _q}}& {{\rm{ - }}{\alpha _q}^{{N_q} - 3}{\beta _q}}& {{\rm{ - }}{\alpha _q}^{{N_q} - 4}{\beta _q}} & \cdots &{{\alpha _q}}\end{array}} \right]\end{equation}

Substituting (49) into (43) gives the state variables and output variables of the system at time i after time k as follows

(53) \begin{equation}\left\{ {\begin{array}{l}{{x_u}(k + i|k) = {A_u}^i{x_u}(k) + \sum\limits_{p = 0}^{i - 1} {{A_u}^{i - p - 1}} {B_u}L{{(p)}^T}\eta = {A_u}^i{x_u}(k) + \phi {{(i)}^T}\eta }\\[9pt]{{y_m}(k + i|k){\rm{ = }}{C_u}{x_u}(k + i|k){\rm{ = }}{C_u}{A_u}^i{x_u}(k) + {C_u}\phi {{(i)}^T}\eta }\end{array}} \right.\end{equation}

where $\phi {(i)^T} = \sum\limits_{p = 0}^{i - 1} {{A_u}^{i - p - 1}} {B_u}L{(p)^T}{\rm{, }}(i \in [1,{N_p}])$ .

By replacing $\Delta u(k + i)$ with ${L^T}(i)\eta $ , the parameter $\eta $ becomes the only optimised parameter vector to replace $\Delta U$ . Due to the smaller dimension of parameter $\eta $ than $\Delta U$ , the calculation amount is reduced and the speed is improved, which is the advantage of LF-MPC over traditional MPC in life extension.

According to the LF definition, substituting (49) into (45) yields.

(54) \begin{equation}J(k) = \sum\limits_{i = 1}^{{N_p}} {{{({y_d}(k + i) - {y_m}(k + i))}^T}} {Q_i}({y_d}(k + i) - {y_m}(k + i)){\rm{ + }}{\eta ^T}R\eta \end{equation}

Substituting Equation (53) into (54), retaining only the terms associated with $\eta $ , and obtain

(55) \begin{equation}J(k) = {\eta ^T}\left(\sum\limits_{i = 1}^{{N_p}} \phi (i)\Omega \phi (i) + R\right)\eta + 2{\eta ^T}\left(\sum\limits_{i = 1}^{{N_p}} {\phi (i)\Omega {A_u}^i{x_u}(k)} \right){\rm{ - }}2{\eta ^T}\left(\sum\limits_{i = 1}^{{N_p}} {\phi (i){C_u}^T{Q_i}{y_d}(k + i|k)} \right)\end{equation}

Let $\frac{{\partial J}}{{\partial \eta }}{\rm{ = }}0$ , then the parameter vector $\eta $ that optimises the performance index can be obtained as

(56) \begin{equation}{\eta ^ * } = \left(\sum\limits_{i = 1}^{{N_p}} \phi (i)\Omega \phi {{(i)}^T} + R\right)^{ - 1} \times \left(\sum\limits_{i = 1}^{{N_p}} \phi (i)\Xi ({y_d}(k + i) - {C_u}{A_u}^i{x_u}(k))\right) \end{equation}

where $\Omega = {C_u}^T{Q_i}{C_u}$ , $\Xi = {C_u}^T{Q_i}$ .

According to the life extension mechanism, by decreasing the element value of $Q$ and increasing the element value of $R$ , it is possible to relax the constraint of error performance indicators, reduce actuator execution pressure, and delay degradation, thereby achieving the goal of prolonging the RUL for the quadrotor UAV.

Assuming that the expected safe working time of a quadrotor UAV is ${t_d}$ , in order to achieve the expected working time limit ${t_d}$ , the shortest expected life extension time ${t_{ex}}$ should be calculated based on the RUL prediction result ${t_{lk}}$ at time ${t_k}$ . When ${t_{ex}} \gt 0$ , the system needs to be maintained autonomously.

(57) \begin{equation}{t_{ex}} = {t_d} - ({t_{lk}} + {t_k}) = {t_d} - {t_{sf}}\end{equation}

4.2 Adaptive autonomous PdM strategy of LF-MPC based on HD

From the above analysis, it can be seen that there is a close relationship between actuator degradation and system performance. When the degradation is small, due to the effect of feedback control, the degradation has a relatively small impact on system performance, while when the degradation increases, the system performance will deteriorate. Due to the definition of the Mahalanobis distance HD, which comprehensively considers the impact of actuator degradation on various states of the system, can provide a more comprehensive evaluation of the impact of actuator degradation on the system. Therefore, when the quadrotor UAV does not meet the requirements of a working time limit, autonomous maintenance can adjust the weight matrices ${Q_i}$ and $R$ in real-time based on HD, to achieve a balance between steady-state performance and control capability of the UAV.

The UAV initially has no degradation, with $HD = 1$ . As the hidden degradation of the actuator increases, the airframe state deviates from the expected value, and the health level decreases. The more severe the degradation, the smaller the health level. Combining the HD evaluation results with maintenance requirements, degradation has a relatively small impact on the UAV during $HD \gt H{D_n}$ and does not require correction; when $H{D_{sf}} \le HD \le H{D_{\rm{n}}}$ does not meet the working time limit ${t_d}$ requirements, adjust the weight matrix. The greater the change in $HD$ , the greater the maintenance requirement. Increase the ${Q_i}$ and $R$ adjustment scale, and the smaller the change in $HD$ , the smaller the maintenance requirement, and fine-tune ${Q_i}$ and $R$ ; however, when $HD \lt H{D_{sf}}$ , the correction of the weight matrices of ${Q_i}$ and $R$ should be stopped. In summary, the adjustment trend of the weight matrix integrating HD and maintenance demand in Fig. 4 can be obtained.

Figure 4. Schematic diagram of weight matrix $\Delta q$ and $\Delta r$ adjustment based on HD.

It can be seen that the changes of $\Delta q$ and $\Delta r$ with $HD$ roughly follow an exponential law. Therefore, the matrices ${Q_i}$ and $R$ , and correction quantities $\Delta q$ and $\Delta r$ based on HD can be expressed as:

(58) \begin{equation}\left\{ {\begin{array}{l}{q(k + 1) = q(k) - \Delta q,\Delta q = \dfrac{1}{{{K_1} + {e^{HD}}}}}\\[10pt]{r(k + 1) = r(k) + \Delta r,\Delta r = \dfrac{1}{{{K_2} + {e^{HD}}}}}\end{array}} \right.\end{equation}

where, $q(t) \geqslant {q_{min}}$ , $r(t) \le {r_{max}}$ , the parameters can change the range of adjustment parameters $\Delta q$ and $\Delta r$ .

4.3 Autonomous PdM process of quadrotor UAV based on HD

For the quadrotor UAV with multi-actuator degradation, the TSKF method is used to estimate the degradation of each actuator and the UAV state in real-time, the entropy weight method is used to obtain the comprehensive variables and the system HD is calculated in real-time. Considering the system performance and deadline constraints, the online autonomous maintenance strategy of a quadrotor UAV is implemented as in Fig. 5.

5.1 Experimental description

To verify the effectiveness of the proposed approach in the article, the object parameters were selected for simulation experiments as shown in Table 1. The LF-MPC control parameters are taken as follows: ${N_c} = 30$ , ${N_p} = 5$ . The Laguerre polynomial parameter is taken as $\alpha = [0.7,0.7,0.7,0.7]$ , ${N_q} = [7,7,7,7]$ . The initial values of the error weight and control quantity weight matrix are taken as ${q_0} = 2$ and ${r_0} = 0.01$ , respectively, and the constraints ${q_{min}} = 1.5$ and ${r_{max}} = 1.0$ are satisfied, and the sampling period is $\Delta t = 0.1$ .

Table 1. Physical parameters of quadrotor UAV

Assuming that the HDs of the quadrotor UAV is $H{D_n} = 0.9$ and $H{D_{sf}} = 0.7$ respectively when the airframe performance starts to degrade and does not meet the requirements, and the expected working time limit is ${t_d} = 500min$ .

5.2 RUL prediction of a quadrotor UAV

Set the initial value of the filter to

\begin{equation*}\tilde x = \hat x = {0_{12 \times 1}},\,\hat \gamma = {0_{4 \times 1}},\,\,{P_0}^x = {\tilde P_0}^x = {E_{12}}\end{equation*}

\begin{equation*}{P_0}^\gamma = {\tilde P_0}^\gamma = {E_4},\,Q_{}^\gamma = {10^{{\rm{ - }}2}} \times {I_{4 \times 4}}\end{equation*}

\begin{align*}Q_{}^x = \left[ {\begin{array}{c@{\quad}c}{{{10}^{{\rm{ - }}10}} \times {I_{6 \times 6}}} & {{0_{6 \times 6}}}\\[4pt]{{0_{6 \times 6}}} &{{{10}^{{\rm{ - }}12}} \times {I_{6 \times 6}}}\end{array}} \right],\,\,{Q^v} = \left[ {\begin{array}{c@{\quad}c}{{{10}^{{\rm{ - }}12}} \times {I_{6 \times 6}}} & {{0_{6 \times 6}}}\\[4pt]{{0_{6 \times 6}}} &{{{10}^{{\rm{ - }}14}} \times {I_{6 \times 6}}}\end{array}} \right]\end{align*}

Considering that the four actuators of the quadrotor UAV undergo varying degrees of degradation, the TSKF algorithm is applied to obtain the estimated value of the degradation process shown in Fig. 6. It can be seen that TSKF can effectively estimate the true values and reflect the actual degradation process.

Figure 5. Flow chart of autonomous PdM algorithm for quadrotor UAV based on HD.

Figure 6. True value and estimated value of actuator degradation.

The curve of the comprehensive degradation quantity obtained by using the aforementioned entropy weight method is shown in Fig. 7. From the graph, it can be seen that the comprehensive degradation $ds(t)$ has a high consistency with the degradation trend of the severely degraded actuator 1, so the comprehensive degradation also follows the Wiener process, and the parameter estimation uses the ‘fminsearch’ function, which is obtained

Figure 7. Degradation curves for both comprehensive and individual actuator.

$\hat \alpha {\rm{ = }}0.0021,\,\hat \beta {\rm{ = }}1.3373,\,\,{\hat \sigma _{\rm{B}}}{\rm{ = }}0.0098.$

The red line in Fig. 8 shows the position and attitude changes of the quadrotor UAV during the degradation process of the four actuators before implementing the maintenance strategy. As an important indicator of autonomous maintenance, an airframe HD is an important basis for determining maintenance timing and maintenance of a quadrotor UAV.

Figure 8. Response curve of quadrotor UAV with actuator degradation.

The red line in Fig. 9 shows the change curve of airframe HD before maintenance.

Figure 9. Change curve of HDs.

For the health failure threshold $H{D_{sf}} = 0.7$ , according to (37), the comprehensive failure threshold ${l_{sf}} = 0.3587$ is obtained, and the RUL distribution of the quadrotor UAV is calculated using (40). Then the RUL distribution is shown in Fig. 10. When the quadrotor UAV exceeds the desired HD constraint, the life of the UAV is ${t_{sf}} = 399\,min$ , which does not meet the desired working time requirements. It can be seen that the airframe performance is constrained based on the HD and the comprehensive degradation perception quantity is combined to make the airframe RUL prediction more direct, convenient and accurate.

Figure 10. Distribution of RUL of quadrotor UAV.

5.3 Autonomous PdM effect of quadrotor UAV

According to the above autonomous maintenance strategy of the quadrotor UAV, when $H{D_n} = 0.9$ , the UAV enters a degraded state. If the minimum expected life extension time ${t_{{\rm{ex}}}} \gt 0$ , it is considered to adaptively adjust $\Delta q$ and $\Delta r$ based on (58), and carry out autonomous maintenance by reducing the $Q$ matrix and increasing the weight value of each diagonal element in the $R$ matrix.

Therefore, maintenance begins after the airframe enters a degraded state. Take ${K_1} = 120$ and ${K_2} = 100$ to obtain the matrix $Q$ and $R$ correction values $\Delta q$ and $\Delta r$ adjustment curves in Fig. 11, as well as the weight matrix $Q$ and $R$ change curves after adjustment in Fig. 12.

Figure 11. Adjustment curve of $\Delta q$ and $\Delta r$ .

Figure 12. Adjustment curve of Q and R for life extension control.

Based on this maintenance strategy, the position, attitude and HD changes of the UAV after maintenance are obtained as shown in Figs 8 and 9 with blue curves.

Figure 13 shows the variation curve of the control voltage u of the actuator before and after autonomous maintenance. It can be seen that the output pressure of the actuator is effectively relieved, the degradation process is delayed, and the HD of the UAV is also improved after maintenance, which can provide more time for the normal work of the UAV.

Figure 13. Control quantity curve before and after actuator maintenance.

After adopting the autonomous PdM method proposed in this paper, the RUL distribution of the UAV is shown in Fig. 14, and the airframe life is extended to ${t_{sf}} = 551\,min$ , which meets the requirement of the HD and working time of the UAV at the same time. This further verifies the effectiveness of the method of adaptively adjusting the weight matrix based on the HD for the autonomous maintenance of the UAV.

To verify the superiority of the adaptive maintenance strategy based on HD in this paper, it is compared with the fixed factor weight correction method and the MPC HD maintenance method. The adjustment factors of $Q$ and $R$ are ${\delta _1} = 0.008$ and ${\delta _2} = 0.009$ . Figure 9 shows the change curves of the HD of the three methods, and the advantages of the proposed method can be seen.

Table 2 shows the failure thresholds and airframe life after no maintenance and different maintenance strategies. It can be seen that there are two advantages to using the LF-MPC method. Firstly, compared to MPC control, under the same control time domain conditions, the LF-MPC calculation efficiency is improved and the control speed is increased. Therefore, the airframe life before the extension of life is increased by 12 min compared to MPC, which improves the control effect; secondly, from the perspective of maintenance strategy, MPC HD maintenance is better than fixed value maintenance, and LF-MPC HD maintenance is better than MPC health maintenance. The reason is also that the adaptive maintenance strategy based on LF-MPC HD is far less than the control time domain after the adoption of LF, resulting in a relatively reduced amount of computation in the solution of quadratic programming, improved real-time performance and more targeted maintenance.

Table 2. Failure threshold and airframe life under different maintenance strategies

Figure 14. Distribution of the RUL of the quadrotor UAV after maintenance.