# Actuator Fault Detection for Discrete-Time Descriptor Systems via a Convex Unknown Input Observer with Unknown Scheduling Variables

Authors: Javier Martinez; Braulio Aguiar; Víctor Estrada-Manzo; Miguel Bernal
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

This paper presents actuator fault detection of discrete-time nonlinear descriptor systems by means of nonlinear unknown input observers. The approach is based on the exact factorization of the estimation error in order to overcome the well-known problem of unmeasurable scheduling variables within the observation of convex models, thus avoiding the use of Lipschitz constants, differential mean value theorem, or robust techniques. As a result, the designing conditions are cast in terms of linear matrix inequalities and efficiently solved via commercially available software. Numerical as well as academic setups are provided to illustrate the advantages and performance of the proposal.


---

## Body

## 1. Introduction

The use of state observers [1] is important for many tasks within control theory. In particular, unknown input (UI) observers whose task is to estimate both the state and unknown inputs (they may be disturbances) have been developed in [2], and they are of particular interest also for fault detection and isolation [3, 4] or fault tolerant control schemes [5, 6]; they can be also applied in other fields [7]. For linear systems, there are several UI observers, e.g., [8–10]; for nonlinear systems, linear methods are mostly chosen because of its simplicity [11]. Nonlinear techniques such as sliding mode [12], adaptive schemes [13], Lipschitz approaches [14], high-gain observers [15], or combinations of them require certain structure of the model or performing nonlinear transformations.On the other hand, convex models (a convex model is a collection of linear models interconnected by scalar functions (also known as scheduling functions), which are nonlinear and hold the convex sum property in a region [16, 17]. If the convex model is the result of the sector nonlinearity [18], it is an exact representation of the nonlinear system) [19, 20] have been directly combined with the direct Lyapunov method; thus obtaining conditions in terms of linear matrix inequalities (LMIs) [21]. Conditions in the form of LMIs are preferred since they are numerically solvable via convex optimization techniques [22]. Within this context, there exist several state observers for both continuous [23, 24] and discrete-time systems [25–27]. In the case of UI observers, there are some works concerning proportional-integral setups [28, 29] or non-Luenberger forms [30–34].Nevertheless, within the convex framework, the observer design presents an open problem: if the scheduling variables do not exclusively depend on measurable/available signals, the designing conditions get involve and difficult to cast as LMIs. Recent works have intended to tackle this issue by employing Lipschitz constraints [24, 35], the differential mean value theorem (DMVT) [36, 37], and robust H∞ approaches to mitigate the influence of the unknown scheduling parameters [38]; more recently, in [39, 40], a transformation that enlarges the size of the state is proposed. In this work, we follow the ideas of [41] for solving this problem: algebraic rearrangements in order to properly factorize the error signal.A generalized representation for standard state-space models is given by descriptor systems [42]; within this setup, a special case considers that the descriptor matrix is invertible [43]. For this class of descriptor systems, in [27], a convex observer design has been developed, and it considers only available scheduling variables; in [44, 45], an UI observer has been proposed, and it can handle unmeasurable scheduling variables by robust argumentations.The estimation of the states and parameters makes possible the fault diagnosis, which is divided into fault detection, fault isolation, and fault identification [46, 47]. The fault detection (FD) is used to identify when there is a malfunction in the system and determine the moment when the fault occurs [48]. In FD, there exist some results using the UI observer for actuator fault, for example, [13] where an adaptive approach is developed for an aircraft actuator fault. In [49], UI observers are employed for fault detection and isolation.Contributions: a novel convex UI observer scheme is used for discrete-time descriptor models with unmeasurable scheduling variables; the scheme makes use of algebraic manipulations instead of Lipschitz constraints, the differential mean value theorem, or robust techniques in order to obtain an adequate estimation error dynamics, thus relaxing the results. Actuator fault detection and estimation is performed by means of the proposed UI observer.The rest of the paper is organized as follows: Section2 provides the background for further developments and notation. Section 3 states the LMI conditions for the nonlinear UI observer design via convex models. In Section 4, the UI observer is applied to actuator fault detection and estimation of nonlinear systems. Section 5 illustrates the proposal via the numerical example and train system. Section 6 concludes the paper by giving some final remarks.


## 2. Problem Statement

Consider the following discrete-time nonlinear descriptor model:(1)Eykxk+1=Axkxk+Bykuk+Dykdk,yk=Cxkxk,where xk∈ℝn is the state vector, uk∈ℝm is the input vector, and dk∈ℝq is the unknown input vector, yk∈ℝo is the output vector; Axk, Byk, Cxk, Dyk, and Eyk are matrix functions whose entries are smooth and bounded in a region Ωx⊂ℝn including the origin. This work only considers the particular case when E⋅ is full rank (the case when descriptor matrix is not invertible has been recently treated in [50, 51], this case is also referred as differential-algebraic-equation (DAE) systems or singular ones [52], and it is out of the scope of this work. Fault diagnosis schemes for this type of systems have been addressed in [53]) for xk∈Ωx; that is, from (1), it is always possible to obtain a standard state-space representation:(2)xk+1=E−1ykAxkxk+Bykuk+Dykdk=fxk,uk,dk.In what follows, arguments will be omitted when their meaning can be inferred from the context.The approach is based on designing an observer for the estimation of both the statexk and the unknown input dk; to this end, an augmented vector is employed [29, 44], that is, χk=xkTdkTT, dk=dkTdk+1T…dk+pTT∈ℝp+1q, where p is such that z−1p+1dk≈0 as proposed in [45]; for instance, for p=1, we have(3)dk+1=Sdk,withdk+1=dk+1dk+2,S=01−12.Hence, an augmented system yields(4)Eykχk+1=Axk,ykχk+ℬykuk,yk=Cxkχk,where χ∈ℝn+p+1q, S∈ℝpq+q×pq+q is a known matrix, and(5)In the literature, most of the observer design approaches deal with systems in standard form (2); for instance, without the unknown input, xk+1=fxk,uk and yk=Cxk; they consider special cases xk+1=Axk+Bfyk,uk+ϕxk; thus, the task is to stabilize an error system with the form ek+1=A−LCek+ϕxk−ϕx^k, where ek=xk−x^k is the estimation error, x^k is the estimated state vector, L is the observer gain to be designed, and the function ϕ⋅ is assumed to hold Lipschitz bounds, i.e., ϕx−ϕx^≤ℒx−x^, ℒ>0 being a Lipschitz constant [54, 55]. In the context of convex models, most of the works only consider that the scheduling variables are available and then ϕxk−ϕx^k=0 [25], which, in practice, is not realistic. For the general case, when ϕxk−ϕx^k≠0, this expression is treated by means for Lipschitz bounds [24], as a perturbation via H∞ approach [38], as an uncertainty via robust approaches [56]; other works employ the differential mean value theorem [37], Jacobian [57], or transformations that enlarge the size of the state [40]. Nevertheless, these approaches are conservative approximations, or increase the computational complexity of the problem, or are only valid for particular cases [58]; additionally, none of them consider observers for systems of form (1). Next section presents a methodology that avoids the use of previous ones.


