Convex Fault Diagnosis of a Three-Degree-of-Freedom Mechanical Crane

Abstract

This paper presents a fault detection and estimation method based on a proportional-integral observer applied to a three-degree-of-freedom mechanical crane. Faults are common in this system and can provoke oscillations that generate a loss of performance and stability. A convex linear parameter varying approach is proposed to stabilize the crane and detect and isolate actuator faults to guarantee the crane’s performance. The linear matrix inequalities obtained from candidate Lyapunov functions give sufficient conditions to guarantee the fault estimation method. Finally, numerical simulations are proposed to illustrate the method’s performance and applicability.

1. Introduction

Cranes are essential in modern industry and are used in various applications, including construction, mining, and heavy-load transport. These machines are found in factories, warehouses, construction sites, and shipyards worldwide. However, due to their high complexity, it is not uncommon for cranes to experience faults, leading to severe consequences, such as accidents and material damage. Developing effective strategies to minimize load oscillations, address sensor noise, and reduce maintenance costs is crucial to mitigate these risks [1]. To address this problem, different methodologies for fault detection have been studied. These techniques span from data-based to model-based approaches. Data-based approaches consider neuronal networks that learn from experience to reduce the load oscillation [2,3], fuzzy logic with if-then rules [4,5], machine learning to estimate the payload-mass lifted [6], machine learning under conditions of strong coastal winds [7], and the genetic algorithm optimization model [8], among others. On the other hand, model-based methods consider dynamical models obtained from physical principles [9,10], where the use of state observers has been fundamental in fault diagnosis and fault-tolerant control schemes [11,12,13]. It is essential to consider that a crane is a nonlinear and underactuated system, which means that it has more degrees of freedom (DOF) than actuators, making the design of model-based fault diagnosis methods more complex. Generally, the analysis is based on linear model approximations that reduce complexity but also reduce representativity and, in consequence, robustness [14].

Crane mathematical models are usually represented by nonlinear differential equations based on the pendulum dynamics [15] and recently also focused on tower cranes [16]. Many reports consider the load oscillation based on these models, for example, by considering a back-stepping approach [17,18], nonlinear optimal control for sway control [19], second-order sliding-mode control [20,21], and disturbance-observer-based nonlinear control [22,23], among others. However, the design of these controllers is challenging and cannot be generalized due to their complexity. In this scenario, linear parameter varying (LPV) models have recently been considered to represent complex dynamics of nonlinear systems. LPV systems are composed of a set of linear time-invariant models that are interpolated by a set of scheduling functions; as a result, LPV models represent nonlinear systems with high accuracy and lower complexity than nonlinear models. In addition, some techniques initially designed for linear systems, such as linear matrix inequalities, can be extended to the LPV case, which reduces the conservatism of the controller solutions and increases its applicability [24].

Few works have been reported considering LPV approaches for studying the performance of cranes. The work by González et al. [25] uses parameter-dependent Lyapunov functions to improve the controller and reduce oscillations. Similarly, Aktas et al. [26] introduce a method that combines tuning the derivative, integral, and proportional parameters with LPV modeling. Furthermore, few researchers have proposed various methods for fault detection in this type of crane, such as the work of Chen and Saif [27]; the method is based on a new input/output relationship derived from the considered nonlinear systems and robust high-order sliding mode differentiators. Meanwhile, Zheng and Zhao [28] utilize historical fault data to build a comprehensive spreader fault tree with three layers of fault phenomena, classification, and causes. These approaches demonstrate the effectiveness of different techniques in fault detection for 3DOF cranes and provide valuable insights for developing reliable and efficient crane systems. The study of using state observers to detect crane faults is addressed by Sjöberg [29]. The work presents a linearized observer that can generate residuals, which can then be used to estimate potential faults in the crane. However, the mathematical model of the 3DOF crane is composed of a nonlinear set of differential equations, so designing a fault diagnosis system based on a linear model limits the applicability of the design to a reduced operating region. Additionally, in these works [27,28,29], the solutions are based on either historical data or the generation of residuals where the violation of a threshold determines if there is a failure, and the best that can be achieved is the isolation of the fault. Unfortunately, these approaches do not provide a means to estimate the magnitude of the fault.

This paper introduces a methodology focused on fault diagnosis for nonlinear systems using qLPV-based state observers. The observer used adopts a proportional integral structure and has as its main characteristic that its representation covers both the states of the system and its possible faults. The application of this methodology is exemplified in a 3DOF crane system, which provides a reliable solution for the optimal operation of this equipment. This choice is significant because cranes are complex and critical systems in various industrial fields. The proposed approach presents an exciting alternative for detecting failures in this type of system and providing a reliable and accurate way to guarantee the correct operation of this type of equipment, avoiding faults that could be expensive or dangerous. The highlight of this methodology lies in its ability to identify the presence of faults and provide the underlying dynamics of these faults. This is valuable to significantly prevent further damage to the crane and enable safety operations with the aid of a fault-tolerant control algorithm (not discussed here). The contribution of this work is the synergy achieved by combining a qLPV representation, which allows variations in system behavior to be captured as a function of multiple parameters, with the design of a PI observer. The designed observer plays a central role in the fault diagnosis scheme by enabling accurate estimation of system states and faults present.

