<figure><img src="image_1.png" /></figure>

HAL
open science

# Set-membership estimation for linear time-varying descriptor systems

Wentao Tang, Zhenhua Wang, Qinghua Zhang, Yi Shen

<figure><img src="image_2.png" /></figure>

## To cite this version:

Wentao Tang, Zhenhua Wang, Qinghua Zhang, Yi Shen. Set-membership estimation for linear time-varying descriptor systems. Automatica, 2020, 10.1016/j.automatica.2020.108867. hal-02970676

HAL Id: hal-02970676

https://inria.hal.science/hal-02970676

Submitted on 18 Oct 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Set-membership estimation for linear time-varying descriptor systems*

Wentao Tang<sup>a</sup>, Zhenhua Wang<sup>a</sup>, Qinghua Zhang<sup>b</sup>, Yi Shen<sup>a</sup>

<sup>a</sup>School of Astronautics, Harbin Institute of Technology, Harbin 150001, P. R. China

<sup>b</sup>Inria/IFSTTAR, Campus de Beaulieu, 35042, Rennes Cedex, France

## Abstract

This paper considers the problem of set-membership estimation for discrete-time linear time-varying descriptor systems subject to unknown but bounded disturbance and noise. We propose a set-membership estimation method based on a descriptor system observer and a zonotopic estimator of the observer error bounds. The observer parameters are optimized in order to minimize the sizes of the zonotopes enclosing all admissible state trajectories. Finally, two simulation results are provided to demonstrate the effectiveness of the proposed method.

**Key words:** Set-membership estimation; Linear time-varying; Descriptor systems; Zonotopes.

## 1 Introduction

State estimation is important in control theory and applications. However, the ubiquitous uncertainties influence the performance of state estimation. Many existing state estimation methods assume that the uncertainties have known probability distribution functions, which may not be true in practice (Alamo, Bravo & Camacho, 2005). Based on a general assumption that the uncertainties are unknown but bounded, set-membership estimation requires no information about the probability distributions of the uncertainties. Set-membership estimation methods have been applied to various fields such as bioprocess monitoring (Gouzé, Rapaport & Hadj-Sadok, 2000; Moisan, Bernard & Gouzé, 2009), robust control (Canale, Fagiano & Milanese, 2009; Efimov, Raïssi & Zolghadri, 2013b; Zhang, Yang, Han & Zhu, 2019) and fault diagnosis (Xu, Tan, Wang, Wang, Liang & Yuan, 2019). Different kinds of geometrical sets such as ellipsoids (Fogel & Huang, 1982; Chernousko, 2005), polytopes (Shamma & Tu, 1999; Blanchini &

Miani, 2008) and intervals (Raïssi, Efimov & Zolghadri, 2012; Efimov, Raïssi, Chebotarev & Zolghadri, 2013a; Wang, Lim & Shen, 2018b) have been used to design set-membership estimation methods. Among these sets, zonotopes can achieve a good tradeoff between estimation accuracy and computation complexity. Recently, the set-membership estimation methods based on zonotopes have received much attention (Le, Stoica, Alamo, Camacho & Dumur, 2013; Combastel, 2015; Scott, Raimondo, Marseglia & Braatz, 2016; Tang, Wang, Wang, Raïssi & Shen, 2019).

Most works on set-membership estimation focus on regular systems. Descriptor systems are more general and can describe many practical systems. Traditional state estimation methods for descriptor systems have been studied extensively. However, only few results study set-membership estimation for descriptor systems. Wang, Puig & Cembrano (2018a) and Wang, Olaru, Valmorobida, Puig & Cembrano (2019a) studied set-membership estimation methods for descriptor systems based on zonotopes. However, Wang et al. (2018a) and Wang et al. (2019a) only considered time-invariant systems.

In this paper, we propose a new set-membership estimation method for discrete-time linear time-varying (LTV) descriptor systems subject to unknown but bounded disturbances and measurement noises. The main contributions of this paper reside in the following aspects:

<footnote>* This work was supported by National Natural Science Foundation of China (Grant No. 61973098, 61773145) and the Key Laboratory Opening Funds of Harbin Institute of Technology under grant HIT.KLOF.2018.073. This paper was not presented at any IFAC meeting. Corresponding author: Zhenhua Wang. Tel. +86-451-86413411-8602.

Email addresses: zhenhua.wang@hit.edu.cn (Zhenhua Wang), qinghua.zhang@inria.fr (Qinghua Zhang).</footnote>
1) A novel set-membership estimation method for discrete-time LTV descriptor systems based on a simple observer structure, which is easy to design and implement.

2) The optimization of all the parameters of the proposed set-membership estimator, with a closed-form solution.

3) A rigorous proof of the optimality minimizing the zonotope size criterion, not simply relying on the first order condition.

Set-membership estimation for LTV descriptor systems has been studied in Wang, Wang, Puig & Cembrano (2019b), but the estimator presented in Wang et al. (2019b) involves coupled parameters, which are only partly optimized. It is also limited to descriptor systems with time-invariant output equations. In contrast, the proposed method is applicable to fully LTV descriptor systems and can optimize all the estimator parameters. Moreover, sufficient optimality conditions are established for the proposed estimator, whereas in some similar results, only the necessary first order condition is considered, like in Combastel (2015).