### 2.1. An Amenable Descriptor Error Form

In order to obtain an amenable error system for Lyapunov analysis, in [41], a methodology that avoids the use of Lipschitz-like bounds has been presented; it is based on algebraic operations in order to factorize the estimation error at the left side of error dynamic difference equation. Thus, the variation of the Lyapunov function along the trajectories of the error system can be written as ΔVe=ekTQ⋅ek and ΔVe<0 if Q⋅<0, and the latter is guaranteed via convex models and linear matrix inequalities. Motivated by these ideas, the following observer structure is adopted:(6)Eykχ^k+1=Ax^k,ykχ^k+ℬykuk+Lx^k,ykyk−y^k,y^k=Cx^kχ^k.with(7)and the nonlinear observer gain Lx^k,yk∈ℝn+pq+q×o depending only on all the available signals; it should be designed such that the estimation error(8)ek=χk−χ^k=xk−x^kdk−d^ksatisfies limk⟶∞ek=0. Thus, in [41], it has been proven that, under mild assumptions, it is always possible to write the error dynamics as(9)Eykek+1=A¯xk,x^k−Lx^k,ykC¯xk,x^kek,where A¯xk,x^kek=Axk,ykχk−Ax^k,ykχ^k and C¯xk,x^kek=Cxχk−Cx^χ^k have bounded entries in Ωx×Ωx^. For instance, consider a polynomial expression px−px^ with px=x1x2 and px^=x^1x^2; following [41], we have px−px^=0.5x2+x^2e1+0.5x1+x^1e2, e1=x1−x^1, e2=x2−x^2. Now, consider a nonpolynomial expression p˜x−p˜x^, with p˜x=sinx1 and p˜x^=sinx^1; a third-order Taylor approximation will give p˜x−p˜x^≈x1−x13/6−x^1+x^13/6, and then we can apply similar procedure as before.Now, error system (9) has an amenable form for Lyapunov-based analysis. However, if the aim is to obtain LMI conditions, (9) should be expressed as an exact convex model, and this is the matter of the following section.


### 2.2. Convex Expressions

The sector nonlinearity [18] is employed to express bounded nonconstant terms z⋅∈z0,z1 as a convex sums of its bounds, that is, z⋅=w0zz0+w1zz1, where z0 and z1 are the minimum and maximum of z⋅ in a region; the functions w0=z1−z⋅/z1−z0 and w1=1−w0 hold the convex sum property in the same modeling region, i.e., w0z+w1z=1 and w0,w1∈0,1.Note that, in (9), matrices Eyk, A¯xk,x^k, and C¯xk,x^k contain nonconstant terms depending on xk, x^k, and yk; clearly, all the state variables are not fully available; thus, a useful convex rewriting must take this into account. The following steps extend the sector nonlinearity to our case:Step 1: identify all the nonconstant terms, also known as scheduling variables, inEyk, A¯xk,x^k, and C¯xk,x^k depending exclusively on available signals, and capture them in the vector zx^k,yk∈ℝs while all the rest of the terms should be grouped in ζxk,x^k,yk∈ℝσ. Each entry is assumed to be bounded in Ωx×Ωx^, i.e., zix^k,yk∈zi0,zi1, i∈1,2,…,s, and ζjxk,x^k,yk∈ζj0,ζj1, j∈1,2,…,σ.Step 2: construct, for eachzix^k,yk, i∈1,2,…,s, and ζjxk,x^k,yk, j∈1,2,…,σ a pair of scalar convex functions as follows:(10)w0ix^,y=zi1−zix^,yzi1−zi0,w1ix^,y=1−w0ix^,y,ω0jx,x^,y=ζj1−ζjx,x^,yζj1−ζj0,ω1jx,x^,y=1−ω0jx,x^,y;by construction, these functions hold the convex sum property inΩx×Ωx^, i.e., w0ix^,y+w1ix^,y=1, w0ix^,y,w1ix^,y∈0,1, ω0jx,x^,y+ω1jx,x^,y=1, ω0jx,x^,y,ω1jx,x^,y∈0,1.Step 3: define the scheduling functions as(11)wiz=wi1jz1wi22z2,…,wisszs,ωjζ=ωj11ζ1ωj22ζ2,…,ωjσσζσ,withi∈1,2,…,r, j∈1,2,…,ρij∈0,1, r=2s, ρ=2σ; moreover, the sets of indexes i1i2,…,is and j1j2,…,jσ are a s-digit and σ-digit binary representation of i−1 and j−1, respectively. The scheduling functions also hold the convex sum property in Ωx×Ωx^, that is, ∑i=1rwiz=1, wiz∈0,1, ∑j=1ρωjζ=1, and ωjζ∈0,1.Step 4: compute the vertex modelsEi=Eyk|wi=1, A¯ij=A¯xk,x⌢k|wiωj=1, C¯ij=C¯xk,x^k|wiωj=1, i∈1,2,…,r, j∈1,2,…,ρ.Thus, an exact convex representation of (9) is(12)∑i=1rwizkEiek+1=∑i=1r∑j=1ρwizkωjζkA¯ij−Lx^k,ykC¯ijek.The nonlinear observer gain isLx^k,yk and will be defined later on.