The structure of this document is as follows. Section 2 addresses the description of the nonlinear model of the 3DOF Crane and its qLPV representation. The methodology used to stabilize the system is detailed in Section 3. Likewise, Section 4 presents the methodology for designing a PI observer, where the gains of this observer are calculated through LMIs. This approach allows us not only to estimate the state of the system but also to identify its faults. The implementation of the fault diagnosis scheme in the nonlinear model and the results obtained are described in Section 5. Finally, conclusions derived from this work are presented in Section 6.

2. Mathematical Model

Let us consider the rigid-body diagram of a 3DOF crane as shown in Figure 1. The system comprises a cart (trolley) for horizontal movement and a hoisting/lowering mechanism of the payload of mass m supported by a rope of length $l\left(t\right)$. In the particular case of the illustrated crane, three different movements represent its degree of freedom (DOF). The first DOF is the motor that generates the force $F_x\left(t\right)$ that moves the cart along the x-axis. The motor that applies the necessary force $F_l\left(t\right)$ to ascend/descend the load represents the second DOF. Finally, the third DOF is represented by the angle $\theta \left(t\right)$ of the load oscillations on the x-axis. It is important to note that because there are more degrees of freedom than controller actuators, the system is classified as an underactuated system, which makes controlling the load oscillation challenging. Furthermore, in this work, we assume that the actuators can be affected by faults that risk the crane’s safe operation. Therefore, the design of control methods for overhead crane systems represents a significant challenge from theoretical and practical points of view. The following ordinary differential equations give the nonlinear model [30]:

$$M\left(q\right)\ddot{q}+D\dot{q}+C\left(q,\dot{q}\right)\dot{q}+G\left(q\right)=F,$$

(1)

where $M\left(q\right)$ is a symmetric mass matrix, D is the damping matrix, $C(q,\dot{q})$ is the centrifugal forces matrix, $G\left(q\right)$ is the vector of gravitational forces, q is a generalized coordinate vector and F is vector of control inputs, with

$$\begin{array}{cc}\hfill & M\left(q\right)=\left[\begin{array}{ccc}{M}_x+m& m{S}_{\theta }& ml\left(t\right)C_{\theta }\\ m{S}_{\theta }& M_l+m& 0\\ ml\left(t\right)C_{\theta }& 0& ml{\left(t\right)}^2\end{array}\right];D=\left[\begin{array}{ccc}{D}_x& 0& 0\\ 0& D_l& 0\\ 0& 0& 0\end{array}\right];\hfill \\ \hfill & C\left(q,\dot{q}\right)=\left[\begin{array}{ccc}0& 2m\dot{\theta }\left(t\right)C_{\theta }& -ml\left(t\right)\dot{\theta }\left(t\right)S_{\theta }\\ 0& 0& -ml\left(t\right)\dot{\theta }\left(t\right)\\ 0& 2ml\left(t\right)\dot{\theta }\left(t\right)& 0\end{array}\right];G\left(q\right)=\left[\begin{array}{c}0\\ mg-mg{C}_{\theta }\\ mgl\left(t\right)S_{\theta }\end{array}\right];\hfill \\ \hfill & F=\left[\begin{array}{c}{F}_x\left(t\right)\\ F_l\left(t\right)\\ 0\end{array}\right];q=\left[\begin{array}{c}{x}_c\left(t\right)\\ l\left(t\right)\\ \theta \left(t\right)\end{array}\right];\hfill \end{array}$$

where $D_x$ and $D_y$ are the viscous damping coefficients associated with the horizontal and vertical axis, respectively; $M_x$, $M_l$ are the $x_c$ (traveling), l (hoisting) components of the crane mass and the equivalent masses of the rotating parts, such as the motors and their drive trains given by $M_x=m_c+m$ and $M_l=m$; m and $m_c$ are the load mass and cart mass, respectively; g is the gravitational acceleration, and for simplicity, the shorthand notation $S_{\theta }=sin\theta \left(t\right)$, $C_{\theta }=cos\theta \left(t\right)$ is used.

Equation (1) can be solved for $\ddot{q}$ as

$$\ddot{q}=-{M\left(q\right)}^{-1}D\dot{q}-{M\left(q\right)}^{-1}C\left(q,\dot{q}\right)\dot{q}-{M\left(q\right)}^{-1}G\left(q\right)+{M\left(q\right)}^{-1}F.$$

(2)

Then, by performing the algebraic products in (2), the mathematical model is computed as