**Notations:** $\mathbb{R}^n$ and $\mathbb{R}^{m \times n}$ denote the $n$ and $m \times n$ dimensional Euclidean spaces, respectively. $I_n$ denotes the identity matrix with dimensions of $n \times n$. For a matrix $A \in \mathbb{R}^{m \times n}$, $A^T$ denotes its transpose and $A^\dagger$ represents its pseudo-inverse. $\text{vec}(A)$ denotes the vectorization of $A$ by concatenating the columns of $A$ into a single column vector. $\sigma_{\min}(A)$ denotes the smallest singular value of $A$. For a square matrix $A$, $\text{tr}(A)$ is its trace. For two matrices $X$ and $Y$, $X \otimes Y$ denotes their Kronecker product. In this paper, the symbols $\geq$, $>$, $\leq$, $<$ and the absolute value operator $| \cdot |$ should be understood elementwise. $P > 0$ ($P < 0$) indicates that $P$ is a positive (negative) definite matrix. $P \geq 0$ ($P \leq 0$) indicates that $P$ is a positive (negative) semidefinite matrix.

## 2 Preliminaries

**Definition 1** An $m$-order zonotope $\mathcal{Z} \subset \mathbb{R}^n$ is the affine transformation of a hypercube $\mathbf{B}^m = [-1, 1]^m$:

$$ \mathcal{Z} = p \oplus H\mathbf{B}^m = \{p + Hz : z \in \mathbf{B}^m\}, $$

where $\oplus$ denotes the Minkowski sum, $p \in \mathbb{R}^n$ is the center of $\mathcal{Z}$, and $H \in \mathbb{R}^{n \times m}$ is called the generator matrix of $\mathcal{Z}$. For simplicity, we also denote $\mathcal{Z} = \langle p, H \rangle$.

**Property 1** (Combastel, 2015) For zonotopes, the following properties hold:

$$ \langle p_1, H_1 \rangle \oplus \langle p_2, H_2 \rangle = \langle p_1 + p_2, [H_1 \quad H_2] \rangle, \quad (1a) $$

$$ L\langle p, H \rangle = \langle Lp, LH \rangle, \quad (1b) $$

$$ \langle p, H \rangle \subseteq \langle p, \text{rs}(H) \rangle, \quad (1c) $$

where $p, p_1, p_2 \in \mathbb{R}^n$, $H, H_1, H_2 \in \mathbb{R}^{n \times m}$ and $L \in \mathbb{R}^{l \times n}$. $\text{rs}(H) = \text{diag}([h_1, \dots, h_n])$ where $h_i = \sum_{j=1}^m |H_{i,j}|$, $i = 1, \dots, n$. Note that $\langle p, \text{rs}(H) \rangle$ is also a box, i.e., an interval vector, which represents a box outer bound of $\langle p, H \rangle$.

A high-order zonotope can be bounded by a lower one via the reduction operation based on the following property.

**Property 2** (Alamo et al., 2005) Given the zonotope $\mathcal{Z} = p \oplus H\mathbf{B}^q \subset \mathbb{R}^n$ and a chosen integer $s$ ($n < s < q$), denote the matrix obtained by reordering the columns of $H$ in decreasing Euclidean norm as $\hat{H}$. Then $\mathcal{Z} \subset \langle p, \mathcal{R}_s(H) \rangle$, where $\mathcal{R}_s(H) = [H^a \quad H^b]$, $H^a$ is composed of the first $s-n$ columns of $\hat{H}$ and $H^b$ is a diagonal matrix with

$$ H_{i,i}^b = \sum_{j=s-n+1}^{q} |\hat{H}_{i,j}|, \quad i = 1, \dots, n. $$

**Property 3** (Golub & Van Loan, 1996) For matrices $X$, $A$, $B$, $C$ with appropriate dimensions, the following equations hold:

$$ \text{tr}(A) = \text{tr}(A^T), \quad (2a) $$

$$ \frac{\partial}{\partial X} \text{tr}(AX^T B) = BA, \quad (2b) $$

$$ \frac{\partial}{\partial X} \text{tr}(AXB^T C) = A^T C^T X B^T + C A X B. \quad (2c) $$

**Lemma 1** (Wang et al., 2018b) Given matrices $\mathcal{W} \in \mathbb{R}^{p \times n}$ and $\mathcal{Y} \in \mathbb{R}^{m \times n}$, if $\text{rank}(\mathcal{W}) = n$, then the general solution to $\mathcal{X}\mathcal{W} = \mathcal{Y}$ exists and is given by

$$ \mathcal{X} = \mathcal{Y}\mathcal{W}^\dagger + S(I_p - \mathcal{W}\mathcal{W}^\dagger) $$

where $S \in \mathbb{R}^{m \times p}$ is an arbitrary matrix.

## 3 Problem statement

Consider the following system

