# Interval Estimation for Time-Varying Descriptor Systems via Simultaneous Optimizations of Multiple Interval Widths

Zhenhua Wang, Youdao Ma, Qinghua Zhang, Wentao Tang, Yi Shen

► To cite this version:

Zhenhua Wang, Youdao Ma, Qinghua Zhang, Wentao Tang, Yi Shen. Interval Estimation for Time-Varying Descriptor Systems via Simultaneous Optimizations of Multiple Interval Widths. IEEE Transactions on Systems, Man, and Cybernetics: Systems, 2024, 54 (6), pp.3774-3782. 10.1109/TSMC.2024.3371677. [https://hal-04842496v1](https://hal.science/hal-04842496v1)

HAL Id: hal-04842496

Submitted on 13 Jan 2025

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Interval Estimation for Time-Varying Descriptor Systems via Simultaneous Optimizations of Multiple Interval Widths

Zhenhua Wang, Member, IEEE, Youdao Ma, Qinghua Zhang, Wentao Tang, Yi Shen, Member, IEEE

*Abstract*—This paper investigates the state interval estimation problem for discrete-time linear time-varying descriptor systems subject to unknown but bounded system uncertainties. We propose a zonotope-based interval estimation method in an optimization framework. First, we present a novel zonotope-based interval estimator structure, in which the estimated interval bounds of each state component have design parameters independent of those of the other state components. Then the widths of the jointly estimated intervals enclosing every state component are simultaneously minimized by solving parallel $L_1$ optimization problems via linear programming. Finally, a simulation study shows the effectiveness and higher accuracy of the proposed method compared with existing methods.

*Index Terms*—Linear time-varying descriptor systems, interval estimation, optimization, zonotope

## I. INTRODUCTION

State estimation has been widely investigated in control theory and applications [1]. In practice, due to modeling and measurement uncertainties, single trajectory estimation would not converge to the true state. Traditionally, such uncertainties are handled with stochastic methods [2], [3]. These methods require that the probability distributions of the system uncertainties are available. However, this assumption may not be satisfied in practice. Alternatively, set membership estimation methods assume bounded uncertainties, without any probabilistic assumption [4]. Since this requirement is more practical, the set-membership estimation methods have attracted considerable attention and have been applied to various fields including controller design [5], cyber attack analysis [6] and fault diagnosis [7]–[9].

Set-membership estimation methods aim to compute a set sequence enclosing all the possible state trajectories [10]. Among the existing set-membership methods, interval estimation is widely used because it directly provides the range of each state component, and consequently, it is very easy to

check if any point of the state space is within the state uncertainty range [11]–[14]. In such results, the intervals bounding different state components are expressed independently (though they are often jointly estimated), the interval bounds constitute in fact a box (i.e., hyperrectangles) in the state space enclosing the state at every time instant. Generally speaking, interval estimation methods can be classified into two types. The first one is interval observer, which is developed based on the monotone theory [15]. In an interval observer, two sub-observers are constructed to estimate the upper and lower interval bounds, respectively. Due to its simple concept and high computational efficiency, interval observer has received much attention in the past two decades [16], [17]. However, the design of interval observer has a stringent requirement that the state estimation error dynamics are both stable and cooperative, which limits its scope of application [18]. Although many novel interval observer structures have been presented to relax this requirement, the performance is still below that of the method in [19]. In [19], the explicit solution of the estimation error is derived, based on which the reachable set of the state vector is tightly bounded by boxes. Unfortunately, this method is only suitable for linear time-invariant systems. While interval observes bound each state component independently, corresponding to the use of boxes enveloping the state trajectory, other set-membership estimation methods are based on more complex geometrical forms, for instance, polytopes [20], paralleltopes [21], and ellipsoids [22]–[24]. Though more accurate uncertainty sets are obtained with more complex geometrical sets, the much higher computation complexity and cost limit the use of such methods in practice. Remarkably, among such set-membership estimation methods, zonotope-based methods constitute a good balance between computational efficiency and estimation accuracy. Therefore, zonotope-based methods benefit from an increasing interest over the past decades [4], [25]–[29]. Nevertheless, these zonotope-based methods have a serious drawback: they adopt a zonotope representation for efficient propagation of uncertainties, but it is very costly to check if any given point of the state space belongs to such a zonotope [30]. For this reason, typically such zonotopes are converted to boxes when they are used in practice [4], [25]–[29], [31], [32]. In other words, zonotope-based set-membership estimation is used as an intermediate step of interval estimation.

It may appear superfluous to use zonotopes in this two-step interval estimation, because other methods directly estimate intervals, notably the aforementioned interval observers. It

This work was supported by National Natural Science Foundation of China under Grant No. 62373125 and by the Natural Science Foundation of Heilongjiang Province under Grant No. YQ2021F008.

Zhenhua Wang, Youdao Ma and Yi Shen are with the Department of Control Science and Engineering, Harbin Institute of Technology, Harbin, 150001 P. R. China (e-mail: zhenhua.wang@hit.edu.cn, shen@hit.edu.cn, youdao.ma@stu.hit.edu.cn).

Qinghua Zhang is with Inria, Université Gustave Eiffel, 35042 Rennes Cedex, France (e-mail: Qinghua.Zhang@inria.fr).

Wentao Tang is with Faculty of Electronic Information and Electrical Engineering, Dalian University Of Technology, Dalian 116024, PR China (e-mail: tang.wentao@dlut.edu.cn).

Manuscript received April 19, 2021; revised August 16, 2021.
is also possible to apply set-based methods directly using intervals (or boxes), but it would cumulate conservative approximation errors, resulting in quite over-estimated intervals. The zonotope-based two-step methods produce more accurate interval estimates, at the price of higher computational cost. Such zonotope-based methods can be designed by minimizing the error zonotope in some sense [4], [29], [31], [32]. Notably, the minimization of the Frobenius norm (F-norm) of the zonotope shape matrix leads to a numerically efficient algorithm similar to the Kalman filter [4], [28], [31]. However, the geometrical interpretation of this minimized criterion is not clear. On the contrary, the geometrical volume of zonotopes can be minimized [4], but it suffers from a high computation complexity. None of these existing design methods is based on the minimization of the obtained state interval bounds.

While most set membership-based state estimation methods are about state-space systems, few studies have considered descriptor systems, despite the fact that descriptor equations are a more convenient framework for modeling engineering systems from existing software libraries of subsystem models [33], [34]. Zonotope-based interval estimation for time-invariant descriptor systems has been studied in [35]. Although [28] and [29] have respectively developed interval estimation methods for linear time-varying descriptor systems and linear parameter-varying descriptor systems, their adopted optimization criteria, the F-norm of the generator matrices and the P-radius of zonotopes respectively, have no explicit geometrical interpretation.

This paper proposes a new zonotope-based interval estimation method for discrete-time linear time-varying descriptor systems, with the following contribution highlights: 1. To our knowledge, it is the first zonotope-based interval estimator designed by minimizing the estimated interval widths. 2. This interval estimator simultaneously and independently minimizes the widths of the intervals enclosing each state component. Though the intervals of different state components are jointly estimated, the parameters for different state components are independent from each other. 3. The involved $L_1$ optimization is based on linear programming, which is numerically efficient. 4. This new method considerably improves the interval estimation accuracy when evaluated with simulation examples against existing methods applicable to time-varying descriptor systems.

The remainder of this paper is structured as follows: In Section II, we give some preliminaries and the studied problem. In Section III, we provide the main results of the proposed zonotope-based interval estimation method which is designed in an optimization framework. In Section IV, we utilize a numerical simulation to demonstrate the effectiveness and superiority of our method. Finally, in Section V, we draw some conclusions.

## II. PRELIMINARIES AND PROBLEM STATEMENT

The notations used throughout this paper are standard. $\mathbb{R}$, $\mathbb{R}^n$, and $\mathbb{R}^{m \times n}$ denote the real numbers, the set of all $n$-dimensional vector, and the set of all $m \times n$-dimensional matrices, respectively. $I_n$ represents the $n$-dimensional identity

matrix. The matrix of all zeros and the vector of all ones are denoted by $0$ and $1$, respectively, where the dimension is unambiguously interpretable from the context. Consider a matrix $A \in \mathbb{R}^{m \times n}$, $A^T$ is its transpose, $A(i,j)$ is its $(i,j)$th element, $A(i,:)$ is its $i$th row, $rank(A)$ is its rank, and $A^\perp$ is a matrix satisfying $A^\perp A = 0$. Consider a vector $x \in \mathbb{R}^n$, $x^T$ and $x(i)$ are its transpose and $i$th element, respectively. The symbols $\geq$, $\leq$, and $| \cdot |$ on vectors or matrices should be understood element-wise. We use $\oplus$ to denote the Minkowski sum of two sets $\mathcal{X} \subseteq \mathbb{R}^n$ and $\mathcal{Y} \subseteq \mathbb{R}^n$, i.e.

$$ \mathcal{X} \oplus \mathcal{Y} = \{x + y : x \in \mathcal{X}, y \in \mathcal{Y}\}. \quad (1) $$

Herein, $\mathcal{X}$ and $\mathcal{Y}$ each can be a singleton set. For a set $\mathcal{X} \subseteq \mathbb{R}^n$, its linear image associated with a matrix $L \in \mathbb{R}^{m \times n}$ is denoted as

$$ L\mathcal{X} = \{Lx : x \in \mathcal{X}\}. \quad (2) $$

The following definitions, properties, and lemma are essential in this paper.

**Definition 1:** For a vector $x \in \mathbb{R}^n$, its 1-norm is defined as:

$$ \|x\|_1 = \sum_{i=1}^{n} |x(i)|. \quad (3) $$

**Definition 2:** A box (or interval vector) $\mathcal{I} \subset \mathbb{R}^n$ is defined, for a given pair $a \in \mathbb{R}^n$ and $b \in \mathbb{R}^n$, as:

$$ \mathcal{I} = [a, b] = \{x : x \in \mathbb{R}^n, a \leq x \leq b\}, \quad (4) $$

where the inequalities are component-wise, and the $i$th element of $\mathcal{I}$ is the interval $[a(i), b(i)]$.

**Definition 3:** [4] An $m$-order zonotope $\mathcal{Z} \subset \mathbb{R}^n$ is defined as:

$$ \mathcal{Z} = \langle p, H \rangle = \{p + Hz : z \in \mathcal{B}^m\}, \quad (5) $$

where $p \in \mathbb{R}^n$ is the center of $\mathcal{Z}$, $H \in \mathbb{R}^{n \times m}$ is the generator matrix of $\mathcal{Z}$, $\mathcal{B}^m \subset \mathbb{R}^m$ is the unit box with all its elements equal to $[-1, 1]$.

**Property 1:** [4], [31] For zonotopes, the following equations 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 (6) $$

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

$$ \langle p, H \rangle \subseteq \langle p, \mathcal{R}_q(H) \rangle, \quad (8) $$

where $p_1 \in \mathbb{R}^n$, $p_2 \in \mathbb{R}^n$ and $p \in \mathbb{R}^n$ are center vectors. $H_1 \in \mathbb{R}^{n \times m_1}$, $H_2 \in \mathbb{R}^{n \times m_2}$ and $H \in \mathbb{R}^{n \times m}$ are generator matrices. $L \in \mathbb{R}^{l \times n}$ is a linear image associated matrix. $\mathcal{R}_q(\cdot)$ is the order reduction operator as defined in [31] reducing the column number of $H$ from $m$ to $q$.

**Remark 1:** With the order reduction operator introduced in [31], $\mathcal{R}_q(H)$ in Property 1 is computed as

$$ \mathcal{R}_q(H) = \begin{cases} H, & q \geq m \\ [H_a \quad H_b], & n < q < m \end{cases}. \quad (9) $$

Herein, $H_a \in \mathbb{R}^{n \times (q-n)}$ is composed of the $q-n$ largest columns of $H$ in $L_2$ norm, and $H_b \in \mathbb{R}^{n \times n}$ is a diagonal
matrix with its diagonal entries

$$H_b(i, i) = \sum_{j=1}^{m-q+n} |\hat{H}(i, j)|, \quad i = 1, \dots, n, \quad (10)$$

where $\hat{H}$ is composed of the $m-q+n$ smallest columns of $H$ in $L_2$ norm.

**Property 2:** [4] Given an $m$-order zonotope $\mathcal{Z} = \langle p, H \rangle \subset \mathbb{R}^n$, the smallest box enclosing it is given by:

$$\mathbf{Box}(\mathcal{Z}) = [p - \tilde{H}, p + \tilde{H}], \quad (11)$$

where

$$\tilde{H} = \begin{bmatrix} \sum_{j=1}^{m} |H(1, j)| \\ \sum_{j=1}^{m} |H(2, j)| \\ \vdots \\ \sum_{j=1}^{m} |H(n, j)| \end{bmatrix}. \quad (12)$$

Based on Definition 1, $\tilde{H}$ can also be rewritten as:

$$\tilde{H} = \begin{bmatrix} \|H(1, : )\|_1 \\ \|H(2, : )\|_1 \\ \vdots \\ \|H(n, : )\|_1 \end{bmatrix}. \quad (13)$$

**Lemma 1:** [36] Given matrices $Q \in \mathbb{R}^{b \times c}$, and $R \in \mathbb{R}^{a \times c}$, if

$$\mathrm{rank}(Q) = \mathrm{rank} \left( \begin{bmatrix} Q \\ R \end{bmatrix} \right), \quad (14)$$

then the general solution of the underdetermined equation

$$PQ = R \quad (15)$$

for the unknown $P \in \mathbb{R}^{a \times b}$ is

$$P = RQ^{\dagger} + \Sigma(I_b - QQ^{\dagger}), \quad (16)$$

where $Q^{\dagger}$ is the Moore-Penrose pseudo inverse of $Q$, and $\Sigma$ is an arbitrary matrix of appropriate size.

Consider the discrete-time descriptor linear time-varying system formulated as follows:

$$Ex_{k+1} = A_k x_k + B_k u_k + D_k w_k, \quad (17a)$$

$$y_k = C x_k + F v_k, \quad (17b)$$

where (17a) and (17b) are the state equation and measurement equation of the considered system, respectively. $x_k \in \mathbb{R}^{n_x}$, $u_k \in \mathbb{R}^{n_u}$ and $y_k \in \mathbb{R}^{n_y}$ represent the system state, the control input and the measurement output, respectively. $w_k \in \mathbb{R}^{n_w}$ and $v_k \in \mathbb{R}^{n_v}$ are system uncertainties wherein $w_k$ denotes the process disturbances and $v_k$ denotes the measurement noise. $A_k \in \mathbb{R}^{n_x \times n_x}$, $B_k \in \mathbb{R}^{n_x \times n_u}$, $C \in \mathbb{R}^{n_y \times n_x}$, $D_k \in \mathbb{R}^{n_x \times n_w}$, $F \in \mathbb{R}^{n_y \times n_v}$ and $E \in \mathbb{R}^{n_x \times n_x}$ are known system matrices.

The aim of this paper is to compute an interval estimation of the system state $x_k$ from the descriptor model (17), the input-output signals and the uncertainty bounds. i.e. to find the upper bound $\bar{x}_k$ and the lower bound $\underline{x}_k$ such that

$$\underline{x}_k \le x_k \le \bar{x}_k, \quad \forall k \ge 0. \quad (18)$$

To this end, the following assumptions are made.

**Assumption 1:** The system matrices satisfy the following condition.

$$\mathrm{rank} \left( \begin{bmatrix} E \\ C \end{bmatrix} \right) = n_x. \quad (19)$$

**Assumption 2:** The initial system state $x_0$, and the system uncertainties $w_k$ and $v_k$ are unknown but bounded by zonotopes as follows:

$$x_0 \in \langle c_0, H_0 \rangle, \quad w_k \in \langle 0, I_{n_w} \rangle, \quad v_k \in \langle 0, I_{n_v} \rangle, \quad (20)$$

where $c_0 \in \mathbb{R}^{n_x}$ and $H_0 \in \mathbb{R}^{n_x \times n_x}$ are known.

### III. MAIN RESULTS

#### A. Zonotope-Based Interval Estimation

In order to design a new zonotope-based interval estimator, let us first make some simple transformations of the descriptor system model (17) as follows:

From the state equation (17a), we obtain

$$-A_{k-1}x_{k-1} + Ex_k = B_{k-1}u_{k-1} + D_{k-1}w_{k-1}, \quad (21)$$

and

$$x_{k-1} = x_{k-1}. \quad (22)$$

Then, according to the measurement equation (17b), we obtain

$$Cx_k = y_k - Fv_k, \quad (23)$$

and

$$Cx_{k-1} = y_{k-1} - Fv_{k-1}. \quad (24)$$

Define the augmented state vector as

$$z_k = \begin{bmatrix} x_{k-1} \\ x_k \end{bmatrix}. \quad (25)$$

For any $k \ge 1$, it is constrained as follows based on (21)-(24).

$$\Theta_k z_k = \begin{bmatrix} x_{k-1} \\ y_{k-1} - F v_{k-1} \\ B_{k-1} u_{k-1} + D_{k-1} w_{k-1} \\ y_k - F v_k \end{bmatrix}, \quad (26a)$$

$$x_k = M z_k, \quad (26b)$$

where

$$\Theta_k = \begin{bmatrix} I_{n_x} & 0 \\ C & 0 \\ -A_{k-1} & E \\ 0 & C \end{bmatrix}, \quad M = \begin{bmatrix} 0 & I_{n_x} \end{bmatrix}. \quad (27)$$

**Remark 2:** In this paper, the condition (19) is used to guarantee that $\Theta_k$ is column full rank. In some literature on descriptor systems, a weaker condition is used [37], as follows:

$$\mathrm{rank} \left( \begin{bmatrix} E \\ E^{\perp} A_k \\ C \end{bmatrix} \right) = n_x. \quad (28)$$

It is worth pointing out that this weaker condition is still available for our adopted state dynamics if we introduce the following equation to the dynamics.

$$E^{\perp} A_k x_k = -E^{\perp} B_k u_k - E^{\perp} D_k w_k, \quad (29)$$
where $E^\perp$ is a matrix satisfying $E^\perp E = 0$. For the sake of simplicity, we only consider the dynamics (26) and the condition (19) later in this paper.

To build a recursive zonotope-based estimator enclosing the state trajectory based on the constraints expressed in (26), at instant $k$, choose a matrix $S_k \in \mathbb{R}^{n_x \times (2n_x+2n_y)}$ such that

$$S_k \Theta_k = M, \quad (30)$$

and therefore

$$x_k = S_k \Theta_k z_k, \quad (31)$$

due to (26b). The existence and choice of $S_k$ satisfying (30) will be discussed later. Then straightforwardly it follows from (26a) that

$$x_k = S_k \begin{bmatrix} x_{k-1} \\ y_{k-1} - Fv_{k-1} \\ B_{k-1}u_{k-1} + D_{k-1}w_{k-1} \\ y_k - Fv_k \end{bmatrix}. \quad (32)$$

Assume that a zonotope enclosing $x_{k-1}$ was already computed at instant $k-1$. The uncertainty terms $v_{k-1}$, $w_{k-1}$ and $v_k$ are also bounded by zonotopes as stated in Assumption 2. Then the right-hand side of (32) corresponds to a zonotope, which in turn encloses the state vector $x_k$ according to (32). It is then clear that equation (32) contains essential elements for the design of new zonotope-based state estimators.

**Theorem 1:** Given an integer $q > n_x$ used in the order reduction operator $\mathcal{R}_q$ (see Property 1) and a matrix $S_k \in \mathbb{R}^{n_x \times (2n_x+2n_y)}$ satisfying (30) with $\Theta$ and $M$ as defined in (27), then the system state is bounded by a zonotope as:

$$x_k \in \langle c_k, H_k \rangle, \quad (33)$$

where $c_k$ and $H_k$ are recursively computed as

$$c_k = S_k \begin{bmatrix} c_{k-1} \\ y_{k-1} \\ B_{k-1}u_{k-1} \\ y_k \end{bmatrix}, \quad H_k = S_k Z_k, \quad (34)$$

in which

$$Z_k = \begin{bmatrix} \mathcal{R}_q(H_{k-1}) & 0 & 0 & 0 \\ 0 & -F & 0 & 0 \\ 0 & 0 & D_{k-1} & 0 \\ 0 & 0 & 0 & -F \end{bmatrix}. \quad (35)$$

Moreover, the interval bounds $\bar{x}_k$ and $\underline{x}_k$ of $x_k$ such that

$$\underline{x}_k \le x_k \le \bar{x}_k \quad (36)$$

are given by

$$\bar{x}_k = c_k + \begin{bmatrix} \|H_k(1, :)\|_1 \\ \|H_k(2, :)\|_1 \\ \vdots \\ \|H_k(n_x, :)\|_1 \end{bmatrix}, \quad \underline{x}_k = c_k - \begin{bmatrix} \|H_k(1, :)\|_1 \\ \|H_k(2, :)\|_1 \\ \vdots \\ \|H_k(n_x, :)\|_1 \end{bmatrix}. \quad (37)$$

*Proof:* From (20), we have $x_0 \in \langle c_0, H_0 \rangle$. Therefore, if we can prove that $x_{k-1} \in \langle c_{k-1}, H_{k-1} \rangle$ implies $x_k \in \langle c_k, H_k \rangle$, we will conclude that (33) holds for all $k \ge 0$.

Now we prove that $x_{k-1} \in \langle c_{k-1}, H_{k-1} \rangle$ implies $x_k \in$

$\langle c_k, H_k \rangle$. Based on Property 1, we obtain

$$\langle c_{k-1}, H_{k-1} \rangle \subseteq \langle c_{k-1}, \mathcal{R}_q(H_{k-1}) \rangle. \quad (38)$$

Therefore, $x_{k-1} \in \langle c_{k-1}, H_{k-1} \rangle$ leads to

$$x_{k-1} \in \langle c_{k-1}, \mathcal{R}_q(H_{k-1}) \rangle. \quad (39)$$

Meanwhile, using the unknown but bounded conditions in Assumption 2, the following inclusion holds

$$y_{k-1} - Fv_{k-1} \in y_{k-1} \oplus -F\langle 0, I_{n_v} \rangle. \quad (40)$$

It follows from Property 1 and (40) that

$$y_{k-1} - Fv_{k-1} \in \langle y_{k-1}, -F \rangle. \quad (41)$$

Similarly, we obtain

$$B_{k-1}u_{k-1} + D_{k-1}w_{k-1} \in \langle B_{k-1}u_{k-1}, D_{k-1} \rangle, \quad (42)$$

$$y_k - Fv_k \in \langle y_k, -F \rangle. \quad (43)$$

Notice that $S_k$ satisfies (30). Therefore, for system (17), following the derivations from (21) to (32), we can use (32) to establish a relationship between $x_k$ and $x_{k-1}$, and to iteratively compute the bounds of $x_k$. Substituting (39), (41)-(43) into (32) yields

$$x_k \in S_k \langle \zeta_k, Z_k \rangle = \langle S_k \zeta_k, S_k Z_k \rangle, \quad (44)$$

where $Z_k$ is computed by (35) and $\zeta_k$ is computed by

$$\zeta_k = \begin{bmatrix} c_{k-1} \\ y_{k-1} \\ B_{k-1}u_{k-1} \\ y_k \end{bmatrix}. \quad (45)$$

The result stated in (44) can be rewritten as

$$x_k \in \langle c_k, H_k \rangle, \quad (46)$$

with $c_k$ and $H_k$ as defined in (34). Therefore, the fact that $x_{k-1} \in \langle c_{k-1}, H_{k-1} \rangle$ implies $x_k \in \langle c_k, H_k \rangle$ is verified, and we can conclude (33) holds for all $k \ge 0$ from the above analysis.

Then, since $x_k$ belongs to the zonotope $\langle c_k, H_k \rangle$, we can compute the smallest box enclosing $\langle c_k, H_k \rangle$ based on Property 2 and use this box as the interval bounds. By substituting $p = c_k$ and $H = H_k$ into (11)-(13), the intervals enclosing the components of $x_k$ are computed as in (36) and (37). ■

**Remark 3:** The existing zonotope-based estimators presented in [28] and [29] are special forms of the proposed estimator. If we select $S_k$ as

$$S_k = [T_k A_{k-1} \quad L_k \quad T_k \quad N_k] \quad (47)$$

where $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}$ satisfying $T_k E + N_k C = I_{n_x}$ are gain matrices to be designed in [28], in this case $S_k \Theta_k = M$ exists, the proposed estimator will reduce to the estimator in [28]. Similarly, if we select $S_k$ as