#### 2.2.1. Notation

For convex expressions, the following shorthand notation will be employed throughout the manuscript: single convex sumsΥwzk=∑i=1rwizkΥi and its inverse Υwzk−1=∑i=1rwizkΥi−1, with delayed scheduling functions Υwzk+1=∑m=1rwmzk+1Υm, or depending on nonavailable variables Υωζk=∑j=1ρωjζkΥj, and so on. Additionally, A>0<0 means that A∈ℝn×n is positive (negative) definite. An asterisk ∗ will be used in matrix expressions to denote the transpose of the symmetric element; for in-line expressions, it will denote the transpose of the terms on its left side:(13)ABTBC=A∗BCandA+B+AT+BT+C=A+B+∗+C.Hence, system (12) is shortly written as(14)Ewzkek+1=A¯wzkωζk−Lx^k,ykC¯wzkωζkek.The following lemmas are useful in order to derive LMI conditions for the design ofLx^k,yk. The first one concerns a sum-relaxation scheme based on [59]; the second one allows avoiding the computation of E−1yk while adding slack variables [27].Lemma 1.
(see [59]). Let Υiljm=ΥiljmT, i,l,m∈1,2,…,r3, j∈1,2,…,ρ, be matrices of adequate dimensions. Then, Υwzkwzkwzk+1ωζk<0 holds if(15)2r−1Υiijm+Υiljm+Υlijm<0,for all i,l,m∈1,2,…,r3, j∈1,2,…,ρ.Lemma 2.
(see [60]). Let ξ∈ℝn, Q=QT∈ℝn×n, and ℬ∈ℝm×n, rankℬ<n; then, the following statements are equivalent:(i)
ξTQξ<0, ∀ℬξ=0, ξ≠0(ii)
∃Z∈ℝn×m:Q+Zℬ+ℬTZT<0
Now, we are ready to establish LMI conditions for the stabilization of error system (9) at the origin e=0 via its exact convex representation (12).


## 2.1. An Amenable Descriptor Error Form

In order to obtain an amenable error system for Lyapunov analysis, in [41], a methodology that avoids the use of Lipschitz-like bounds has been presented; it is based on algebraic operations in order to factorize the estimation error at the left side of error dynamic difference equation. Thus, the variation of the Lyapunov function along the trajectories of the error system can be written as ΔVe=ekTQ⋅ek and ΔVe<0 if Q⋅<0, and the latter is guaranteed via convex models and linear matrix inequalities. Motivated by these ideas, the following observer structure is adopted:(6)Eykχ^k+1=Ax^k,ykχ^k+ℬykuk+Lx^k,ykyk−y^k,y^k=Cx^kχ^k.with(7)and the nonlinear observer gain Lx^k,yk∈ℝn+pq+q×o depending only on all the available signals; it should be designed such that the estimation error(8)ek=χk−χ^k=xk−x^kdk−d^ksatisfies limk⟶∞ek=0. Thus, in [41], it has been proven that, under mild assumptions, it is always possible to write the error dynamics as(9)Eykek+1=A¯xk,x^k−Lx^k,ykC¯xk,x^kek,where A¯xk,x^kek=Axk,ykχk−Ax^k,ykχ^k and C¯xk,x^kek=Cxχk−Cx^χ^k have bounded entries in Ωx×Ωx^. For instance, consider a polynomial expression px−px^ with px=x1x2 and px^=x^1x^2; following [41], we have px−px^=0.5x2+x^2e1+0.5x1+x^1e2, e1=x1−x^1, e2=x2−x^2. Now, consider a nonpolynomial expression p˜x−p˜x^, with p˜x=sinx1 and p˜x^=sinx^1; a third-order Taylor approximation will give p˜x−p˜x^≈x1−x13/6−x^1+x^13/6, and then we can apply similar procedure as before.Now, error system (9) has an amenable form for Lyapunov-based analysis. However, if the aim is to obtain LMI conditions, (9) should be expressed as an exact convex model, and this is the matter of the following section.


## 2.2. Convex Expressions