$$\begin{array}{cc}\hfill \ddot{q}=& -\frac{1}{\xi \left(t\right)}\left[\begin{array}{ccc}\left(M_l+m\right)D_x& -m{S}_{\theta }{D}_l& 0\\ -m{S}_{\theta }{D}_x& \left(M_x+m{S}_{\theta }^2\right)D_l& 0\\ -\frac{C\left(M_l + m\right)D_x}{l\left(t\right)}& \frac{m{C}_{\theta }{S}_{\theta }{D}_l}{l\left(t\right)}& 0\end{array}\right]\dot{q}-\frac{1}{\xi \left(t\right)}\left[\begin{array}{ccc}0& 0& -m{M}_{l}l\left(t\right)S_{\theta }\dot{\theta }\left(t\right)\\ 0& 0& -m{M}_{x}l\left(t\right)\dot{\theta }\left(t\right)\\ 0& \frac{2\dot{\theta }\left(t\right)\xi \left(t\right)}{l\left(t\right)}& m{M}_l{C}_{\theta }{S}_{\theta }\dot{\theta }\left(t\right)\end{array}\right]\dot{q}\hfill \\ \hfill & -\frac{1}{\xi \left(t\right)}\left[\begin{array}{c}-mg\left(M_l{C}_{\theta }+m\right)S_{\theta }\\ \left(M_x\left(1-C_{\theta }\right)+m{S}_{\theta }^2\right)mg\\ \frac{g{S}_{\theta }\left(m^2{C}_{\theta } + M_x{M}_l + M_{x}m + M_{l}m\right)}{l\left(t\right)}\end{array}\right]\hfill \\ \hfill & +\frac{1}{\xi \left(t\right)}\left[\begin{array}{ccc}\left(M_l+m\right)& -m{S}_{\theta }& -\frac{C_{\theta }\left(M_l + m\right)}{l\left(t\right)}\\ -m{S}_{\theta }& \left(M_x+m{S}_{\theta }^2\right)& \frac{m{C}_{\theta }{S}_{\theta }}{l\left(t\right)}\\ -\frac{C_{\theta }\left(M_l + m\right)}{l\left(t\right)}& \frac{m{C}_{\theta }{S}_{\theta }}{l\left(t\right)}& \frac{M_x{M}_l + m\left(M_l + M_x\right) + m^2{C}_{\theta }^2}{ml{\left(t\right)}^2}\end{array}\right]\left[\begin{array}{c}{F}_x\left(t\right)\\ F_l\left(t\right)\\ 0\end{array}\right],\hfill \end{array}$$

(3)

with:

$$\begin{array}{c}\hfill \xi \left(t\right)=M_x{M}_l+M_{x}m+M_{l}m{S}_{\theta }^2.\end{array}$$

Setting $x\left(t\right)={[x_1\left(t\right),x_2\left(t\right),x_3\left(t\right),x_4\left(t\right),x_5\left(t\right),x_6\left(t\right)]}^T={[x_c\left(t\right),{\dot{x}}_c\left(t\right),l\left(t\right),\dot{l}\left(t\right),\theta \left(t\right),\dot{\theta }\left(t\right)]}^T$, $u\left(t\right)={[F_x\left(t\right),F_l\left(t\right)]}^T$ and considering that the designed control law will keep the load oscillations small makes it possible to assume that $\theta \approx 0$, i.e., $S_{\theta }\approx \theta$, $C_{\theta }\approx 1$, ${\theta }^2\left(t\right)\approx 0$, ${\dot{\theta }}^2\left(t\right)\approx 0$; moreover, $\theta$ must be considered small to reduce the number of nonlinearities involved in the system, and then, model (3) can be represented in a state-space nonlinear form as

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{x}}_1\left(t\right)\\ {\dot{x}}_2\left(t\right)\\ {\dot{x}}_3\left(t\right)\\ {\dot{x}}_4\left(t\right)\\ {\dot{x}}_5\left(t\right)\\ {\dot{x}}_6\left(t\right)\end{array}\right]=& \left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& -{m_{2}D}_x& 0& m_1{D}_l{x}_5\left(t\right)& m_{4}g& 0\\ 0& 0& 0& 1& 0& 0\\ 0& m_1{D}_x{x}_5\left(t\right)& 0& -{m_{3}D}_l& 0& 0\\ 0& 0& 0& 0& 0& 1\\ 0& \frac{m_2{D}_x}{x_3\left(t\right)}& 0& -\left(\frac{m_1{D}_l{x}_5\left(t\right)}{x_3\left(t\right)}+\frac{2x_6\left(t\right)}{x_3\left(t\right)}\right)& -\frac{m_{5}g}{x_3\left(t\right)}& 0\end{array}\right]\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ x_3\left(t\right)\\ x_4\left(t\right)\\ x_5\left(t\right)\\ x_6\left(t\right)\end{array}\right]\hfill \\ & +\left[\begin{array}{cc}0& 0\\ m_2& -m_1{x}_5\left(t\right)\\ 0& 0\\ -m_1{x}_5\left(t\right)& m_3\\ 0& 0\\ -\frac{m_2}{x_3\left(t\right)}& \frac{m_1{x}_5\left(t\right)}{x_3\left(t\right)}\end{array}\right]u\left(t\right);\hfill \end{array}$$

(4)

with

$$\begin{array}{cc}\hfill m_1=& \frac{m}{M_l{M}_x+M_{x}m},m_2=\frac{\left(M_l+m\right)}{M_l{M}_x+M_{x}m},m_3=\frac{M_x}{M_l{M}_x+M_{x}m},\hfill \\ \hfill m_4=& \frac{m\left(M_l+m\right)}{M_l{M}_x+M_{x}m},m_5=\frac{m^2+M_x{M}_l+M_{x}m+M_{l}m}{M_l{M}_x+M_{x}m}.\hfill \end{array}$$

Convex Linear Parameter Varying Model

To obtain a LPV representation, System (4) can be rewritten as

$$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& -{m_{2}D}_x& 0& {m_{1}D}_l{\rho }_1& m_{4}g& 0\\ 0& 0& 0& 0& 1& 0\\ 0& m_1{D}_x{\rho }_1& 0& -{m_{3}D}_l& 0& 0\\ 0& 0& 0& 0& 0& 1\\ 0& m_2{D}_x{\rho }_2& 0& -\left(m_1{D}_l{\rho }_3+2{\rho }_4\right)& -m_{5}g{\rho }_2& 0\end{array}\right]x\left(t\right)+\left[\begin{array}{cc}0& 0\\ m_2& -m_1{\rho }_1\\ 0& 0\\ -m_1{\rho }_1& m_3\\ 0& 0\\ -m_2{\rho }_2& m_1{\rho }_3\end{array}\right]u\left(t\right),\hfill \end{array}$$