$$S_k = \begin{bmatrix} (I_{n_x} - \Lambda_k C) T A_{k-1} & 0 & (I_{n_x} - \Lambda_k C) T \\ (I_{n_x} - \Lambda_k C) N + \Lambda_k & 0 \end{bmatrix} \quad (48)$$

where $T \in \mathbb{R}^{n_x \times n_x}$, $N \in \mathbb{R}^{n_x \times n_y}$ and $\Lambda_k \in \mathbb{R}^{n_x \times n_y}$
satisfying $TE + NC = I_{n_x}$ are gain matrices to be designed in [29], in this case $S_k\Theta_k = M$ exists, then the proposed estimator will reduce to the estimator in [29]. Compared with the estimators presented in [28] and [29], the proposed estimator has more degrees of design freedom, which can be used to improve the estimation accuracy by properly optimizing the parameters of the estimator. Therefore, the form of estimator (32) is more general than these existing estimators. More importantly, the form of (32) with design parameters in $S_k$ is particularly suitable for the simultaneous optimizations of multiple interval widths.

**Remark 4:** Notice that, from Assumption 1, the block-triangular matrix $\Theta_k$ in (27) is column full rank, hence it is straightforward to show that

$$ \mathrm{rank}(\Theta_k) = \mathrm{rank}\left(\begin{bmatrix} \Theta_k \\ M \end{bmatrix}\right). \quad (49) $$