The sector nonlinearity [18] is employed to express bounded nonconstant terms z⋅∈z0,z1 as a convex sums of its bounds, that is, z⋅=w0zz0+w1zz1, where z0 and z1 are the minimum and maximum of z⋅ in a region; the functions w0=z1−z⋅/z1−z0 and w1=1−w0 hold the convex sum property in the same modeling region, i.e., w0z+w1z=1 and w0,w1∈0,1.Note that, in (9), matrices Eyk, A¯xk,x^k, and C¯xk,x^k contain nonconstant terms depending on xk, x^k, and yk; clearly, all the state variables are not fully available; thus, a useful convex rewriting must take this into account. The following steps extend the sector nonlinearity to our case:Step 1: identify all the nonconstant terms, also known as scheduling variables, inEyk, A¯xk,x^k, and C¯xk,x^k depending exclusively on available signals, and capture them in the vector zx^k,yk∈ℝs while all the rest of the terms should be grouped in ζxk,x^k,yk∈ℝσ. Each entry is assumed to be bounded in Ωx×Ωx^, i.e., zix^k,yk∈zi0,zi1, i∈1,2,…,s, and ζjxk,x^k,yk∈ζj0,ζj1, j∈1,2,…,σ.Step 2: construct, for eachzix^k,yk, i∈1,2,…,s, and ζjxk,x^k,yk, j∈1,2,…,σ a pair of scalar convex functions as follows:(10)w0ix^,y=zi1−zix^,yzi1−zi0,w1ix^,y=1−w0ix^,y,ω0jx,x^,y=ζj1−ζjx,x^,yζj1−ζj0,ω1jx,x^,y=1−ω0jx,x^,y;by construction, these functions hold the convex sum property inΩx×Ωx^, i.e., w0ix^,y+w1ix^,y=1, w0ix^,y,w1ix^,y∈0,1, ω0jx,x^,y+ω1jx,x^,y=1, ω0jx,x^,y,ω1jx,x^,y∈0,1.Step 3: define the scheduling functions as(11)wiz=wi1jz1wi22z2,…,wisszs,ωjζ=ωj11ζ1ωj22ζ2,…,ωjσσζσ,withi∈1,2,…,r, j∈1,2,…,ρij∈0,1, r=2s, ρ=2σ; moreover, the sets of indexes i1i2,…,is and j1j2,…,jσ are a s-digit and σ-digit binary representation of i−1 and j−1, respectively. The scheduling functions also hold the convex sum property in Ωx×Ωx^, that is, ∑i=1rwiz=1, wiz∈0,1, ∑j=1ρωjζ=1, and ωjζ∈0,1.Step 4: compute the vertex modelsEi=Eyk|wi=1, A¯ij=A¯xk,x⌢k|wiωj=1, C¯ij=C¯xk,x^k|wiωj=1, i∈1,2,…,r, j∈1,2,…,ρ.Thus, an exact convex representation of (9) is(12)∑i=1rwizkEiek+1=∑i=1r∑j=1ρwizkωjζkA¯ij−Lx^k,ykC¯ijek.The nonlinear observer gain isLx^k,yk and will be defined later on.


### 2.2.1. Notation

For convex expressions, the following shorthand notation will be employed throughout the manuscript: single convex sumsΥwzk=∑i=1rwizkΥi and its inverse Υwzk−1=∑i=1rwizkΥi−1, with delayed scheduling functions Υwzk+1=∑m=1rwmzk+1Υm, or depending on nonavailable variables Υωζk=∑j=1ρωjζkΥj, and so on. Additionally, A>0<0 means that A∈ℝn×n is positive (negative) definite. An asterisk ∗ will be used in matrix expressions to denote the transpose of the symmetric element; for in-line expressions, it will denote the transpose of the terms on its left side:(13)ABTBC=A∗BCandA+B+AT+BT+C=A+B+∗+C.Hence, system (12) is shortly written as(14)Ewzkek+1=A¯wzkωζk−Lx^k,ykC¯wzkωζkek.The following lemmas are useful in order to derive LMI conditions for the design ofLx^k,yk. The first one concerns a sum-relaxation scheme based on [59]; the second one allows avoiding the computation of E−1yk while adding slack variables [27].Lemma 1.
(see [59]). Let Υiljm=ΥiljmT, i,l,m∈1,2,…,r3, j∈1,2,…,ρ, be matrices of adequate dimensions. Then, Υwzkwzkwzk+1ωζk<0 holds if(15)2r−1Υiijm+Υiljm+Υlijm<0,for all i,l,m∈1,2,…,r3, j∈1,2,…,ρ.Lemma 2.
(see [60]). Let ξ∈ℝn, Q=QT∈ℝn×n, and ℬ∈ℝm×n, rankℬ<n; then, the following statements are equivalent:(i)
ξTQξ<0, ∀ℬξ=0, ξ≠0(ii)
∃Z∈ℝn×m:Q+Zℬ+ℬTZT<0
Now, we are ready to establish LMI conditions for the stabilization of error system (9) at the origin e=0 via its exact convex representation (12).


## 2.2.1. Notation

For convex expressions, the following shorthand notation will be employed throughout the manuscript: single convex sumsΥwzk=∑i=1rwizkΥi and its inverse Υwzk−1=∑i=1rwizkΥi−1, with delayed scheduling functions Υwzk+1=∑m=1rwmzk+1Υm, or depending on nonavailable variables Υωζk=∑j=1ρωjζkΥj, and so on. Additionally, A>0<0 means that A∈ℝn×n is positive (negative) definite. An asterisk ∗ will be used in matrix expressions to denote the transpose of the symmetric element; for in-line expressions, it will denote the transpose of the terms on its left side:(13)ABTBC=A∗BCandA+B+AT+BT+C=A+B+∗+C.Hence, system (12) is shortly written as(14)Ewzkek+1=A¯wzkωζk−Lx^k,ykC¯wzkωζkek.The following lemmas are useful in order to derive LMI conditions for the design ofLx^k,yk. The first one concerns a sum-relaxation scheme based on [59]; the second one allows avoiding the computation of E−1yk while adding slack variables [27].Lemma 1.
(see [59]). Let Υiljm=ΥiljmT, i,l,m∈1,2,…,r3, j∈1,2,…,ρ, be matrices of adequate dimensions. Then, Υwzkwzkwzk+1ωζk<0 holds if(15)2r−1Υiijm+Υiljm+Υlijm<0,for all i,l,m∈1,2,…,r3, j∈1,2,…,ρ.Lemma 2.
(see [60]). Let ξ∈ℝn, Q=QT∈ℝn×n, and ℬ∈ℝm×n, rankℬ<n; then, the following statements are equivalent:(i)
ξTQξ<0, ∀ℬξ=0, ξ≠0(ii)
∃Z∈ℝn×m:Q+Zℬ+ℬTZT<0
Now, we are ready to establish LMI conditions for the stabilization of error system (9) at the origin e=0 via its exact convex representation (12).