(5)

where

$$\begin{array}{c}\hfill \rho =\left[{\rho }_1,{\rho }_2,{\rho }_3,{\rho }_4\right]=\left[x_5,\frac{1}{x_3},\frac{x_5}{x_3},\frac{x_6}{x_3}\right]\end{array}$$

(6)

are the nonlinear terms that will be considered as the scheduling variables. These nonlinear terms are chosen, such as ${\underline{\rho }}_k<{\rho }_k<{\overline{\rho }}_k,k=1,2,3,4$, where ${\rho }_k$ are the non-constants elements in system (5). The bounds of these scheduling variables are chosen according to physical limits or experimental constraints. In this case, these values are selected as $x_3\in [0.1,0.72]$ m, $x_5\in [-0.35,0.35]$ rad, $x_6\in [-3.467,3.467]$ rad/s. Note that these limits correspond to the rope’s length, the payload’s maximum/minimum expected oscillation, and the angular velocity of an experimental 3DOF crane system. Therefore, the scheduling variables are limited as

$$\begin{array}{cc}\hfill -0.35\le & {\rho }_1\le 0.35,\hfill \\ \hfill 1\le & {\rho }_2\le 10,\hfill \\ \hfill -3.5\le & {\rho }_3\le 3.5,\hfill \\ \hfill -34.67\le & {\rho }_4\le 34.67.\hfill \end{array}$$

For each ${\rho }_k$, two local scheduling functions are constructed as

$$w_0^k\left({\rho }_k\right)=\frac{{\overline{\rho }}_k-{\rho }_k}{{\overline{\rho }}_k-{\underline{\rho }}_k},w_1^k=1-w_0^k,k=1,2,3,4.$$

(7)

Therefore, for $k=4$, $i=2^k=16$ scheduling functions are computed as the product of the weighting functions that correspond to each local model:

$$h_i\left(\rho \left(t\right)\right)=\prod _{k=1}^4{\omega }_{ik}\left({\rho }_i\right).$$

(8)

The scheduling functions are convex which means that $h_i\left(\rho \left(t\right)\right)\ge 0$, ${\sum }_{i=1}^{16}{h}_i\left(\rho \left(t\right)\right)=1$. Finally, by considering the scheduling variables on (5), a linear parameter varying model in polytopic form is obtained as

$$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =\sum _{i=1}^{16}{h}_i\left(\rho \left(t\right)\right)\left[A_{i}x\left(t\right)+B_{i}u\left(t\right)\right];\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill y\left(t\right)& =Cx\left(t\right).\hfill \end{array}$$

(10)

where $A_i,B_i$, and C are constant matrices of appropriate dimensions obtained by evaluating (5) on the limits of ${\rho }_k$. The resulting linear models are not shown here due to space limitations, but they can be easily computed by considering the values of ${\rho }_k$ given above.

3. Stabilizing Controller

The overall proposed scheme for control and fault diagnosis is shown in Figure 2. In order to achieve system stabilization, an LPV state feedback controller is used, whose control law is defined by

$$u\left(t\right)=-\sum _{i=1}^{16}{h}_i\left(\rho \left(t\right)\right)K_{i}x\left(t\right),$$

(11)

by replacing the control law (11) in System (9), the closed-loop is represented as

$$\dot{x}\left(t\right)=\sum _{i=1}^{16}{h}_i\left(\rho \left(t\right)\right)\sum _{j=1}^{16}{h}_j\left(\rho \left(t\right)\right)\left(A_i-B_i{K}_j\right)x\left(t\right).$$

(12)

To guarantee asymptotic convergence of the closed-loop system, a Lyapunov candidate function is considered such as $V\left(x\left(t\right)\right)=x{\left(t\right)}^{T}Px\left(t\right)>0$, with $P\in {\mathbb{R}}^{n\times n},P>0,P=P^T$, and $V\left(0\right)=0$. The following theorem, obtained by solving the Lyapunov function, establishes sufficient conditions to guarantee the stabilization of the system.

Theorem 1.

The control law (11), with gains $K_i$, stabilizes the qLPV system (12), if and only if there are matrices $Q=Q^T>0$ and $W_j\forall j\in 1,2,\dots ,16$, such that the following LMI holds:

$$A_i{Q}^T-B_i{W}_j+Q{A}_i^T-W_j^T{B}_i^T<0.$$

(13)

Proof. 

The derivative of $V\left(x\right)$ along the trajectories of x is obtained as

$$\dot{V}\left(x\left(t\right)\right)=x{\left(t\right)}^{T}P\dot{x}\left(t\right)+\dot{x}{\left(t\right)}^{T}Px\left(t\right).$$

(14)

Then, by substituting (12) into (14) yields