Based on Lemma 1, there exists $S_k \in \mathbb{R}^{n_x \times (2n_x+2n_y)}$ satisfying the equation constraint (30).

Equation (30) counts $n_x \times 2n_x$ constraints on $S_k$, which contains $n_x \times (2n_x + 2n_y)$ unknowns. There are fewer equations than unknowns. The remaining degrees of freedom in $S_k$ can be explored to optimize some criterion as

$$ \min_{S_k} \mathcal{S}, \quad \text{s.t.} \quad S_k\Theta_k = M, \quad (50) $$

where $\mathcal{S}$ is a chosen criterion. Since the purpose of this paper is to obtain state interval bounds, it will dictate the criterion choice in the next subsection.

## B. An Optimization Method Based on 1-Norm

For the interval bounds given by Theorem 1 in (36), the interval widths are

$$ \xi_k = \bar{x}_k - \underline{x}_k. \quad (51) $$

Notice that $\xi_k \in \mathbb{R}^{n_x}$ is a vector of the same dimension as the state vector $x_k$, as well as its upper and low bounds $\bar{x}_k$ and $\underline{x}_k$. To optimize the estimated interval bounds, we choose the criterion as $\mathcal{S} = \xi_k$ and design the optimal $S_k$ by solving the following optimization problem.