$$ \left\{ \begin{array}{l} E_k x_k = A_k x_{k-1} + B_k u_{k-1} + w_{k-1} \\ y_k = C_k x_k + v_k \end{array} \right. \quad (3) $$

where $x_k \in \mathbb{R}^{n_x}$, $u_{k-1} \in \mathbb{R}^{n_u}$ and $y_k \in \mathbb{R}^{n_y}$ are the vectors of state, input and measurement output, respectively. $E_k \in \mathbb{R}^{n_x \times n_x}$, $A_k \in \mathbb{R}^{n_x \times n_x}$, $B_k \in \mathbb{R}^{n_x \times n_u}$ and $C_k \in \mathbb{R}^{n_y \times n_x}$ are known matrices. $E_k$ may be singular. $w_{k-1} \in \mathbb{R}^{n_x}$ and $v_k \in \mathbb{R}^{n_y}$ denote the unknown disturbance and the measurement noise, respectively.

For system (3), we have the following assumptions.
**Assumption 1.** $x_0$, $w_k$ and $v_k$ are unknown but bounded as follows:

$$x_0 \in \langle p_0, H_0 \rangle, \quad w_k \in \langle 0, D_k \rangle, \quad v_k \in \langle 0, F_k \rangle, \quad \forall k \ge 0,$$

where $p_0 \in \mathbb{R}^{n_x}$, $H_0 \in \mathbb{R}^{n_x \times n_x}$, $D_k \in \mathbb{R}^{n_x \times n_x}$ and $F_k \in \mathbb{R}^{n_y \times n_y}$ are known vector and matrices.

**Assumption 2.** We assume that

$$\sigma_{\min}(D_k) \ge \epsilon_1, \quad \sigma_{\min}(F_k) \ge \epsilon_2, \quad \forall k \ge 0 \quad (4)$$

where $\epsilon_1$ and $\epsilon_2$ are positive constants.

**Assumption 3.** The parameters of system (3) are assumed to satisfy the following condition:

$$\mathrm{rank} \begin{bmatrix} E_k \\ C_k \end{bmatrix} = n_x, \quad \forall k \ge 0. \quad (5)$$

**Remark 1.** Assumption 3 is commonly used in descriptor system estimation (e.g. Nikoukhah, Willsky & Levy, 1990; Hamdi, Rodrigues, Mechmeche & Braiek, 2012; Zhang, Swain & Nguang, 2014).

In this paper, we aim to develop an optimal set-membership estimation method for (3) to obtain the minimized zonotope sets that can tightly enclose all admissible state trajectories. In addition, we can further obtain the upper and lower bounds of $x_k$, i.e., an interval vector $[\underline{x}_k, \bar{x}_k]$ satisfying $\underline{x}_k \le x_k \le \bar{x}_k$.

# 4 Main results

## 4.1 Set-membership estimator

We propose a set-membership estimator for system (3) to obtain zonotopes that can enclose possible state trajectories. The proposed set-membership estimator will consist of two parts: a single state trajectory estimator and a bound estimator of the estimation error of this trajectory estimator. The first part is based on the observer structure proposed in Wang, Shen, Zhang & Wang (2012), which is in the form of

$$\hat{x}_k = T_k A_k \hat{x}_{k-1} + T_k B_k u_{k-1} + N_k y_k + L_k (y_{k-1} - C_{k-1} \hat{x}_{k-1}), \quad (6)$$

where $\hat{x}_k \in \mathbb{R}^{n_x}$ is the single trajectory estimation, $T_k \in \mathbb{R}^{n_x \times n_x}$, $N_k \in \mathbb{R}^{n_x \times n_y}$ and $L_k \in \mathbb{R}^{n_x \times n_y}$ are the matrices to be designed. $T_k$ and $N_k$ should satisfy

$$T_k E_k + N_k C_k = I_{n_x}. \quad (7)$$

**Remark 2.** According to Lemma 1, if Assumption 3 holds, then there exist $T_k$ and $N_k$ such that (7) holds.

We use observer (6) to generate the single trajectory estimation. Then, the following theorem is used to obtain the set-membership estimation.

**Theorem 1** Given matrices $T_k$, $N_k$ and $L_k$, if $x_{k-1} \in \langle \hat{x}_{k-1}, H_{k-1} \rangle \subset \mathbb{R}^{n_x}$, then $x_k \in \langle \hat{x}_k, H_k \rangle$, where $\hat{x}_k$ is obtained from (6) and $H_k$ satisfies

$$H_k = [(T_k A_k - L_k C_{k-1}) \tilde{H}_{k-1} G_k], \quad (8)$$

where

$$G_k = \begin{bmatrix} T_k D_{k-1} - L_k F_{k-1} - N_k F_k \\ \tilde{H}_{k-1} = \begin{cases} \mathcal{R}_s(H_{k-1}), & q > s; \\ H_{k-1}, & q \le s. \end{cases} \end{bmatrix}, \quad (9)$$

Herein, $q$ is the number of columns of $H_{k-1}$ and $s$ is a fixed integer.

**PROOF.** By combining (3) with (7), we have

$$x_k = T_k A_k x_{k-1} + T_k B_k u_{k-1} + T_k w_{k-1} + N_k y_k - N_k v_k. \quad (10)$$

Define the estimation error as $e_k = x_k - \hat{x}_k$. Since $x_{k-1} \in \langle \hat{x}_{k-1}, H_{k-1} \rangle$, it follows that $e_{k-1} \in \langle 0, H_{k-1} \rangle$.

According to Property 2, (9) implies $\langle 0, H_{k-1} \rangle \subseteq \langle 0, \tilde{H}_{k-1} \rangle$. Then, we have $e_{k-1} \in \langle 0, \tilde{H}_{k-1} \rangle$.

Subtracting (6) from (10) yields

$$e_k = (T_k A_k - L_k C_{k-1}) e_{k-1} + T_k w_{k-1} - L_k v_{k-1} - N_k v_k. \quad (11)$$

According to Assumption 1, (11) implies

$$e_k \in (T_k A_k - L_k C_{k-1}) \langle 0, \tilde{H}_{k-1} \rangle \oplus T_k \langle 0, D_{k-1} \rangle \oplus (-L_k \langle 0, F_{k-1} \rangle) \oplus (-N_k \langle 0, F_k \rangle)$$

Then, by using Property 1, we have $e_k \in \langle 0, H_k \rangle$.

Since $x_k = \hat{x}_k + e_k$, it follows that $x_k \in \langle \hat{x}_k, H_k \rangle$. $\square$

**Remark 3.** According to the property in (1c), we can obtain the upper and lower bounds of $x_k$ as follows:

$$\begin{cases} \bar{x}_k = \hat{x}_k + \mathrm{rs}(H_k), \\ \underline{x}_k = \hat{x}_k - \mathrm{rs}(H_k). \end{cases}$$
## 4.2 Optimization design

To obtain accurate estimation, the matrices $T_k, N_k$ and $L_k$ should be designed to minimize the size of the constructed zonotope set, $\langle \hat{x}_k, H_k \rangle$. In this paper, we chose

$$J_k(T_k, N_k, L_k) = \|H_k(T_k, N_k, L_k)\|_F^2 = \text{tr}(H_k^T H_k) \quad (12)$$

as the size criterion of $\langle \hat{x}_k, H_k \rangle$. It can be interpreted as the size of the zonotope generator segments, as mentioned in Alamo et al. (2005) and Combastel (2015), or more precisely, the sum of the squares of the generator segment half lengths. It is chosen in this paper due to its simplicity in analysis and implementation.

From (8), $J_k$ is a quadratic form in terms of $T_k, N_k$ and $L_k$. In addition, $T_k$ and $N_k$ are required to satisfy the equality constraint (7). Therefore, we aim to solve a constrained quadratic optimization problem as follows:

$$\min_{T_k E_k + N_k C_k = I_{n_x}} J_k(T_k, N_k, L_k). \quad (13)$$

According to Lemma 1, $T_k$ and $N_k$ can be chosen as

$$\begin{aligned} T_k &= \Theta_k^\dagger \alpha_1 + S_k \Omega_k \alpha_1, \\ N_k &= \Theta_k^\dagger \alpha_2 + S_k \Omega_k \alpha_2, \end{aligned} \quad (14)$$

where

$$\begin{aligned} \Theta_k &= \begin{bmatrix} E_k \\ C_k \end{bmatrix}, \quad \alpha_1 = \begin{bmatrix} I_{n_x} \\ 0 \end{bmatrix}, \quad \alpha_2 = \begin{bmatrix} 0 \\ I_{n_y} \end{bmatrix}, \\ \Omega_k &= I_{n_x+n_y} - \Theta_k \Theta_k^\dagger. \end{aligned}$$

Note that it is difficult to minimize the size criterion $J_k$ in terms of $S_k$, because the matrix $\Omega_k$ is always rank deficient. In Wang et al. (2019b), $S_k$ is chosen based on personal experience. In fact, if $S_k$ is given, we can obtain the optimal $L_k$ to minimize the size criterion.

**Theorem 2** Given any matrix pair $T_k$ and $N_k$ of appropriate dimensions satisfying (7) (possibly chosen following (14)), then the $L_k$ obtained from (15) is the optimal one minimizing the size criterion $J_k$.

$$\begin{aligned} V_{k-1} &= C_{k-1} P_{k-1} C_{k-1}^T + R_{k-1}, \\ L_k &= T_k A_k P_{k-1} C_{k-1}^T V_{k-1}^{-1}, \end{aligned} \quad (15)$$

where

$$P_{k-1} = \tilde{H}_{k-1} \tilde{H}_{k-1}^T, \quad R_{k-1} = F_{k-1} F_{k-1}^T. \quad (16)$$

## PROOF.

$$\begin{aligned} J_k = \text{tr}((T_k A_k - L_k C_{k-1}) P_{k-1} (T_k A_k - L_k C_{k-1})^T) \\ + \text{tr}(T_k Q_{k-1} T_k^T) + \text{tr}(L_k R_{k-1} L_k^T) \\ + \text{tr}(N_k R_k N_k^T), \end{aligned}$$

where $Q_{k-1} = D_{k-1} D_{k-1}^T$. By using the matrix calculus properties in (2), we can obtain the partial derivative of $J_k$ for $L_k$ as follows:

$$\frac{\partial J_k}{\partial L_k} = 2L_k(C_{k-1}P_{k-1}C_{k-1}^T + R_{k-1}) - 2T_k A_k P_{k-1}C_{k-1}^T.$$

The $L_k$ obtained from (15) satisfies $\frac{\partial J}{\partial L_k} = 0$, which means that it is a stationary point.

From (4) and (16), we have $R_{k-1} \succ 0$. Since $V_{k-1} \succeq R_{k-1}$, it follows that $V_{k-1}$ is invertible. Therefore, the $L_k$ obtained from (15) is the unique stationary point of $J_k$. In addition, $J_k$ is a quadratic function of $L_k$ and by definition $J_k \ge 0$ for any $L_k$, then $J_k$ is a convex function of $L_k$. Therefore, the $L_k$ obtained from (15) is the unique minimum point of $J_k$. □

Although $T_k$ and $N_k$ can be obtained from (14), it is difficult to find the optimal $S_k$ minimizing the size criterion $J_k$. A freely chosen $S_k$ may cause large conservatism.

Instead of optimizing $T_k$ and $N_k$ with the parametrization in terms of $S_k$, we propose the following theorem directly optimizing $T_k, N_k$ and $L_k$ to minimize the size criterion $J_k$, under the equality constraint (7).

**Theorem 3** The constrained optimization problem (13) has a unique optimum in terms of $T_k, N_k$ and $L_k$, which is explicitly obtained from

$$V_{k-1} = C_{k-1} P_{k-1} C_{k-1}^T + R_{k-1}, \quad (17a)$$

$$M_k = A_k P_{k-1} C_{k-1}^T V_{k-1}^{-1} C_{k-1} P_{k-1} A_k^T, \quad (17b)$$

$$W_k = A_k P_{k-1} A_k^T + Q_{k-1}, \quad (17c)$$

$$O_k = 2(M_k - W_k), \quad (17d)$$

$$\Psi_k = -2R_k, \quad (17e)$$

$$\Lambda_k = (E_k^T O_k^{-1} E_k + C_k^T \Psi_k^{-1} C_k)^{-1}, \quad (17f)$$

$$T_k = \Lambda_k E_k^T O_k^{-1}, \quad (17g)$$

$$N_k = \Lambda_k C_k^T \Psi_k^{-1}, \quad (17h)$$

$$L_k = T_k A_k P_{k-1} C_{k-1}^T V_{k-1}^{-1}. \quad (17i)$$

## PROOF.
To prove that the solution $(T_k, N_k, L_k)$ obtained from (17) is the unique optimum of the constrained optimization problem (13), usually the first order and the second order optimality conditions should be checked. However, this is a non-trivial task due to the fact that the constraint (7) and the size criterion
$J_k$ are formulated in terms of parameter matrices, instead of the vector form in standard quadratic optimization problems. To overcome this difficulty, we conduct our proof in two steps. First, we prove that the solution ($T_k, N_k, L_k, \Lambda_k$) obtained from (17) is the unique stationary point of the Lagrange function of the constrained optimization problem. Second, we verify that the constrained optimization problem (13) satisfies the condition ensuring that the unique stationary point is the unique optimum of the considered problem.

First, define the following Lagrange function:

$$Z_k(T_k, N_k, L_k, \Lambda_k) = J_k(T_k, N_k, L_k) + \text{tr}(\Lambda_k \Phi_k),$$

where $\Lambda_k \in \mathbb{R}^{n_x \times n_x}$ and $\Phi_k = E_k^T T_k^T + C_k^T N_k^T - I_{n_x}$. According to the matrix calculus properties in (2), the partial derivatives of $Z_k$ for $T_k, N_k, L_k$ and $\Lambda_k$ are

$$\frac{\partial Z_k}{\partial T_k} = 2T_k W_k - 2L_k C_{k-1} P_{k-1} A_k^T + \Lambda_k E_k^T,$$

$$\frac{\partial Z_k}{\partial N_k} = 2N_k R_k + \Lambda_k C_k^T,$$

$$\frac{\partial Z_k}{\partial L_k} = 2L_k V_{k-1} - 2T_k A_k P_{k-1} C_{k-1}^T,$$

$$\frac{\partial Z_k}{\partial \Lambda_k} = T_k E_k + N_k C_k - I_{n_x}.$$

Let all the partial derivatives equal to 0, we have:

$$\left\{ \begin{array}{l} X_k \Sigma_k + \Lambda_k \Pi_k^T = 0, \\ X_k \Pi_k = I_{n_x}, \end{array} \right. \qquad (18a)$$

where

$$X_k = [T_k \ N_k \ L_k], \quad \Pi_k = \begin{bmatrix} E_k \\ C_k \\ 0 \end{bmatrix}, \qquad (19) \\ \Sigma_k = 2 \begin{bmatrix} W_k & 0 & -A_k P_{k-1} C_{k-1}^T \\ 0 & R_k & 0 \\ -C_{k-1} P_{k-1} A_k^T & 0 & V_{k-1} \end{bmatrix}.$$

From (17a) and (17c), the matrix $\Sigma_k$ can be rewritten as

$$\Sigma_k = \Xi_k + \Gamma_k, \qquad (20)$$

where

$$\Xi_k = 2 \begin{bmatrix} Q_{k-1} & 0 & 0 \\ 0 & R_k & 0 \\ 0 & 0 & R_{k-1} \end{bmatrix}, \\ \Gamma_k = 2 \begin{bmatrix} A_k \\ 0 \\ -C_{k-1} \end{bmatrix} P_{k-1} \begin{bmatrix} A_k \\ 0 \\ -C_{k-1} \end{bmatrix}^T.$$

From (4), we have $Q_k \succ 0$ and $R_k \succ 0$ for any $k \ge 0$. It follows that $\Xi_k \succ 0$. In addition, $\Gamma_k \succeq 0$ since $P_{k-1} \succeq 0$. Therefore, from (20), we have $\Sigma_k \succ 0$, which implies that $\Sigma_k$ is invertible. Then from (18a), we have

$$X_k = -\Lambda_k \Pi_k^T \Sigma_k^{-1}. \qquad (21)$$

Substituting it into (18b) yields

$$\Lambda_k(-\Pi_k^T \Sigma_k^{-1} \Pi_k) = I_{n_x}. \qquad (22)$$

From (5), we have rank $\Pi_k = n_x$. Then, $\Pi^T$ is row full rank. Since $\Sigma_k^{-1}$ is invertible, it follows that $-\Pi_k^T \Sigma_k^{-1} \Pi_k$ is invertible. Therefore, $\Lambda_k$ can be uniquely determined by (22). Then, $X_k = [T_k \ N_k \ L_k]$ is uniquely determined by (21). Therefore, (18) has one unique solution. It is easy to check that ($T_k, N_k, L_k, \Lambda_k$) obtained from (17) is the solution to (18). Therefore, ($T_k, N_k, L_k, \Lambda_k$) obtained from (17) is the unique stationary point of the Lagrange function.

Next, we will verify that the unique stationary point is the unique optimum of the considered problem.

Define $x = \text{vec}(X_k)$. The size criterion $J_k$ defined in (12) is a quadratic form in $x$ (a homogeneous polynomial containing only second degree terms of $x$). Moreover, this quadratic form is positive semidefinite, since it can never be negative for any value of $x$.

From the constraint (7), we have $X_k \Pi_k = I_{n_x}$, which is equivalent to

$$A x = \text{vec}(I_{n_x}),$$

where $A = \Pi_k^T \otimes I_{n_x}$. Since $\Pi_k^T$ is row full rank, it follows that $A$ is row full rank. Then, according to Markowitz (1956), the solution obtained from (17) is the unique optimum of the considered problem. $\square$

**Remark 4.** Assumption 2 has been introduced to establish, in Theorem 3, the optimality of the proposed solution. This assumption can be relaxed if the optimal solution is simply deduced by the first order optimality condition, like in Combastel (2015), which is not sufficient to ensure the optimality. In fact, the optimality
proof of Theorem 3 requires that the matrix $\Sigma_k$ is invertible, which is a weaker condition than the conditions formulated in (4). However, in practice, it is more difficult to check the large matrix $\Sigma_k$ defined in (19), and easier to check the matrices $D_k$ and $F_k$. Assumption 2 means that uncertainties exist in all the directions of the state space and of the output space. In practice, if there is no uncertainty in some directions, a simple solution is to slightly modify $D_k$ or $F_k$ so that (4) is satisfied with a small value of $\epsilon_1$ or $\epsilon_2$. This approach is often used for parameter estimation with methods related to the Kalman filter. Alternatively, for a more accurate solution, the dimension of the state equation or of the output equation could be reduced so that Assumption 2 becomes true, but this approach still requires non-trivial work in future studies.

## 5 Simulation results

A numerical example is used to demonstrate the performance of the proposed method, which has parameters as follows:

$$ 
\begin{aligned}
E_k &= \begin{bmatrix} 1 + 0.2 \sin(0.1k) & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \\
A_k &= \begin{bmatrix} 0.6 - 0.1 \sin(0.1k) & 0 & 0.3 \\ -0.2 \sin(0.2k) & 0.4 & 0 \\ 0 & 1 & 1 \end{bmatrix}, \\
B_k &= \begin{bmatrix} 1 \\ \sin(0.1k) \\ 0 \end{bmatrix}, \\
C_k &= \begin{bmatrix} 0 & 1 + 0.2 \sin(0.2k) & 0.5 \end{bmatrix}.
\end{aligned}
$$

In the simulation study, the disturbance $w_k$ and measurement noise $v_k$ are bounded as $w_k \in \langle 0, 0.04I_3 \rangle$ and $v_k \in \langle 0, 0.01 \rangle$. They are randomly generated with uniform distributions corresponding to the assumed uncertainty sets. The initial state $x_0 = [0 \ 0 \ -0.0301]^T$. We set $\hat{x}_0 = [0 \ 0 \ 0]^T$ and $H_0 = 0.2I_3$. The input is set as $u_k = 0.5 \sin(0.15k)$. The maximal zonotope order is chosen as 20. The methods based on Theorem 2 and Theorem 3 are compared in the simulation. For the method based on Theorem 2, we set $S_k = 0$.

The simulation results are shown in Fig. 1. Thanks to sufficient optimization, the estimation results obtained by the method based on Theorem 3 are more accurate than those by the method based on Theorem 2.

To better demonstrate the improvement of Theorem 3, some quantitative results are provided in Table 1, which shows that the method based on Theorem 3 can obtain more accurate estimation bounds for each state.

<figure><img src="image_3.png" /><figcaption>Fig. 1. The interval estimation results of state.</figcaption></figure>

<table><caption>Table 1
The average estimation bounds for each state</caption><thead><tr><th>Method</th><th>x<sub>k</sub>(1)</th><th>x<sub>k</sub>(2)</th><th>x<sub>k</sub>(3)</th></tr></thead><tbody><tr><td>Theorem 2</td><td>0.5778</td><td>0.2440</td><td>0.5498</td></tr><tr><td>Theorem 3</td><td>0.3508</td><td>0.1786</td><td>0.4024</td></tr></tbody></table>

To more explicitly demonstrate the improvement brought by the sufficient optimization in Theorem 3, the estimation zonotope sets in a few time instants are depicted in Fig. 2, where the green zonotopes are the estimation results obtained by the method based on Theorem 2 and the red ones are those by the method based on Theorem 3. At instant $k = 0$, the two initial zonotopes are identical. They are drawn in gray instead of green or red. Fig. 2 also shows that the method based on Theorem 3 can obtain more accurate estimation results than those by the method based on Theorem 2.

<figure><img src="image_4.png" /><figcaption>Fig. 2. The estimated zonotopes in a few time instants.</figcaption></figure>

Note that the method in Wang et al. (2019b) cannot apply to the above example since the measurement matrix is time-varying. To better demonstrate the performance of the proposed method based on Theorem 3, an example
ple of a truck-trailer system from Wang et al. (2019b) is used to compare the proposed method with the method in Wang et al. (2019b). The system example in Wang et al. (2019b) is in a polytopic LPV form. It can be reformulated into an equivalent form which is identical to (3). The parameters of the equivalent form are as follows:

$$ 
\begin{aligned}
E &= \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & -1 \\ 0 & 0 & 0 & 0 \end{bmatrix}, \quad A_k = \begin{bmatrix} 1 - \frac{\theta_k T_s}{L} & 0 & 0 & 0 \\ \frac{\theta_k T_s}{L} & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & -\frac{\theta_k T_s}{2} & 0 & 1 \end{bmatrix}, \\
B_k &= \begin{bmatrix} \frac{\theta_k T_s}{l} \\ 0 \\ 0 \\ 0 \end{bmatrix}, \quad C_k = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix},
\end{aligned}
$$

where $\theta_k = -0.9 + 0.3 \sin(0.0386k + 1)$. And $L = 5.5m$ is the length of trailer, $l = 4.8m$ is the length of truck, $T_s = 0.2s$ is the sampling time.

In the simulation, $x_0 = [-0.1745 \ 0.0873 \ 3.0 \ -0.0057]^T$. We set $H_0 = \text{diag}([0.1 \ 0.02 \ 0.1 \ 0.02])$ and $\hat{x}_0 = [-0.1 \ 0.08 \ 3 \ 0]^T$, then $x_0 \in \langle \hat{x}_0, H_0 \rangle$. The disturbance and noise are bounded as $w_k \in \langle 0, D_k \rangle$ and $v_k \in \langle 0, 0.002I_3 \rangle$, where $D_k = \text{diag}([0 \ 0 \ 0.006\theta_k \ 0])$.

The maximal zonotope order is set as 20, which is the same with that in Wang et al. (2019b). $w_k$ and $v_k$ are randomly generated with uniform distributions corresponding to the assumed uncertainty sets. The simulation results are given in Fig. 3, which shows that the estimation results of all state components obtained by the proposed method are more accurate than those by the method in Wang et al. (2019b). Some quantitative results of the accuracy of the two compared methods are given in Table 2, which also illustrates the superiority of the proposed method.

<table><caption>Table 2
The average estimation bounds for each state</caption><thead><tr><th>Method</th><th>x<sub>k</sub>(1)</th><th>x<sub>k</sub>(2)</th><th>x<sub>k</sub>(3)</th><th>x<sub>k</sub>(4)</th></tr></thead><tbody><tr><td>Wang et al. (2019b)</td><td>0.6606</td><td>0.1046</td><td>0.0953</td><td>0.0397</td></tr><tr><td>The proposed one</td><td>0.1331</td><td>0.0084</td><td>0.0063</td><td>0.0057</td></tr></tbody></table>

Note that $D_k$ does not satisfy (4), which has been assumed for Theorem 3. Nevertheless, the matrix $\Sigma_k$ is invertible for all $k$, as shown by its condition number monitored during the simulation, and the optimality result of Theorem 3 holds also under this weaker condition. See Remark 4 in the previous section. The example with a singular $D_k$ is for the purpose of comparison

<figure><img src="image_5.png" /><figcaption>Fig. 3. The estimation results of the two compared methods.</figcaption></figure>

with Wang et al. (2019b), where this example was initially presented. In practice, it is less easy to check the condition on the large matrix $\Sigma_k$ than the stronger conditions (4), which may be weakened in future studies. To simply rely on the conditions expressed in (4), simulations have also been made after modifications of $D_k$ satisfying (4), leading to visually unnoticeable changes in the numerical results.

## 6 Conclusion

In this paper, we propose a zonotopic set-membership estimation method for discrete-time LTV descriptor systems. The proposed set-membership estimator is based on combining a single state trajectory estimator with the reachable set of estimation error. We propose an optimization method to design the parameters of the estimator such that the considered size criterion is minimized. The simulation examples have shown the effectiveness of the proposed method.

## References

Alamo, T., Bravo, J. M., & Camacho, E. F. (2005). Guaranteed state estimation by zonotopes. Automatica, 41, 1035–1043.

Blanchini, F., & Miani, S. (2008). Set-theoretic methods in control. Springer.

Canale, M., Fagiano, L., & Milanese, M. (2009). Set membership approximation theory for fast implementation of model predictive control laws. Automatica, 45, 45–54.

Chernousko, F. L. (2005). Ellipsoidal state estimation for dynamical systems. Nonlinear Analysis: Theory, Methods & Applications, 63, 872–879.

Combastel, C. (2015). Zonotopes and Kalman observers: gain optimality under distinct uncertainty paradigms and robust convergence. Automatica, 55, 265–273.
Efimov, D., Raïssi, T., Chebotarev, S., & Zolghadri, A. (2013a). Interval state observer for nonlinear time varying systems. *Automatica*, **49**, 200–205.

Efimov, D., Raïssi, T., & Zolghadri, A. (2013b). Control of nonlinear and LPV systems: interval observer-based framework. *IEEE Transactions on Automatic Control*, **58**, 773–778.

Fogel, E., & Huang, Y. F. (1982). On the value of information in system identification—bounded noise case. *Automatica*, **18**, 229–238.

Golub, G. H., & Van Loan, C. F. (1996). *Matrix computations*. Baltimore, MD, USA: Johns Hopkins University Press.

Gouzé, J. L., Rapaport, A., & Hadj-Sadok, M. Z. (2000). Interval observers for uncertain biological systems. *Ecological modelling*, **133**, 45–56.

Hamdi, H., Rodrigues, M., Mechmeche, C., & Braiek, N. B. (2012). Robust fault detection and estimation for descriptor systems based on multi-models concept. *International Journal of Control, Automation and Systems*, **10**, 1260–1266.

Le, V. T. H., Stoica, C., Alamo, T., Camacho, E. F., & Dumur, D. (2013). Zonotopic guaranteed state estimation for uncertain systems. *Automatica*, **49**, 3418–3424.

Markowitz, H. (1956). The optimization of a quadratic function subject to linear constraints. *Naval Research Logistics Quarterly*, **3**, 111–133.

Moisan, M., Bernard, O., & Gouzé, J. L. (2009). Near optimal interval observers bundle for uncertain bioreactors. *Automatica*, **45**, 291–295.

Nikoukhah, R., Willsky, A. S., & Levy, B. C. (1990). Kalman filtering and riccati equations for descriptor systems. In *29th IEEE Conference on Decision and Control* (pp. 2886–2891).

Raïssi, T., Efimov, D., & Zolghadri, A. (2012). Interval state estimation for a class of nonlinear systems. *IEEE Transactions on Automatic Control*, **57**, 260–265.

Scott, J. K., Raimondo, D. M., Marseglia, G. R., & Braatz, R. D. (2016). Constrained zonotopes: A new tool for set-based estimation and fault detection. *Automatica*, **69**, 126–136.

Shamma, J. S., & Tu, K. Y. (1999). Set-valued observers and optimal disturbance rejection. *IEEE Transactions on Automatic Control*, **44**, 253–264.

Tang, W., Wang, Z., Wang, Y., Raïssi, T., & Shen, Y. (2019). Interval estimation methods for discrete-time linear time-invariant systems. *IEEE Transactions on Automatic Control*, **64**, 4717–4724.

Wang, Y., Olaru, S., Valmorbida, G., Puig, V., & Cembrano, G. (2019a). Set-invariance characterizations of discrete-time descriptor systems with application to active mode detection. *Automatica*, **107**, 255–263.

Wang, Y., Puig, V., & Cembrano, G. (2018a). Set-membership approach and Kalman observer based on zonotopes for discrete-time descriptor systems. *Automatica*, **93**, 435–443.

Wang, Y., Wang, Z., Puig, V., & Cembrano, G. (2019b). Zonotopic set-membership state estimation for discrete-time descriptor LPV systems. *IEEE Transactions on Automatic Control*, **64**, 2092–2099.

Wang, Z., Lim, C. C., & Shen, Y. (2018b). Interval observer design for uncertain discrete-time linear systems. *Systems & Control Letters*, **116**, 41–46.

Wang, Z., Shen, Y., Zhang, X., & Wang, Q. (2012). Observer design for discrete-time descriptor systems: An LMI approach. *Systems & Control Letters*, **61**, 683–687.

Xu, F., Tan, J., Wang, Y., Wang, X., Liang, B., & Yuan, B. (2019). Robust fault detection and set-theoretic uio for discrete-time LPV systems with state and output equations scheduled by inexact scheduling variables. *IEEE Transactions on Automatic Control*, **64**, 4982–4997.

Zhang, J., Swain, A. K., & Nguang, S. K. (2014). Robust $H_\infty$ adaptive descriptor observer design for fault estimation of uncertain nonlinear systems. *Journal of the Franklin Institute*, **351**, 5162–5181.

Zhang, Y., Yang, F., Han, Q., & Zhu, Y. (2019). A novel set-membership control strategy for discrete-time linear time-varying systems. *IET Control Theory & Applications*, **13**, 3087–3095.