$$\begin{array}{cc}\hfill \dot{V}\left(x\left(t\right)\right)& =\sum _{i=1}^{16}{h}_i\left(\rho \left(t\right)\right)\sum _{j=1}^{16}{h}_j\left(\rho \left(t\right)\right)[x{\left(t\right)}^{T}P{A}_{i}x\left(t\right)-x{\left(t\right)}^{T}P{B}_i{K}_{j}x\left(t\right)\hfill \end{array}$$

$$\begin{array}{cc}& +x{\left(t\right)}^T{A}_i^{T}Px\left(t\right)-x{\left(t\right)}^T{K}_j^T{B}_i^{T}Px\left(t\right)];\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\sum _{i=1}^{16}{h}_i\left(\rho \left(t\right)\right)\sum _{j=1}^{16}{h}_j\left(\rho \left(t\right)\right)x{\left(t\right)}^T\left(P{A}_i-P{B}_i{K}_j+A_i^{T}P-K_j^T{B}_j^{T}P\right)x\left(t\right).\hfill \end{array}$$

(16)

To fulfill that $\dot{V}\left(x\right)<0$, it is enough to prove that $P{A}_i-P{B}_i{K}_j+A_i^{T}P-K_j^T{B}_j^{T}P<0$. Nonetheless, to find a feasible solution in a linear matrix inequality (LMI) form, a substitution $Q\in {\mathbb{R}}^{n\times n},Q=Q^T=P^{-1}$ is performed, so the inequality is pre- and post-multiplied by Q and $Q^T$, respectively. After canceling the identities, it yields

$$\begin{array}{c}\hfill A_i{Q}^T-B_i{K}_j{Q}^T+Q{A}_i^T-Q{K}_j^T{B}_i^T<0.\end{array}$$

(17)

Since the term $B_i{K}_j{Q}^T$ is still a quadratic term, a substitution $W\in {\mathbb{R}}^{m\times n},W_j=K_{j}Q$ is considered, such as

$$A_i{Q}^T-B_i{W}_j+Q{A}_i^T-W_j^T{B}_i^T<0.$$

(18)

Therefore, by solving (13) the controller gains to stabilize the crane are computed. □

4. Proportional Integral Fault Estimation LPV Observer

Proportional-integral (PI) observers have become popular in recent years due to their robustness against disturbances due to the addition of a term that is proportional to the error in the estimation [31]. PI observers allow simultaneous estimation of system states and unknown inputs. In this case, the unknown inputs are actuator faults affecting the LPV model under the assumption that it is of slow variation, e.g., ${\dot{f}}_a\left(t\right)\approx 0$; in practice, it is possible to relax this condition [32,33]. In this case, faults present in the actuators are included in the analysis, with which the LPV representation of the plant is

$$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right)\left(A_{i}x\left(t\right)+B_{i}u\left(t\right)+M{f}_a\left(t\right)\right),\hfill \\ \hfill y\left(t\right)& =\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right)\left(C_{i}x\left(t\right)\right),\hfill \end{array}$$

(19)

where $f_a\left(t\right)\in {\mathbb{R}}^{n_f}$$y\left(t\right)\in {\mathbb{R}}^p$ represent actuator failures and measured outputs, respectively; moreover $A_i$, $B_i$, $C_i$, M are the constant matrices of appropriate dimensions. The structure of the PI observer is expressed as

$$\begin{array}{cc}\hfill \dot{\widehat{x}}\left(t\right)& =\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right)\left(A_i\widehat{x}\left(t\right)+B_{i}u\left(t\right)+L_{\mathcal{P},i}(y\left(t\right)-\widehat{y}\left(t\right))\right),\hfill \\ \hfill {\dot{\widehat{f}}}_a\left(t\right)& =\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right)\left(L_{\mathcal{I},i}(y\left(t\right)-\widehat{y}\left(t\right))\right),\hfill \\ \hfill \widehat{y}\left(t\right)& =\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right)\left(C_i\widehat{x}\left(t\right)\right),\hfill \end{array}$$

(20)

where the second line of this equation represents the fault dynamic, with $L_{P,i}$ and $L_{I,i}$ representing both the proportional and integral gains, respectively. The strategy to follow consists of grouping the states and faults into a single vector $\overline{x}\left(t\right)={\left[\begin{array}{cc}x{\left(t\right)}^T& f_a{\left(t\right)}^T\end{array}\right]}^T$, such that the system is rewritten as

$$\begin{array}{cc}\hfill \dot{\overline{x}}\left(t\right)& =\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right)\left({\overline{A}}_i\overline{x}\left(t\right)+{\overline{B}}_{i}u\left(t\right)\right),\hfill \\ \hfill y\left(t\right)& =\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right){\overline{C}}_{i}x\left(t\right),\hfill \end{array}$$

(21)

where ${\overline{A}}_i=\left[\begin{array}{cc}{A}_i& M\\ 0& 0\end{array}\right]$, ${\overline{B}}_i=\left[\begin{array}{c}{B}_i\\ 0\end{array}\right]$, and ${\overline{C}}_i=\left[\begin{array}{cc}{C}_i& 0\end{array}\right]$. Similarly, we rewrite the PI observer (20) in its augmented form:

$$\begin{array}{cc}\hfill \dot{\widehat{\overline{x}}}\left(t\right)& =\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right)\left({\overline{A}}_i\overline{x}\left(t\right)+{\overline{B}}_{i}u\left(t\right)+{\overline{L}}_i(y\left(t\right)-\widehat{y}\left(t\right))\right),\hfill \\ \hfill y\left(t\right)& =\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right){\overline{C}}_{i}x\left(t\right),\hfill \end{array}$$

(22)