$$ \min_{S_k} \xi_k, \quad \text{s.t.} \quad S_k\Theta_k = M. \quad (52) $$

**Remark 5:** Unlike in usual optimization problems, the criterion $\xi_k$ is an $n_x$ dimensional vector. It means that the optimization problem (52) minimizes simultaneously the $n_x$ components of $\xi_k$. These $n_x$ parallel minimizations will be independently solved based on the facts that, for every $i = 1, 2, \dots, n_x$, the $i$th component of $\xi_k$ depending on the $i$th row of $S_k$ which is independent of the other rows of $S_k$, and moreover, the constraints on the rows of $S_k$ in (52) are independent of each other. These facts will be confirmed in what follows.

The matrix constraint equation (30) can be separated into $n_x$ equations independently constraining each of the $n_x$ rows

of $S_k$:

$$ \begin{aligned} S_k(1, :)\Theta_k &= M(1, :), \\ S_k(2, :)\Theta_k &= M(2, :), \\ \vdots \\ S_k(n_x, :)\Theta_k &= M(n_x, :). \end{aligned} \quad (53) $$

Meanwhile, substituting (37) into (51) yields that

$$ \xi_k = 2 \left[ \|H_k(1, :)_1 \quad \|H_k(2, :)_1 \quad \cdots \quad \|H_k(n_x, :)_1\| \right]^T. \quad (54) $$