## 3. LMI-Based Stabilization of the Error System

This section provides LMI conditions to computeLx^k,yk. The developments are based on a convex Lyapunov function candidate (in the context of TS models, it has been introduced as nonquadratic Lyapunov function [61], in the context of fuzzy systems as fuzzy Lyapunov functions [62], in the context of LPV models, as parameter dependent Lyapunov functions [63]. Here, the name convex Lyapunov function is adopted due to its dependence on the convex functions wizk):(16)Ve=ekTPwzkek,Pwzk=∑i=1rwizkPi,with Pi∈ℝn+pq+q×n+pq+q:Pi>0, i∈1,2,…,r; its variation is(17)ΔVe=ek+1TPwzk+1ek+1−ekTPwzkek,which can be expressed as follows (without substituting the dynamics of estimation error (9)):(18)ΔVe=ekek+1T−Pwzk00Pwzk+1ekek+1.Note that function (16) is convex and it only depends on available signals. With this in mind, the following result states LMI conditions for the design of the nonlinear gain Lx^k,yk.Theorem 1.
The origine=0 of error system (9), with an exact convex representation (12), is asymptotically stable if there exist matrices Pl∈ℝn+pq+q×n+pq+q, Nl∈ℝn+pq+q×o, and Gl∈ℝn+pq+q×n+pq+q, l∈1,2,…,r such that Pl>0 and the LMIs in (15) are satisfied with(19)Υiljm≔−Pl∗GlA¯ij−NlC¯ij−GlEi−EiTGlT+Pm,for all i,l,m∈1,2,…,r3, j∈1,2,…,ρ. Then, the observer gain is computed as Lx^k,yk=Gwzk−1Nwzk. Moreover, any trajectory ek starting in the outermost Lyapunov level set e:Vek≤c⊂Ωx×Ωx^,c>0 goes to zero as time goes to infinity.Proof.
Recall the variation of the Lyapunov function expressed as in (18); thus, error system (9), with exact convex representation (12), is also expressed as(20)A¯wzkωζk−Lx^k,ykC¯wzkωζk−Ewzkekek+1=0.
By Lemma2, (18) and (20) can be written together:(21)Z1Z2A¯wzkωζk−Lx^k,ykC¯wzkωζk−Ewzk+∗+−Pwzk00Pwzk+1<0.
Hence, by choosingZ1=0 and Z2=Gwzk and the definition Lx^k,yk=Gwzk−1Nwzk, Gwzk∈ℝpq+q+n×pq+q+n, Nwzk∈ℝpq+q+n×o, and we have(22)Υwzkwzkwzk+1ωζk≔−Pwzk∗GwzkA¯wzkωζk−NwzkC¯wzkωζk−GwzkEwzk+∗+Pwzk+1<0.
Finally, by means of relaxation Lemma1, the desired result yields.Remark 1.
The speed convergence of observer (6) can be increased if the following condition is verified ΔVe≤α2−1Ve, 0<α≤1 which can be translated into LMIs, that is, solving the LMIs in (15) with(23)Υiljm≔−α2Pl∗GlA¯ij−N¯lij−GlEi−EiTGlT+Pm,for all i,l,m∈1,2,…,r3, j∈1,2,…,ρ.
The following result provides LMI conditions for standard systems and its corresponding observer, i.e., (6) with Eyk=In, where In is the identity matrix. Hence, the next result follows directly from Theorem 1.Corollary 1.
The origin of error system (9), with Eyk=In, is asymptotically stable if there exist Pl∈ℝn+pq+q×n+pq+q, Nl∈ℝn+pq+q×o, and Gl∈ℝn+pq+q×n+pq+q, l∈1,2,…,r such that Pl>0 and LMIs in (15) hold with(24)Υiljm=−Pl∗GlA¯ij−NlC¯ij−Gl−GlT+Pm,for all i,l,m∈1,2,…,r3, j∈1,2,…,ρ. The nonlinear observer gain Lx^k,yk=Gwzk−1Nwzk. Moreover, any trajectory ek starting in the outermost Lyapunov level set e:Ve≤c⊂Ωx×Ωx^,c>0 goes to zero as time goes to infinity.Proof.
It follows directly by consideringEyk=In from the developments in Theorem 1.Remark 2.
The methodology in previous approaches starts by computing a convex model of the given nonlinear one (standard or descriptor), and then, the estimation error is computed; if the problem involves unmeasurable scheduling variables, it is solved by considering the already defined scheduling functions. In contrast with those works, the methodology hereby presented begins by obtaining the nonlinear estimation error dynamics via factorizations; then, it employs the sector nonlinearity approach for writing an equivalent convex representation; thus, a fitter convex model is employed and the observer gain includes all the available signals.Remark 3.
In this work, a convex Lyapunov function such as the one introduced in [61] is employed; naturally delayed scheduling functions can be included similar to [25] or to the generalization [27]. For instance, following [25], the Lyapunov function candidate would be Ve=ekTPwzk−1ek, Pwzk−1=∑m=1rwmzk−1Pm together with the observer gain Lx^k,x^k−1,yk,yk−1=Gwzkwzk−1−1Nwzkwzk−1.Remark 4.
The numerical complexity of Theorem1 can be approximated by log10nd3nl where nd=rn+pq+q0.5n+pq+q+1+n+pq+q+o is the number of decision variables and nl=2r3ρn+pq+q is the number of LMI rows [64]. Moreover, if the number LMI conditions in Theorem 1 is large, one can reduce its complexity by “judiciously eliminating” some convex functions either from the Lyapunov function or the observer gain.
In the case without unknown inputsdk=0 in (1), the states of the resulting system can be estimated by the following nonlinear observer:(25)Eykx^k+1=Ax^k,ykx^k+Bykuk+Lx^k,ykyk−y^k,y^k=Cx^k,ykx^k,where x^k∈ℝn is the observer state and Lx^k,yk∈ℝn×o is the observer gain to be designed. Following the methodology previously presented, the error dynamics is(26)Eykek+1=A¯xk,x^k−Lx^k,ykC¯xk,x^kek,where A¯xk,x^kek=Axkxk−Ax^k,ykx^k and C¯xk,x^kek=Cxkxk−Cx^k,ykx^k. Thus, we have the next result.Corollary 2.
The origin of error system (26) is asymptotically stable if there exist Pl∈ℝn×n, Nl∈ℝn×o, and Gl∈ℝn×n, l∈1,2,…,r such that Pl>0 and LMIs in (15) hold with(27)Υiljm=−Pl∗GlA¯ij−NlC¯ij−GlEi−EiTGlT+Pm,for all i,l,m∈1,2,…,r3, j∈1,2,…,ρ. The nonlinear observer gain Lx^k,yk=Gwzk−1Nwzk. Moreover, any trajectory ek starting in the outermost Lyapunov level set e:Ve≤c⊂Ωx×Ωx^,c>0 goes to zero as time goes to infinity.Proof.
The proof follows a similar path than the one for Theorem1, with a Lyapunov function candidate Ve=ekTPwzkek, Pi>0, Pi∈ℝn×n, i∈1,2,…,r, and the nonlinear error dynamics (26) expressed in a convex form as(28)A¯wzkωζk−Lx^k,ykC¯wzkωζk−Ewzkekek+1=0,and by the use of Finsler’s lemma, with Lx^k,yk=Gwzk−1Nwzk, Z1=0, and Z2=Gwzk, we have(29)Υwzkwzkwzk+1ωζk≔−Pwzk∗GwzkA¯wzkωζk−NwzkC¯wzkωζk−GwzkEwzk+∗+Pwzk+1<0.
Finally, by means of relaxation Lemma1, we conclude the proof.
Recall that most of the previous works do not consider scheduling functions depending on unmeasurable variables [25, 27] nor descriptor systems of form (1) [25].