where ${\overline{L}}_i$ are matrices that are composed of the proportional and integral gains ${\overline{L}}_i=\left[\begin{array}{c}{L}_{\mathcal{P},i}\\ L_{\mathcal{I},i}\end{array}\right]$. The error is calculated with the extended system (21) and the PI observer (22): $\overline{e}\left(t\right)=\overline{x}\left(t\right)-\widehat{\overline{x}}\left(t\right)$, whose dynamic is given by

$$\begin{array}{c}\hfill \dot{\overline{e}}\left(t\right)=\dot{\overline{x}}\left(t\right)-\dot{\widehat{\overline{x}}}\left(t\right).\end{array}$$

(23)

So the dynamics of the error is described as follows:

$$\dot{\overline{e}}\left(t\right)=\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right)\sum _{j=1}^r{h}_j\left(\rho \left(t\right)\right)\left({\overline{A}}_i-{\overline{L}}_j{\overline{C}}_i\right)\overline{e}\left(t\right).$$

(24)

With (24), the problem focuses on finding the appropriate gains ${\overline{L}}_j$ so that the PI observer reaches the behavior of the augmented system. The aim is to obtain LMI conditions to ensure that the error asymptotically converges to zero as ${lim}_{t\to \infty }\overline{e}\left(t\right)\approx 0$. The following theorem gives sufficient conditions to reach this goal.

Theorem 2.

The estimation error (24) is asymptotically stable if there exists a matrix $\overline{P}={\overline{P}}^T>0$ and gains ${\overline{L}}_j$, such that the LMI

$$\frac{2}{r-1}{\mathrm{Y}}_{ii}+{\mathrm{Y}}_{ij}+{\mathrm{Y}}_{ji}<0,$$

(25)

holds for every $(i,j)\in \{1,2,\dots ,r\}$ with

$${\mathrm{Y}}_{ij}:={\overline{A}}_i^T\overline{P}+\overline{P}{\overline{A}}_i-{\overline{W}}_j{\overline{C}}_i-{\overline{C}}_i^T{\overline{W}}_j^T+2\alpha \overline{P}.$$

The observer gains are calculated as ${\overline{L}}_j={\overline{P}}^{-1}{\overline{W}}_j$, $i\in \{1,2,\dots ,r\}$.

Proof. 

Consider a Lyapunov candidate function of the form

$$V\left(\overline{e}\left(t\right)\right)={\overline{e}}^T\left(t\right)\overline{P}\overline{e}\left(t\right),\mathrm{with} \overline{P}={\overline{P}}^T>0.$$

(26)

Then, the derivative of the Lyapunov candidate function is

$$\begin{array}{cc}\hfill \dot{V}\left(\overline{e}\left(t\right)\right)& =\left({\dot{\overline{e}}}^T\left(t\right)\overline{P}\overline{e}\left(t\right)+{\overline{e}}^T\left(t\right)\overline{P}\dot{\overline{e}}\left(t\right)\right)<0,\hfill \end{array}$$

(27)

and substituting the dynamics of the error (24),

$$\begin{array}{cc}\hfill \dot{V}\left(\overline{e}\left(t\right)\right)& ={\left(\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right)\sum _{j=1}^r{h}_j\left(\rho \left(t\right)\right)\left({\overline{A}}_i-{\overline{L}}_j{\overline{C}}_i\right)\overline{e}\left(t\right)\right)}^T\overline{P}\overline{e}\left(t\right)\hfill \\ \hfill & +{\overline{e}}^T\left(t\right)\overline{P}\left(\sum _{i=1}^r{h}_i\left(\rho \left(t\right)\right)\sum _{j=1}^r{h}_j\left(\rho \left(t\right)\right)\left({\overline{A}}_i-{\overline{L}}_j{\overline{C}}_i\right)\overline{e}\left(t\right)\right)<0,\hfill \end{array}$$

(28)

and developing (28), the following is obtained:

$$\dot{V}\left(\overline{e}\left(t\right)\right)\sum _{i=1}^r\sum _{j=1}^r{h}_i\left(\rho \left(t\right)\right)h_j\left(\rho \left(t\right)\right){\overline{e}}^T\left(t\right)\left({\overline{A}}_i^T\overline{P}-{\overline{C}}_i^T{\overline{L}}_j^T\overline{P}+\overline{P}{\overline{A}}_i-\overline{P}{\overline{L}}_j{\overline{C}}_i\right)\overline{e}\left(t\right)<0.$$

(29)

To eliminate the quadratic term, a change of variable $W_j=P{L}_j$ is made, obtaining the following LMI condition [34]:

$$A_i^{T}P+P{A}_i-W_j{C}_i-C_i^T{W}_j^T<0.$$

(30)

Within the context of convex models and qLPV systems, some relationships could assist LMIs to ensure a broader set of solutions. The relaxation lemma proposed by Tuan et al. [35] will be used to achieve this goal. Furthermore, to find the maximum possible for the associated Lyapunov function, a decay rate [36] is established, where it is established that if there exists a scalar $\alpha$ such that $\dot{V}\left(\overline{e}\left(t\right)\right)\le -2\alpha V\left(\overline{e}\left(t\right)\right)$, then the error states converge to the desired trajectories with a decay rate $\alpha$. Finally, the gains ${\overline{L}}_j$ conform as ${\overline{L}}_j=\left[\begin{array}{c}{L}_{\mathcal{P},j}\\ L_{\mathcal{I},j}\end{array}\right]$ providing sufficient LMI conditions to sustain (25), thus ending the proof. □

5. Numerical Results