Then using $H_k(i, :) = S_k(i, :)Z_k$, (54) leads to

$$ \xi_k(i) = 2\|S_k(i, :)Z_k\|_1, \quad i = 1, \dots, n_x. \quad (55) $$

It is important to notice that, in (55) the $i$th component of $\xi_k$ depending on $i$th row of $S_k$ is independent of the other rows of $S_k$. Moreover, in (53) each row of $S_k$ is individually constrained, independently of the other rows. These facts make possible the simultaneous and independent minimization of every component of $\xi$. Then, from (52) and (53), the design of the optimal $S_k(i, :)$ can be formulated as the following $L_1$ optimization problem

$$ \min_{S_k(i,:)} \|S_k(i, :)Z_k\|_1, \quad \text{s.t.} \quad S_k(i, :)\Theta_k = M(i, :). \quad (56) $$

After solving these parallel $L_1$ optimization problems, we can assemble the optimal $S_k$ by concatenating the optimal row vectors $S_k(i, :)$, $i = 1, \dots, n_x$. i.e.

$$ S_k = [(S_k(1, :))^{\mathrm{T}} \quad (S_k(2, :))^{\mathrm{T}} \quad \cdots \quad (S_k(n_x, :))^{\mathrm{T}}]^{\mathrm{T}}. \quad (57) $$

**Remark 6:** Though the rows of $S_k$ are independently optimized, the intervals bounding $n_x$ state components are jointly estimated, due to the fact that the same zonotope is involved in all their computations, and the recursive computation of the zonotope centers $c_k$ made in (34) has interactions between different components of $c_k$.

The constrained $L_1$ optimization in (56) is similar to the basis pursuit problem [38], which considers the case where $Z_k$ is the identity matrix. The following theorem generalizes an important result of the basis pursuit problem to the problem formulated in (56) with a rectangular matrix $Z_k$. By applying this theorem with $\Theta = \Theta_k^T$, $\mu = (M(i, :))^T$ and $s = (S_k(i, :))^T$, it will be established that in the solution of the linear programming problem

$$ (a^*, s^*) = \arg \min_{a \in \mathbb{R}^{n_a}, s \in \mathbb{R}^{n_s}} \mathbf{1}^T a \quad (58a) $$

$$ \text{s.t.} \quad \begin{cases} \Theta_k^T s = M_i^T \\ -a \le Z_k^T s \le a \end{cases}, \quad (58b) $$

with $n_a$ the number of columns of $Z_k$ and $n_s = 2n_x+2n_y$, $s^*$ solves also the $L_1$ optimization problem formulated in (56), i.e.,

$$ s^* = \arg \min_{s \in \mathbb{R}^{n_s}} \|s^T Z_k\|_1 \quad \text{s.t.} \quad s^T \Theta_k = M_i. \quad (59) $$

**Theorem 2:** Given two matrices $\Theta \in \mathbb{R}^{2n_x \times n_s}$, $Z \in \mathbb{R}^{n_a \times n_s}$ and a vector $\mu \in \mathbb{R}^{2n_x}$, let $a^* \in \mathbb{R}^{n_a}$ and $s^* \in \mathbb{R}^{n_s}$ denote the solution of a linear programming problem defined with $\Theta$,
$Z$ and $\mu$, i.e.,

$$ (a^*, s^*) = \arg \min_{a \in \mathbb{R}^{n_a}, s \in \mathbb{R}^{n_s}} \mathbf{1}^T a \quad (60a) $$

$$ \text{s.t.} \quad \left\{ \begin{array}{l} \Theta s = \mu \\ -a \le Zs \le a \end{array} \right., \quad (60b) $$

then $s^*$ solves also the 1-norm minimization problem:

$$ s^* = \arg \min_{s \in \mathbb{R}^{n_s}} \|Zs\|_1 \quad \text{s.t.} \quad \Theta s = \mu. \quad (61) $$

*Proof:* The constrained optimization problem (60) involving both $a$ and $s$ will be decomposed in two steps in the sense that

$$ \min_{x,y} f(x,y) = \min_{y} (g(y)) \quad (62) $$

with a well defined function $g(y) = \min_x f(x, y)$.

The two side constraints in (60b) will be more compactly written as $|Zs| \le a$. Then the linear programming problem (60) becomes

$$ \min_{a \in \mathbb{R}^{n_a}, s \in \mathbb{R}^{n_s}} \mathbf{1}^T a \quad \text{s.t.} \quad \Theta s = \mu, \quad |Zs| \le a \quad (63) $$

$$ = \min_{s \in \mathbb{R}^{n_s}} \left( \min_{a \in \mathbb{R}^{n_a}} \mathbf{1}^T a \right) \quad \text{s.t.} \quad \Theta s = \mu, \quad |Zs| \le a \quad (64) $$

$$ = \min_{s \in \mathbb{R}^{n_s}, \Theta s = \mu} \left( \min_{a \in \mathbb{R}^{n_a}, |Zs| \le a} \mathbf{1}^T a \right). \quad (65) $$

Remark that, for the optimization problem inside the parentheses, $s$ is a given value, hence this optimization is only subject to the second constraint $|Zs| \le a$, which is now written under the min operator. Of course, for the other optimization problem, $s$ must satisfy the first constraint $\Theta s = \mu$, which is written under the corresponding min operator.

First consider the optimization inside the parentheses. For any given $s \in \mathbb{R}^n$, the value of $|Zs|$ is also given, with the matrix $Z$ specified in the statement of Theorem 1. Then the solution of

$$ \min_{a \in \mathbb{R}^{n_a}, |Zs| \le a} \mathbf{1}^T a \quad (66) $$

is trivially given by $a = |Zs|$, and $\mathbf{1}^T a = \mathbf{1}^T |Zs|$ is a well defined function of $s$. Hence

$$ \min_{a \in \mathbb{R}^{n_a}, |Zs| \le a} \mathbf{1}^T a = \mathbf{1}^T |Zs| = \|Zs\|_1. \quad (67) $$

Based on this result, (65) is then continued as

$$ \min_{s \in \mathbb{R}^{n_s}, \Theta s = \mu} \left( \min_{a \in \mathbb{R}^{n_a}, |Zs| \le a} \mathbf{1}^T a \right) = \min_{s \in \mathbb{R}^{n_s}, \Theta s = \mu} \|Zs\|_1. \quad (68) $$

It then follows from equalities (64), (65) and (68) that the solution of the linear programming problem (63) satisfies

$$ s^* = \arg \min_{s \in \mathbb{R}^{n_s}, \Theta s = \mu} \|Zs\|_1, \quad (69a) $$

$$ a^* = |Zs^*|. \quad (69b) $$

The constrained 1-norm minimization problem in (69a) is indeed the same as that in (61). Therefore, solving the linear programming problem (63), which is a shorthand of (60), yields $s^*$ minimizing also the constrained 1-norm (61). ■

To better illustrate the proposed method, an overall block diagram is depicted in Fig.1.

<figure><img src="image_2.png" /><figcaption>Fig. 1. The block diagram of the proposed method.</figcaption></figure>

*Remark 7:* The proposed method has the particularity of simultaneously minimizing individual interval widths for every state component, whereas other methods minimize some scaler criterion involving all the state components together. For this reason, we obtain the tightest interval for every state component, which will be illustrated in Section IV by the comparison study among the proposed method, the method in [28], and the method in [29].

*Remark 8:* The computational complexity of the proposed method mainly depends on solving the linear programming problem in (59). Usually, the linear programming problem can be solved by simplex method or interior point method. The worst computational complexity of the simplex method is exponential, whereas the interior point method has a polynomial computational complexity. Compared with the simplex method, the interior point method is more effective for solving linear programming problems. Therefore, we use the interior point method to solve the linear programming in (59) and then compute $S_k$ via (57). According to [39], the computational complexity of solving a linear programming problem with $n$ decision variables using interior point method is $O(n^{3.5}\ell)$, where $\ell$ is a parameter related to the computation precision. Since there exist $n_s + n_a$ decision variables in (59), the parallel computational complexity of rows of $S_k$ is $O((n_s + n_a)^{3.5}\ell)$, i.e. the computational complexity to compute $S_k$ is $O((n_s + n_a)^{3.5}\ell)$. Then, based on (34) and (35), the computational complexity of $c_k$ and $H_k$ are $n_x n_s + n_x n_u$ and $n_x n_s n_a + q \log(q)$, respectively. Finally, based on (37), the complexity of the simultaneous computation of $\bar{x}_k$ and $\underline{x}_k$ is $n_a + n_x$. Therefore, the computational complexity of the proposed method is $O((n_s + n_a)^{3.5}\ell)$, which is not very high.

## IV. SIMULATION

In this section, an example of a truck-trailer system from [28] is utilized to illustrate the performance of the proposed method. The simulated system is described by (17) with the following matrix parameters.

$$ E = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & -1 \\ 0 & 0 & 0 & 0 \end{bmatrix}, \quad C = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}, $$
$$
A_k = \begin{bmatrix} 1 - \theta_k T_s / l & 0 & 0 & 0 \\ \theta_k T_s / L & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & -\theta_k T_s / 2 & 0 & 1 \end{bmatrix}, B_k = \begin{bmatrix} \theta_k T_s / l \\ 0 \\ 0 \\ 0 \end{bmatrix}, \\
D_k = \operatorname{diag}(\begin{bmatrix} 0 & 0 & 0.006\theta_k & 0 \end{bmatrix}), \quad F = 0.002I_3,
$$