## 4. Actuator Fault Detection and Estimation

In this section, an application of the UI observer developed above is employed in order to estimate and detect actuator faults [46, 47]. Thus, let us now consider nonlinear discrete-time descriptor system (1) with an actuator fault:(30)Eykxk+1=Axkxk+Bykuk+fauk,yk+Dykdk,yk=Cxkxk,where fauk,yk∈ℝm represents the actuator fault vector; if it can be included inside of the unknown input, then (30) can be expressed as (1); in this case, conditions in Theorem 1 can be used. As customary in fault detection schemes, a residual signal based on the observer estimation of the unknown input is generated [65, 66]:(31)rk=d^k.The residual can be filtered to get a clear signal. The estimation of the fault can be done once the fault is considered as an unknown input [67].Remark 5.
It is possible to apply a mix ofH∞ and H_ performances, for instance, H∞ attenuation in order to make residual (31) robust to uncertainties, noise, and the fault and the index H_ index in order to make more sensitive the detection of the fault [68, 69].


## 5. Examples

In this section, two examples are presented. The first one is intended to compare the performance of our proposal in contrast with recent approaches while the second one illustrates the estimation of actuator faults in the train system. The LMI conditions have been implemented in YALMIP [70] with SeDuMi [71], while simulations have been run in Simulink for MATLAB2019b.Example 1.
Consider descriptor system (1) with matrices(32)Eyk=2−1+x12−11+x12−11,Axk=x22−0.5β0.7sinx1,B=0.10.1,D=0.20.2,C=10T,with β>0 is a known parameter. It is assumed that the unknown input has dynamics such as z−12dk≈0, i.e., dk+1=Sdk with S=01−12. In order to estimate both the state xk and the unknown input dk, an extended system of form (4) and its corresponding observer (6) are constructed. Thus, the error dynamics ek=χk−χ^k can be computed as(33)Eykek+1=x1x22−0.5βx2+0.2d10.7x1+x2sinx1+0.2d1d2−d1+2d2−x^1x^22−0.5βx^2+0.2d^10.7x^1+x^2sinx^1+0.2d^1d^2−d^1+2d^2−L⋅y−y^.
Applying the factorization methodology in [41] and considering that y=x1, we have x22−x^22=x2+x^2x2−x^2=x2+x^2e2 leading to final error dynamics(34)2−1+x12−1001+x12−110000100001︸Eykek+1=0x1x2+x^2−0.5β0.200.7sinx10.20000100−12︸A¯xk,x^k−L⋅1000︸C¯ek.
In order to synthesize the nonlinear observer gain via LMIs, it is necessary to express (34) in a convex form (12); to this end, let us consider the compact sets Ωx=x1≤0.5,x2≤0.8 and Ωx^=x^2≤0.8. Thus, the nonconstant terms and their bounds are z1=1+x12−1∈0.8,1, z2=sinx1∈−0.4794,0.4794, z3=x1∈−0.5,0.5, z4=x^2∈−0.8,0.8, and ζ1=x2∈−0.8,0.8 (nonavailable signals). The scalar convex functions are wiz=wi11z1wi22z2wi33z3wi44z4, w0izi=zi1−zi/zi1−zi0, w1izi=1−w0izi, i∈1,2,3,4, and ωjζ=ωj11ζ1, ω01ζ1=ζ11−ζ1/ζ11−ζ10, ω11=1−ω01. The vertex matrices are(35)Ei=2−z1i100z1i110000100001,A¯ij=0z3i3ζ1j1+z4i4−0.5β0.200.7z2i20.20000100−12.
For example, ifi=6 and j=1, we have w6z=w01z1w12z2w03z3w14z4, ω1ζ=ω01ζ1, and(36)E6=2−z1000z1010000100001=2−0.8000.810000100001,A¯6,1=0z30ζ10+z41−0.5β0.200.7z210.20000100−12=0−0.5β0.200.70.47940.20000100−12.
Several numerical tests have been performed in order to show the effectiveness of the proposal in contrast with recent works:(i)
LMIs in Theorem1 have been run as well as the ones in [45], Theorem 2, with a Lyapunov function candidate (16) and decay rate α=0.95 (see Remark 1), seeking feasibility for the largest β>0. It results that our proposal is feasible up to β=0.863, while the one in [45], Theorem 2, is only up to β=0.794. Comparisons have been done with the same Lyapunov function (16) and under the same relaxation scheme.(ii)
In terms of the number of decision variables, for this example, we have that numerical complexity for Theorem1 is 12.86 while for [45] is 10.24.(iii)
A UI observer uses the generalization in [27], and although it is for discrete-time descriptor systems, it cannot be applied as it only considers available scheduling variables.(iv)
If UI observer schemes for standard state-space models are to be applied [33, 39], then the inverse of the matrix Eyk has to be computed. For instance, using conditions in Corollary 1 yields numerical problems because 65568 LMIs have to be solved.
For illustration purposes, the LMI conditions in form (1) have been found feasible for β=0.85 and a decay rate α=0.95, and some of the computed matrices are(37)N1=1×10−40.35020.44550.0786−0.0245,G1=1×10−40.45280.4482−0.1117−0.0147−0.41970.68390.0704−0.07720.0656−0.30330.7189−0.5508−0.07960.1505−0.57510.4700,N5=1×10−40.49360.52950.0875−0.0429,G5=1×10−40.37100.5399−0.22360.0594−0.44580.63160.0105−0.04170.0960−0.21680.6811−0.5264−0.09610.0893−0.54220.4492,N16=1×10−4−0.13300.75690.1866−0.0911,G16=1×10−40.8122−0.32600.0464−0.0970−0.08140.9592−0.20180.0460−0.0773−0.10130.6744−0.5174−0.00470.0364−0.54190.4463.
A first simulation has been performed without unknown inputs, i.e.,dk=0 and the known input uk=0.4sin0.6k. Figure 1 shows evolution of the error signals converging asymptotically to the origin. A second simulation considers uk=0, and the unknown input is defined as(38)dk=0.4k,if3≤k<5,0.8,if5≤k<12,0,otherwise,and initial conditions x0=0.4−0.6T, d0=00T, x^0=00T, and d^0=00T. Figure 2 shows the state x2 being adequately estimated while in Figure 3, it can be seen that the observer effectively reconstructs the unknown input.
The next example provides an application of the UI observer hereby proposed for the task of actuator fault detection in a train system.Figure 1
Error signalek in the absence of unknown inputs.Figure 2
Statex2 and its estimation x^2.Figure 3
Unknown inputdk and its estimation d^k.Example 2.
This example provides a fault detection scheme via the proposed UI observer for a train system, and see Figure4 for a diagram. Let us consider first the fault-free case, and thus, a continuous-time model is [72–74](39)Mv˙t=bTvtuTt−bBvtuBt−wvt−gθp,p˙t=vt,where vt is the speed of the train, pt is the position of the train, bBv is the maximum coefficient of the braking force, uBt is the relative braking force, bTv is the maximum coefficient of the traction force, uTt is the relative traction force, gθ is the force of declivity or external force, wv is the friction force, and M is the mass of the train. In practical cases, the declivity force is considered as gθp=Mgsinθp, where θp is the slope angle on the position pt and g is the gravity force [75]. In [76], the friction force wv can be estimated via David’s equation:(40)wv=w0+w1v+w2v2,where w0, w1, and w2 are real coefficients that depend on the train characteristics and rail type; their values are unknown in the practice, but there exist some methods for their estimation [76]. By means of the Euler approximation x˙t≈xk+1−xk/Ts, where Ts is the sampling time. Therefore, a discrete-time train system of (23) is(41)Mvk+1=bTvkuTk−bBvkuBk−wvk−gθpkTs+Mvk,pk+1=vkTs+pk.
The fault studied in this example is the jamming fault, and it may occur only during the braking. This fault depends on the adherence conditions between the wheel and the track, producing that when the brake locks the wheel, the wheel slides on the track [68, 74]; in the real setup, the fault can occur depending on weather conditions. The main issue of a jamming fault is the impact directly on the measurements, producing a wrong estimate of both position and speed of the train. In Figure 5, it can be seen that the speed vk under the fault is lower, and this is because the wheel is locked during the fault and the sensor is not able to provide the right measure and thus both the speed and position are erroneous; for example, the measured position differs 4.10 meters from the real position, see Figure 6. Therefore, this fault can be considered as an actuator fault or exogenous input, when it occurs, the control is inhibited [72–74]. The actuator fault to be considered is(42)fauTk,uBk=−bTvkuTk−bBvkuBk+ξkfk,where fk∈0,1 represents the grade of the fault; if fk=0, then the system is fault-free; if fk=1, a total fault occurs and the control is inhibited. The term ξk∈ℝ represents a resistive force produced by the fault, and it is unknown [74].
As the previous examples, the unknown input is assumed to comply withdk+1=Sdk; therefore, we have an augmented system of form (3) with χ=vkpkdkdk+1T and matrices(43)E=M000010000100001,Ayk=−w1−w2vkTs+MTsgθpkpkTs0Ts100000100−12,ℬ=b000,C=1000T,and the parameter b is a constant value for bBv and bTv, producing a unique control signal uk, where uk>0 and uk<0 correspond to uTk and uBk, respectively. Let us consider the braking period, i.e., when the train is arriving to the station and there are no big slopes, we have gθpk=Mgsinθpk, with θpk=βpk∈−1,1, where β is a constant.
Thus, with UI observer (4), the error dynamics ek=χk−χ^k is(44)Eek+1=−w1vk−w2vk2+dkTs+Mvk−TsMgsinβpkTsvk+pkdk+1−dk+2dk+1−−w1v^k−w2v^k2+d^kTs+Mv^k−TsMgsinβp^kTsv^k+p^kd^k+1−d^k+2d^k+1−L⋅y−y^.
Applying the factorization in [41] and considering that yk=vk (the only available signal is the velocity vk), we have vk2−v^k2=vk+v^ke1=2yke1, e1=vk−v^k, as for sinβpk−sinβp^k, two terms of the series Taylor are taken; then, the error dynamics yields:(45)Eek+1=A¯yk,χk,χ^k−LykC¯ek,with matrices(46)A¯yk,χk,χ^k=−w1+2w2ykTs+MTsMg−1−β36pk2+pkp^k+p^k2Ts0Ts100000100−12,C¯=1000T.
Therefore, considering thatvk∈0,30m/s and pk,p^k∈0,600m, we have that the nonconstant terms are z=vk∈0,30 (available), ζ1=pk2∈0,360000, ζ2=pk∈0,600, ζ3=p^k∈0,600, and ζ4=p^k2∈0,360000 (nonavailable signals). The scheduling functions are defined as ωjζ=ωj11ζ1ωj22ζ2ωj33ζ3ωj44ζ4, ω0jζj=ζj1−ζj/ζj1−ζj0, ω1j=1−ω0jζj, and w1z=z1−z/z1−z0, w2z=1−w1z. The vertex matrices are(47)A¯i,j=−w1+2w2ziTs+MTsMg−1−β36ζ1j1+ζ2j2ζ3j3+ζ4j4Ts0Ts100001000−12,where i=1,2 and j=1,2,3,4. The LMI conditions in Theorem 1 have been found feasible with a decay rate α=0.7 (see Remark 1), and the parameters are M=408.5 tons, w0=3.525(N), w1=2.98×10−3(N/m), w2=4.575×10−4N2/m2, and Ts=0.001(s). The computed matrices are(48)N1=1×10−5139.522950.00157−0.000230.00003,N2=1×10−5151.497200.00172−0.000220.00003,G1=1×10−50.341520.00284−0.000700.000350.0000070.0000001−0.000000030.00000002−0.000001−0.000000030.00000001−0.000000010.00000080.00000002−0.000000010.000000008,G2=1×10−50.369920.00305−0.000710.000350.0000070.0000001−0.000000030.00000002−0.000001−0.000000030.00000001−0.000000010.00000080.00000002−0.000000010.000000008,P1=1×10−5125.554220.00267−0.000650.000330.002670.00000009−0.000000030.00000002−0.00065−0.000000030.00000001−0.000000010.000330.00000002−0.000000010.0000000075,P2=1×10−5129.303790.00278−0.000680.000350.002780.0000001−0.000000030.00000002−0.00068−0.000000030.00000001−0.000000010.000350.00000002−0.000000010.000000007.
A simulation has been performed for initial conditionsv0=28, p0=0, d0=00T, v^0=27.9, p^0=0.05, and d^0=0.050.01T, with a constant braking force uk=−0.5, β=0.0016, and the unknown input dk=fauk, with ξk=2M and the fault occurrence(49)fk=1,if4≤k≤4.2,1,if6≤k≤6.2,1,if8≤k≤8.2,0,otherwise.
As the previous example, first simulation has been performed without unknown inputs, i.e.,dk=0. Figure 7 shows the speed being adequately estimated while Figure 8 plots the Lyapunov function, whose signal is always positive and monotonously decreases to zero.
In the case wheredk=fauk, in order to detect the fault, the unknown input estimation is used as residual rk=d^k with a single threshold of ϵ=M/2. If the residual rk is bigger than ϵ, then the fault is occurring; otherwise, the fault has finished. Thus, we obtain a delay average of 0.003 seconds to detect when the fault begins and 0.006 to detect when the fault ends. As we can see from Figure 9, our proposal is enough to detect when the fault occurs.Figure 4
Fault detection scheme.Figure 5
Comparison between the time evolution of the speedvk with fault and without fault.Figure 6
Comparison between the time evolution of the positionpk with fault and without fault.Figure 7
The time evolution of the speedvk and its estimation v^k.Figure 8
The time evolution of the Lyapunov functionVe.Figure 9
The time evolution of the residualrk and the fault detection.


## 6. Conclusions

It has been presented a methodology for the design of unknown input observers for discrete-time nonlinear descriptor systems. It is based on algebraic rearrangements and allows overcoming the problem of unmeasurable premise variables as to get a compatible error dynamic system with the direct Lyapunov method; thus, sufficient conditions in terms of LMIs have been obtained. The resulting UI observer is proven to be less conservative than those in the literature. Moreover, the proposed scheme has been employed to solve the detection and estimation of actuator faults. The advantages of the proposal have been illustrated via numerical example and applied to the fault detection in train systems.



---

*Source: 8825609-2021-03-30.xml*