This section is dedicated to illustrating the fault diagnosis method performance. For such a purpose, consider the crane parameters in

Table 1. Physical parameters of the crane.

Symbol	Value
g	9.81 m/s${}^2$
m	1 kg
$M_x$	3.49 kg
$M_l$	1 kg
$D_x$	100 Ns/m
$D_l$	82 Ns/m

. Since the convex qLPV model is just a convex rewriting of the nonlinear model, it is not necessary to perform a validation test as described by Bernal et al. [37].

The LMIs given in (13) are solved using SEDUMI [38] and YALMIP [39] Matlab toolboxes to compute the controller gains. A similar procedure is performed to compute the fault diagnosis observer gains within the LMIs given in (25), with $\alpha =1$, resulting in the following $\overline{P}$ matrix:

$$\overline{P}=\left[\begin{array}{ccccccc}0.3605& 0.0107& 0.278& 0.0788& 0.0514& -0.1484& -0.0252\\ 0.0107& 0.06374& 0.0888& 0.0264& -0.0224& 0& -0.0101\\ 0.278& 0.0888& 0.23488& 0.05229& 0.3254& 0& -0.2049\\ 0.0788& 0.0264& 0.05229& 0.0267& 0& 0& 0\\ 0.0514& -0.0224& 0.3254& 0& 0.40369& -0.2149& 0\\ -0.1484& 0& 0& 0& -0.2149& 1.6734& 0\\ -0.0252& -0.0101& -0.2049& 0& 0& 0& 2.2463\end{array}\right];$$

in addition, the necessary gains $L_i$ that guarantee the convergence of the observer have the following values:

$$\begin{array}{cc}\hfill {\overline{L}}_1& =\left[\begin{array}{ccc}29& 10.075& 63.757\\ 1934.3& 1213.7& 9296.9\\ -3.0068& 7.7705& -65.733\\ -80.936& 202.21& 916.79\\ 87.738& 57.218& 597.38\\ 15398& 10033& 103900\\ 12200& 8215.8& 50183\end{array}\right],{\overline{L}}_2=\left[\begin{array}{ccc}24.061& 8.1288& -7.7007\\ 1914.9& 1206.1& 9016.9\\ 0.13525& 9.0086& -20.278\\ 105.44& 275.67& 3613.6\\ 88.756& 57.62& 612.14\\ 15540& 10089& 105960\\ 12294& 8252.8& 51544\end{array}\right],\hfill \\ \hfill {\overline{L}}_3& =\left[\begin{array}{ccc}26.072& 8.921& 21.39\\ 1922.8& 1209.2& 9130.9\\ -1.1436& 8.5047& -38.784\\ 29.577& 245.77& 2515.9\\ 88.34& 57.456& 606.13\\ 15482& 10066& 105120\\ 12256& 8237.7& 50990\end{array}\right],{\overline{L}}_4=\left[\begin{array}{ccc}21.133& 6.9746& -50.075\\ 1903.4& 1201.5& 8850.8\\ 1.9985& 9.7429& 6.6774\\ 215.97& 319.22& 5212.8\\ 89.36& 57.858& 620.9\\ 15625& 10122& 107190\\ 12350& 8274.8& 52351\end{array}\right],\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\overline{L}}_5& =\left[\begin{array}{ccc}26.258& 8.9944& 24.073\\ 1629.5& 1093.6& 4869\\ -1.591& 8.3284& -45.052\\ -44.536& 216.56& 1439\\ 82.026& 54.968& 513.38\\ 14548& 9697.9& 91304\\ 10697& 7623.3& 28344\end{array}\right],{\overline{L}}_6=\left[\begin{array}{ccc}21.318& 7.0479& -47.407\\ 1610.2& 1086& 4588.9\\ 1.5508& 9.5666& 0.42363\\ 141.87& 290.02& 4136\\ 83.051& 55.37& 528.14\\ 14691& 9754.1& 93369\\ 10791& 7660.3& 29704\end{array}\right],\hfill \\ \hfill {\overline{L}}_7& =\left[\begin{array}{ccc}23.329& 7.8402& -18.303\\ 1618& 1089.1& 4703\\ 0.27219& 9.0626& -18.095\\ 65.988& 260.12& 3038.2\\ 82.631& 55.206& 522.13\\ 14632& 9731.2& 92528\\ 10752& 7645.3& 29151\end{array}\right],{\overline{L}}_8=\left[\begin{array}{ccc}18.39& 5.8938& -89.767\\ 1598.7& 1081.5& 4422.9\\ 3.4142& 10.301& 27.366\\ 252.38& 333.57& 5735.1\\ 83.651& 55.608& 536.89\\ 14775& 9787.5& 94593\\ 10846& 7682.3& 30511\end{array}\right],\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\overline{L}}_9& =\left[\begin{array}{ccc}12.82& -3.7051& 68.159\\ 2419.1& 1564.2& 9312.7\\ 7.8688& 16.967& -68.533\\ 653.91& 813.93& 750.09\\ 102.91& 68.79& 596.44\\ 17613& 11720& 103770\\ 15405& 10615& 50091\end{array}\right],{\overline{L}}_{10}=\left[\begin{array}{ccc}7.8815& -5.6515& -3.2993\\ 2399.7& 1556.5& 9032.6\\ 11.011& 18.205& -23.077\\ 840.29& 887.38& 3446.9\\ 103.93& 69.192& 611.2\\ 17756& 11777& 105830\\ 15499& 10652& 51452\end{array}\right],\hfill \\ \hfill {\overline{L}}_{11}& =\left[\begin{array}{ccc}9.8925& -4.8593& 25.801\\ 2407.6& 1559.6& 9146.6\\ 9.7322& 17.701& -41.592\\ 764.42& 857.48& 2349.2\\ 103.51& 69.028& 605.19\\ 17697& 11754& 104990\\ 15461& 10637& 50898\end{array}\right],{\overline{L}}_{12}=\left[\begin{array}{ccc}4.953& -6.8057& -45.67\\ 2388.2& 1552& 8866.6\\ 12.874& 18.939& 3.8756\\ 950.81& 930.94& 5046\\ 104.53& 69.43& 619.95\\ 17840& 11810& 107060\\ 15555& 10674& 52259\end{array}\right],\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\overline{L}}_{13}& =\left[\begin{array}{ccc}10.079& -4.7859& 28.487\\ 2114.3& 1444.1& 4884.8\\ 9.2849& 17.525& -47.863\\ 690.3& 828.28& 1272.3\\ 97.193& 66.54& 512.43\\ 16763& 11385& 91172\\ 13902& 10023& 28252\end{array}\right],{\overline{L}}_{14}=\left[\begin{array}{ccc}5.1393& -6.7323& -42.981\\ 2095& 1436.4& 4604.8\\ 12.427& 18.763& -2.3978\\ 876.69& 901.73& 3969.1\\ 98.215& 66.942& 527.2\\ 16905& 11442& 93236\\ 13996& 10060& 29613\end{array}\right],\hfill \\ \hfill {\overline{L}}_{15}& =\left[\begin{array}{ccc}7.1496& -5.9401& -13.895\\ 2102.9& 1439.5& 4718.7\\ 11.148& 18.259& -20.9\\ 800.83& 871.83& 2871.5\\ 97.8& 66.778& 521.19\\ 16847& 11419& 92396\\ 13957& 10045& 29059\end{array}\right],{\overline{L}}_{16}=\left[\begin{array}{ccc}2.2105& -7.8866& -85.358\\ 2083.5& 1431.9& 4438.7\\ 14.29& 19.497& 24.56\\ 987.22& 945.28& 5568.3\\ 98.82& 67.18& 535.95\\ 16990& 11475& 94461\\ 14051& 10082& 30419\end{array}\right].\hfill \end{array}$$