where $T_s = 0.2s$ is the sampling time, $\theta_k = -0.9 + 0.3\sin(0.0386k + 1)$, $L = 5.5m$ and $l = 4.8m$ are the length of trailer and truck, respectively. The control input is set as

$$u_k = K y_k + \frac{K_d L (y_k - y_{k-1})}{\theta_k T_s},$$

where

$$K = \begin{bmatrix} 3.29 & -0.65 & 0 \end{bmatrix}, \quad K_d = \begin{bmatrix} 3.11 & 0 & 0 \end{bmatrix}.$$

In this simulation, the initial system state is $x_0 = [-0.1745 \ 0.0873 \ 3.0 \ -0.0057]^T$. Hence it is bounded as $x_0 \in \langle c_0, H_0 \rangle$ wherein $c_0 = [-0.1 \ 0.08 \ 3 \ 0]^T$ and $H_0 = \operatorname{diag}([0.2 \ 0.05 \ 0.1 \ 0.05])$. The system uncertainties are random signals bounded as $w_k \in \langle 0, I_4 \rangle$ and $v_k \in \langle 0, I_3 \rangle$.

Based on these parameters, we adopt a comparison among the proposed method and the methods in [28] and [29] to show the advantages of the presented $L_1$ optimization method.

The method in [28] uses an observer to obtain a single trajectory $\hat{x}_k$ and analyses the zonotopic bounds $\langle 0, G_k \rangle$ of the estimation error. Then, the interval estimation results are obtained as follows:

$$\bar{x}_k = \hat{x}_k + \begin{bmatrix} \|G_k(1, :)\|_1 \\ \|G_k(2, :)\|_1 \\ \vdots \\ \|G_k(n_x, :)\|_1 \end{bmatrix}, \quad \underline{x}_k = \hat{x}_k - \begin{bmatrix} \|G_k(1, :)\|_1 \\ \|G_k(2, :)\|_1 \\ \vdots \\ \|G_k(n_x, :)\|_1 \end{bmatrix},$$

in which $\hat{x}_k$ and $G_k$ are obtained as

$$
\begin{align*}
\hat{x}_k &= T_k A_{k-1} \hat{x}_{k-1} + T_k B_{k-1} u_{k-1} + N_k y_k \\
&\quad + L_k (y_{k-1} - C \hat{x}_{k-1}), \\
G_k &= \begin{bmatrix} (T_k A_{k-1} - L_k C) \mathcal{R}_q (G_{k-1}) & T_k D_{k-1} \\ -L_k F & -N_k F \end{bmatrix},
\end{align*}
$$

where the optimal $T_k$, $N_k$ and $L_k$ satisfying $T_k E + N_k C = I_{n_x}$ are obtained using an $F$-norm-based optimization.

The method in [29] constructs a parameterized zonotope $\langle p_k, R_k \rangle$ to enclose the intersection zonotope of two zonotopes associated with the system model and measurement outputs. Then, the interval estimation results are obtained as follows:

$$\bar{x}_k = p_k + \begin{bmatrix} \|R_k(1, :)\|_1 \\ \|R_k(2, :)\|_1 \\ \vdots \\ \|R_k(n_x, :)\|_1 \end{bmatrix}, \quad \underline{x}_k = p_k - \begin{bmatrix} \|R_k(1, :)\|_1 \\ \|R_k(2, :)\|_1 \\ \vdots \\ \|R_k(n_x, :)\|_1 \end{bmatrix},$$

in which $p_k$ and $R_k$ are obtained as

$$
\begin{align*}
p_k &= (I_{n_x} - \Lambda_k C) T A_{k-1} p_{k-1} + (I_{n_x} - \Lambda_k C) T B_{k-1} u_{k-1} \\
    &\quad + (N - \Lambda_k C N + \Lambda_k) y_k, \\
R_k &= \begin{bmatrix} (I_{n_x} - \Lambda_k C) T A_{k-1} \mathcal{R}_q (R_{k-1}) & (I_{n_x} - \Lambda_k C) T D_{k-1} \\ -(I_{n_x} - \Lambda_k C) N F & -(I_{n_x} - \Lambda_k C) N F - \Lambda_k F \end{bmatrix},
\end{align*}
$$

where the optimal $T$, $N$ and $\Lambda_k$ satisfying $TE + NC = I_{n_x}$ are computed by solving an offline optimization problem such that the $P$-radius of zonotopes achieve the peak-to-peak performance.