For simulation purposes, the proposed strategy will be applied to the nonlinear dynamics of the system described by Equations (9) and (10). In addition, the following initial conditions are considered: $x\left(0\right)={\left[\begin{array}{cccccc}0& 0.22& 0.22& 0& 0.1& 0.1\end{array}\right]}^T$ and $\widehat{x}\left(0\right)={\left[\begin{array}{cccccc}0.03& 0.2& 0.18& 0& 0.07& 0.11\end{array}\right]}^T$.

Because our primary goal is not to control but estimate the actuator faults, they are considered additive-type faults. The induced fault appears at $t=20\mathrm{s}$ and remains with a magnitude of $1\mathrm{N}$ until $t=40\mathrm{s}$ when it disappears. By monitoring the output signal of the actuator, it is possible to detect this fault from the estimation made by the observer described in (20).

The results are displayed in Figure 3, Figure 4 and Figure 5. Figure 3 and Figure 4 show the state estimation for all six states. The blue solid line represents the true state value, while the dashed red line represents the designed observer state estimation. It can be seen that the observer converges to the real values of the states in a negligible time for $x_1,x_2,x_5,$ and $x_6$. However, it takes a little longer for $x_3$ and $x_4$. Nevertheless, the observer estimates the states adequately despite the different initial conditions between the observer and the crane. Figure 5 shows the estimation of the fault affecting the system. It is relevant to highlight that the presence of an induced fault generates an additional demand on the PI observer, as illustrated in the figure. However, this demand is effectively addressed thanks to the strategy employed to calculate observer gains.

One of the notable strengths of the presented scheme is its ability to estimate fault behavior. The goal of understanding both the magnitude and dynamics of these faults is a significant advantage compared to other approaches. An example of the above is described by Guzmán-Rabasa et al. [40], where a fault detection method is shown based on obtaining residuals. Although this strategy is helpful for alerting about the presence of faults in the system and isolating the fault, it is essential to note that it does not provide the dynamics and magnitudes of these faults.

6. Conclusions

This document proposes a fault detection system based on a PI observer capable of estimating additive faults. The proposed methodology demonstrates its effectiveness in calculating the magnitude and dynamics of the fault, even in the presence of noise. To do this, the observer obtains its gains from a set of linear matrix inequalities, which, together with a decay rate, guarantee the asymptotic convergence. Then, the contribution of this work lies in the design of a PI observer based on a qLPV model for fault diagnosis with convergence guarantees. Measurement noise and parametric uncertainty will be considered for future work to obtain even more precise results when the system is affected by these factors. One of the shortcomings of this work is that the membership functions depend on measurable states, which is not a limitation for the 3DOF crane of this work. However, to improve the method’s applicability, future work will include membership functions based on estimated states. Another limitation to acknowledge is the current framework’s absence of external disturbance rejection mechanisms. Consequently, incorporating disturbance rejection techniques is a prospective avenue for future research as is incorporating a fault-tolerant control algorithm to take advantage of the information provided by the proposed fault diagnosis method.