For a fair comparison, we set $q = 20$ in the order reduction operator, which is the same value as that used in [28] and [29]. The interval estimation results provided by these three methods are depicted in Fig.2. To better analyze the performance of these methods, we consider the following mean interval width sum (MIWS) of each state component, which are relative measures of estimation accuracy

$$\text{MIWS}(i) = \frac{1}{N} \sum_{j=0}^{N} x_j(i), \quad i = 1, \dots, n_x,$$

where MIWS($i$) denotes MIWS of the $i$th state component and $N = 150$ is the total simulation steps.

It can be seen from Fig.2 that for all state components, the proposed method produces the tightest interval estimation results during the whole simulation process. Also from Table I, the estimated interval widths by the proposed method show a significant reduction than the other two methods. Compared with the method in [28], MIWS for each component decrease by 84.42%, 47.73%, 23.19% and 25.81%, respectively. Compared with the method in [29], MIWS for each component decrease by 96.71%, 83.51%, 78.01% and 77.45%, respectively. These simulation results confirm that the simultaneous minimization of the widths of the estimated interval bounds does improve the accuracy of interval estimation when compared with the methods in [28] and [29].

<table><caption>TABLE I
INDICES PROVIDED BY DIFFERENT METHODS.</caption><thead><tr><th>Index</th><th>The Proposed Method</th><th>Method 1</th><th>Method 2</th></tr></thead><tbody><tr><td>MIWS(1)</td><td>0.0214</td><td>0.1319</td><td>0.6506</td></tr><tr><td>MIWS(2)</td><td>0.0046</td><td>0.0088</td><td>0.0279</td></tr><tr><td>MIWS(3)</td><td>0.0053</td><td>0.0069</td><td>0.0241</td></tr><tr><td>MIWS(4)</td><td>0.0046</td><td>0.0062</td><td>0.0204</td></tr></tbody></table>

Method 1: The method in [28], Method 2: The method in [29].

V. CONCLUSION

In this paper, a zonotope-based method has been proposed to produce interval estimation results for discrete-time linear time-varying systems subject to unknown but bounded uncertainties. The proposed method is implemented in a framework of optimization. To produce tight interval estimation, a zonotope-based interval estimator is proposed. Then, in the design of the estimator, we present an $L_1$ optimization method based on 1-norm which has an explicit meaning of the acquired intervals. This $L_1$ optimization can be efficiently solved via linear programming. The performance of the proposed method is illustrated by a widely used example. Direct generalization of the proposed method to systems with parameter uncertainty would involve more conservatism than in the linear case. Extending the proposed method to time-varying descriptor
<figure><img src="image_3.png" /><figcaption>Fig. 2. Interval estimation results provided by the proposed method, the method in [28] and the method in [29].</figcaption></figure>

systems with parameter uncertainty will be done in future
works with important efforts for reducing conservatism.

REFERENCES

[1] Y. Hou, F. Zhu, X. Zhao, and S. Guo, "Observer design and unknown input reconstruction for a class of switched descriptor systems," *IEEE Transactions on Systems, Man, and Cybernetics: Systems*, vol. 48, no. 8, pp. 1411-1419, 2018.

[2] M. Zorzi, "Robust Kalman filtering under model perturbations," *IEEE Transactions on Automatic Control*, vol. 62, no. 6, pp. 2902-2907, 2017.

[3] X.-X. Ren and G.-H. Yang, "Modified Kalman filtering for stochastic nonlinear systems under non-Gaussian-Lévy noise and cyber attacks," *IEEE Transactions on Systems, Man, and Cybernetics: Systems*, vol. 53, no. 2, pp. 1222-1232, 2023.

[4] T. Alamo, J. M. Bravo, and E. F. Camacho, "Guaranteed state estimation by zonotopes," *Automatica*, vol. 41, no. 6, pp. 1035-1043, 2005.

[5] X. Ping, S. Yang, B. Ding, T. Raïssi, and Z. Li, "Observer-based output feedback robust MPC via zonotopic set-membership state estimation for LPV systems with bounded disturbances and noises," *Journal of the Franklin Institute*, vol. 357, no. 11, pp. 7368-7398, 2020.

[6] J. Li, Z. Wang, Y. Shen, and L. Xie, "Security synthesis for Cyber-Physical systems," *IEEE Transactions on Systems, Man, and Cybernetics: Systems*, vol. 53, no. 2, pp. 1027-1037, 2023.

[7] Z.-H. Zhang and G.-H. Yang, "Distributed fault detection and isolation for multiagent systems: An interval observer approach," *IEEE Transactions on Systems, Man, and Cybernetics: Systems*, vol. 50, no. 6, pp. 2220-2230, 2020.

[8] M. Zhou, Z. Cao, M. Zhou, J. Wang, and Z. Wang, "Zonotopic fault estimation for discrete-time LPV systems with bounded parametric uncertainty," *IEEE Transactions on Intelligent Transportation Systems*, vol. 21, no. 2, pp. 690-700, 2020.

[9] Z. Fei, X. Wang, and Z. Wang, "Event-based fault detection for unmanned surface vehicles subject to denial-of-service attacks," *IEEE Transactions on Systems, Man, and Cybernetics: Systems*, vol. 52, no. 5, pp. 3326-3336, 2022.

[10] F. Xu, "Generalized set-theoretic interval observer using element-wise nonnegativity transformation," *Automatica*, vol. 125, p. 109452, 2021.

[11] H. Ethabet, D. Rabehi, D. Efimov, and T. Raïssi, "Interval estimation for continuous-time switched linear systems," *Automatica*, vol. 90, pp. 230-238, 2018.

[12] W. Tang, Z. Wang, Y. Wang, T. Raïssi, and Y. Shen, "Interval estimation methods for discrete-time linear time-invariant systems," *IEEE Transactions on Automatic Control*, vol. 64, no. 11, pp. 4717-4724, 2019.

[13] Z. Wang, T. N. Dinh, Q. Zhang, T. Raïssi, and Y. Shen, "Fast interval estimation for discrete-time linear systems: An L1 optimization method," *Automatica*, vol. 137, p. 110029, 2022.

[14] S. Guo, W. Ren, C. K. Ahn, C. Wen, and H.-K. Lam, "Reachability analysis-based interval estimation for discrete-time Takagi-Sugeno fuzzy systems," *IEEE Transactions on Fuzzy Systems*, vol. 30, no. 6, pp. 1981-1992, 2022.

[15] D. Efimov, T. Raïssi, S. Chebotarev, and A. Zolghadri, "Interval state observer for nonlinear time varying systems," *Automatica*, vol. 49, no. 1, pp. 200-205, 2013.

[16] Y. Wang, D. M. Bevly, and R. Rajamani, "Interval observer design for LPV systems with parametric uncertainty," *Automatica*, vol. 60, pp. 79-85, 2015.

[17] J. Li, Z. Wang, Y. Shen, and Y. Wang, "Interval observer design for discrete-time uncertain Takagi-Sugeno fuzzy systems," *IEEE Transactions on Fuzzy Systems*, vol. 27, no. 4, pp. 816-823, 2019.

[18] T. Raïssi, D. Efimov, and A. Zolghadri, "Interval state estimation for a class of nonlinear systems," *IEEE Transactions on Automatic Control*, vol. 57, no. 1, pp. 260-265, 2012.

[19] W. Tang, Z. Wang, and Y. Shen, "Interval estimation for discrete-time linear systems: A two-step method," *Systems & Control Letters*, vol. 123, pp. 69-74, 2019.

[20] D. D. Huff, M. Fiacchini, and J. M. G. da Silva Jr, "Polyhedral estimates of the region of attraction of the origin of linear systems under aperiodic sampling and input saturation," *Automatica*, vol. 144, p. 110490, 2022.

[21] M. Maiga, N. Ramdani, L. Travé-Massuyès, and C. Combastel, "A comprehensive method for reachability analysis of uncertain nonlinear hybrid systems," *IEEE Transactions on Automatic Control*, vol. 61, no. 9, pp. 2341-2356, 2016.

[22] W. Zhang, Z. Wang, T. Raïssi, and Y. Shen, "Ellipsoid-based interval estimation for Lipschitz nonlinear systems," *IEEE Transactions on Automatic Control*, vol. 67, no. 12, pp. 6802-6809, 2022.

[23] S. Liu, Z. Wang, L. Wang, and G. Wei, "Recursive set-membership state estimation over a Flexray network," *IEEE Transactions on Systems, Man, and Cybernetics: Systems*, vol. 52, no. 6, pp. 3591-3601, 2022.

[24] Z. Wang, Y. Zhang, M. Shen, and Y. Shen, "Ellipsoidal set-membership filtering for discrete-time linear time-varying systems," *IEEE Transactions on Automatic Control*, vol. 68, no. 9, pp. 5767-5777, 2023.

[25] F. Zhu, Y. Tang, and Z. Wang, "Interval-observer-based fault detection and isolation design for TS fuzzy system based on zonotope analysis," *IEEE Transactions on Fuzzy Systems*, vol. 30, no. 4, pp. 945-955, 2022.
[26] J. Li, Z. Wang, C. K. Ahn, and Y. Shen, "Fault detection for Lipschitz nonlinear systems with restricted frequency-domain specifications," *IEEE Transactions on Systems, Man, and Cybernetics: Systems*, vol. 51, no. 12, pp. 7486-7496, 2021.

[27] M. Shen, T. Zhang, J. H. Park, Q.-G. Wang, and L.-W. Li, "Iterative proportional-integral interval estimation of linear discrete-time systems," *IEEE Transactions on Automatic Control*, vol. 68, no. 7, pp. 4249-4256, 2023.

[28] W. Tang, Z. Wang, Q. Zhang, and Y. Shen, "Set-membership estimation for linear time-varying descriptor systems," *Automatica*, vol. 115, p. 108867, 2020.

[29] Y. Wang, Z. Wang, V. Puig, and G. Cembrano, "Zonotopic set-membership state estimation for discrete-time descriptor LPV systems," *IEEE Transactions on Automatic Control*, vol. 64, no. 5, pp. 2092-2099, 2019.

[30] M. Althoff and J. J. Rath, "Comparison of guaranteed state estimators for linear time-invariant systems," *Automatica*, vol. 130, p. 109662, 2021.

[31] C. Combastel, "Zonotopes and Kalman observers: Gain optimality under distinct uncertainty paradigms and robust convergence," *Automatica*, vol. 55, pp. 265-273, 2015.

[32] V. T. H. Le, C. Stoica, T. Alamo, E. F. Camacho, and D. Dumur, "Zonotopic guaranteed state estimation for uncertain systems," *Automatica*, vol. 49, no. 11, pp. 3418-3424, 2013.

[33] Y. Wang, X. Xie, M. Chadli, S. Xie, and Y. Peng, "Sliding-mode control of fuzzy singularly perturbed descriptor systems," *IEEE Transactions on Fuzzy Systems*, vol. 29, no. 8, pp. 2349-2360, 2020.

[34] W. Tang, Q. Zhang, Z. Wang, and Y. Shen, "Set-membership estimation based on ellipsoid bundles for discrete-time LPV descriptor systems," *Automatica*, vol. 145, p. 110580, 2022.

[35] Y. Wang, V. Puig, and G. Cembrano, "Set-membership approach and Kalman observer based on zonotopes for discrete-time descriptor systems," *Automatica*, vol. 93, pp. 435-443, 2018.

[36] A. Ben-Israel and T. N. Greville, *Generalized inverses: theory and applications*. Springer Science & Business Media, 2003, vol. 15.

[37] M. Darouach, "$H_\infty$ unbiased filtering for linear descriptor systems via LMI," *IEEE Transactions on Automatic Control*, vol. 54, no. 8, pp. 1966-1972, 2009.

[38] D. L. Donoho and Y. Tsaig, "Fast solution of $l_1$-norm minimization problems when the solution may be sparse," *IEEE Transactions on Information Theory*, vol. 54, no. 11, pp. 4789-4812, 2008.

[39] M. Padberg, *Linear optimization and extensions*. Springer Science & Business Media, 2013, vol. 12.

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

**Qinghua Zhang** received the B.S. degree from the University of Science and Technology of China in 1986, the D.E.A. (Diplôme d'Etude Approfondie) from Université de Lille 1, France, in 1988, the PhD. and the H.D.R. (Habitation à Diriger des Recherches) both from Université de Rennes 1, France, in 1991 and 1999, respectively. During 1992, he was a post-doctorant at Linköping University, Sweden. Since 1993 he works at Institut National de Recherche en Informatique et en Automatique (Inria), France, as a researcher, and since 2001 as a senior researcher. His main research interests are in nonlinear system identification, fault diagnosis and signal processing.

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

**Wentao Tang** received his B.E. and Ph.D. degrees in control science and engineering from Harbin Institute of Technology, Harbin, China, in 2015 and 2021, respectively. He is an associate research fellow with the School of Control Science and Engineering, Dalian University of Technology. His research interests include fault diagnosis and fault-tolerant control.

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

**Yi Shen** received the B.E., M.Sc., and Ph.D. degrees in control science and engineering from the Harbin Institute of Technology, Harbin, China, in 1985, 1988, and 1995, respectively. He is a Professor at the School of Astronautics, Harbin Institute of Technology. His research interests include instrumentation and measurement, ultrasound signal processing, remote sensing, and modern detection techniques.

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

**Zhenhua Wang** received his B.Sc. degrees from Nanjing University of Aeronautics and Astronautics, Nanjing, China, in 2008 and his M.Sc. and Ph.D. degrees from Harbin Institute of Technology, Harbin, China, in 2010 and 2014, respectively. He is a professor at the Department of Control Science and Engineering at Harbin Institute of Technology. His research interests include fault diagnosis and fault-tolerant control. He is a member of the IFAC Technical Committee 1.1 on Modeling, Identification and Signal Processing.

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

**Youdao Ma** received his B.Sc. degrees from Harbin Institute of Technology, Harbin, China, in 2019. He is currently a Ph.D. candidate with the Department of Control Science and Engineering at Harbin Institute of Technology. His research interests include interval estimation and fault diagnosis.