Set-based state estimation and fault diagnosis of linear discrete-time descriptor systems using constrained zonotopes

Many physical processes such as battery packs, robotic systems with holonomic and nonholonomic constraints, and socioeconomic systems (Janschek, 2011; Yang et al., 2019), exhibit static relations between their internal variables. These processes are known as descriptor systems (or implicit systems), which have generalized dynamic and static behaviors described through differential and algebraic equations, respectively (Wang et al., 2018a). Descriptor systems appear in several contexts, such as linear control (Bender and Laub, 1987), fault-tolerant control (Shi and Patton, 2014), and fault diagnosis (Wang et al., 2019). However, few strategies can deal effectively with state estimation and fault diagnosis of descriptor systems when uncertainties with unknown probability distribution are present (Hamdi et al., 2012).

Set-based strategies for state estimation of discrete-time descriptor systems with unknown-but-bounded uncertainties is a recent subject, often addressed using intervals, zonotopes and ellipsoids (Efimov et al., 2015; Wang et al., 2018a, b; Merhy et al., 2019). Interval methods are used in Efimov et al. (2015) to design a state estimator for time-delay descriptor systems based on linear matrix inequalities and the Luenberger structure. However, interval arithmetics can lead to conservative enclosures due to the wrapping effect. A few methods are proposed in Wang et al. (2018a) for set-valued state estimation of linear descriptor systems by enclosing the intersection of two consistent sets with a zonotope bound. Nevertheless, since the intersection cannot be computed exactly, the resulting bound can be conservative. These strategies are extended in Wang et al. (2018b) for linear parameter-varying descriptor systems, but conservative enclosures are still present since the intersection method is maintained. Ellipsoids are used in Merhy et al. (2019) for state estimation of linear descriptor systems based on Luenberger type observers. Despite being able to provide stable bounds, the ellipsoidal estimation can be conservative since the complexity of the set is fixed and static relations are not directly incorporated. In addition, as in Wang et al. (2018a), restrictive assumptions on the rank of the system matrices are required to be able to design the proposed estimator. Moreover, due to the static constraints, the reachable sets of models of a descriptor system may be asymmetric even if the initial set is symmetric, and therefore methods based on the sets above can provide conservative enclosures.

Fault diagnosis aims to determine exactly which fault a process is subject to. For the case of descriptor systems, this problem has been considered in recent works using zonotopes (Yang et al., 2019; Wang et al., 2019). While Wang et al. (2019) explores the use of unknown input observers for robust passive fault diagnosis limited to additive faults, Yang et al. (2019) proposes a zonotope-based method for active fault diagnosis (AFD) of descriptor systems. The latter is based on the design of an input sequence for separation of the reachable sets. Unfortunately, the generator representation of zonotopes cannot incorporate exactly static relations between the state variables in general linear descriptor systems. Consequently, in practice, this may lead to a more difficult diagnosis.

This manuscript proposes strategies for set-based state estimation and AFD of linear descriptor systems based on constrained zonotopes (CZs), aiming to reduce the conservativeness described above through this asymmetric set representation. In contrast to the aforementioned sets, it is worth noting that linear static constraints, typical of descriptor systems, can be directly incorporated in the mathematical description of CZs. Thanks to this feature, set-valued methods based on CZs can provide less conservative enclosures. The proposed methods extensively explore the basic properties of CZs (Scott et al., 2016). In addition, efficient complexity reduction methods for CZs, available in the literature (Scott et al., 2016; Yang and Scott, 2018), are used in the proposed strategies allowing a direct trade-off between accuracy and computational efficiency. The superiority of the proposed approaches with respect to zonotope-based methods is highlighted in numerical examples.

This section presents a new method for set-based state estimation of system (1). For the sake of clarity, the prediction-update algorithm for standard discrete-time systems is first recalled. Consider (1) with ${\mathbf{E}}={\mathbf{I}}$. Given an initial condition $X_{0}$ and an input ${\mathbf{u}}_{k}$ with $k\geq 0$, let

$$\hat{X}_{0}=\{\mathbf{x}\in X_{0}:{\mathbf{C}}{\mathbf{x}}+{\mathbf{D}}{\mathbf{u}}_{0}+{\mathbf{D}}_{v}{\mathbf{v}}={\mathbf{y}}_{0},\,\mathbf{v}\in V\}.$$

(9)

Then, for $k\geq 1$ the prediction-update algorithm consists in computing sets $\bar{X}_{k}$ and $\hat{X}_{k}$ satisfying (Scott et al., 2016)

	$\displaystyle\bar{X}_{k}$	$\displaystyle=\{{\mathbf{A}}{\mathbf{x}}+{\mathbf{B}}{\mathbf{u}}_{k-1}+{\mathbf{B}}_{w}{\mathbf{w}}:\mathbf{x}\in\hat{X}_{k-1},\,\mathbf{w}\in W\},$		(10)
	$\displaystyle\hat{X}_{k}$	$\displaystyle=\{\mathbf{x}\in\bar{X}_{k}:{\mathbf{C}}{\mathbf{x}}+{\mathbf{D}}{\mathbf{u}}_{k}+{\mathbf{D}}_{v}{\mathbf{v}}={\mathbf{y}}_{k},\,\mathbf{v}\in V\},$		(11)

in which (10) is referred to as the prediction step, and (11) as the update step. The prediction step (10) must be reformulated when ${\mathbf{E}}$ is not invertible. We first consider the following assumption.

The set estimation in Scott et al. (2016) relied on (6)–(8) to compute the steps (10)–(11). In order to extend this method to the case of descriptor systems, when computing $\bar{X}_{k}$ it is necessary to take into account the static constraints arising from the possible singularity of ${\mathbf{E}}$. The proposed method is based on singular value decomposition (SVD). Let ${\mathbf{E}}={\mathbf{U}}{\mathbf{\Sigma}}{\mathbf{V}}^{T}$, where ${\mathbf{U}}$ and ${\mathbf{V}}$ are invertible by construction. Since ${\mathbf{E}}$ is square, then ${\mathbf{\Sigma}}$ is also square. Without loss of generality, let ${\mathbf{\Sigma}}$ be arranged as ${\mathbf{\Sigma}}=\text{blkdiag}(\tilde{{\mathbf{\Sigma}}},{\mathbf{0}})$, where $\tilde{{\mathbf{\Sigma}}}\in\mathbb{R}^{n_{z}\times n_{z}}$ is diagonal with all the $n_{z}=\text{rank}({\mathbf{E}})$ nonzero singular values of ${\mathbf{E}}$. Moreover, let ${\mathbf{z}}_{k}=(\tilde{{\mathbf{z}}}_{k},\check{{\mathbf{z}}}_{k})={{\mathbf{T}}}{\,\!}^{-1}{\mathbf{x}}_{k},\;\tilde{{\mathbf{z}}}_{k}\in\mathbb{R}^{n_{z}},\check{{\mathbf{z}}}_{k}\in\mathbb{R}^{n-n_{z}}$,

	$\displaystyle\begin{bmatrix}\tilde{{\mathbf{A}}}\\ \check{{\mathbf{A}}}\end{bmatrix}$	$\displaystyle=\begin{bmatrix}\tilde{{\mathbf{\Sigma}}}^{-1}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{I}}\end{bmatrix}{{\mathbf{U}}}{\,\!}^{-1}{\mathbf{A}}{\mathbf{T}},$		(12)
	$\displaystyle\begin{bmatrix}\tilde{{\mathbf{B}}}\\ \check{{\mathbf{B}}}\end{bmatrix}$	$\displaystyle=\begin{bmatrix}\tilde{{\mathbf{\Sigma}}}^{-1}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{I}}\end{bmatrix}{{\mathbf{U}}}{\,\!}^{-1}{\mathbf{B}},\begin{bmatrix}\tilde{{\mathbf{B}}}_{w}\\ \check{{\mathbf{B}}}_{w}\end{bmatrix}=\begin{bmatrix}\tilde{{\mathbf{\Sigma}}}^{-1}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{I}}\end{bmatrix}{{\mathbf{U}}}{\,\!}^{-1}{\mathbf{B}}_{w},$

with ${\mathbf{T}}={({\mathbf{V}}^{T})}{\,\!}^{-1}$, $\tilde{{\mathbf{A}}}\in\mathbb{R}^{n_{z}\times n}$, $\tilde{{\mathbf{B}}}\in\mathbb{R}^{n_{z}\times n_{u}}$, $\tilde{{\mathbf{B}}}_{w}\in\mathbb{R}^{n_{z}\times n_{w}}$. Then, system (1) can be rewritten with decoupled dynamics given by (Jonchkeere, 1988)

	$\displaystyle\tilde{{\mathbf{z}}}_{k}$	$\displaystyle=\tilde{{\mathbf{A}}}{\mathbf{z}}_{k-1}+\tilde{{\mathbf{B}}}{\mathbf{u}}_{k-1}+\tilde{{\mathbf{B}}}_{w}{\mathbf{w}}_{k-1},$		(13a)
	$\displaystyle{\mathbf{0}}$	$\displaystyle=\check{{\mathbf{A}}}{\mathbf{z}}_{k-1}+\check{{\mathbf{B}}}{\mathbf{u}}_{k-1}+\check{{\mathbf{B}}}_{w}{\mathbf{w}}_{k-1},$		(13b)
	$\displaystyle{\mathbf{y}}_{k}$	$\displaystyle={\mathbf{C}}{\mathbf{T}}{\mathbf{z}}_{k}+{\mathbf{D}}{\mathbf{u}}_{k}+{\mathbf{D}}_{v}{\mathbf{v}}_{k}.$		(13c)

Consider $W=\{{\mathbf{G}}_{w},{\mathbf{c}}_{w},{\mathbf{A}}_{w},{\mathbf{b}}_{w}\}$, and $\hat{X}_{0}$ given by (9). From (13b), the state ${\mathbf{z}}_{0}$ must satisfy $\check{{\mathbf{A}}}{\mathbf{z}}_{0}+\check{{\mathbf{B}}}{\mathbf{u}}_{0}+\check{{\mathbf{B}}}_{w}{\mathbf{w}}_{0}={\mathbf{0}}$. This constraint can be incorporated in the CG-rep of the initial set $\hat{Z}_{0}$ as follows. Let ${{\mathbf{T}}}{\,\!}^{-1}\hat{X}_{0}\triangleq\{{\mathbf{G}}_{0},{\mathbf{c}}_{0},{\mathbf{A}}_{0},{\mathbf{b}}_{0}\}$. Then, $\hat{Z}_{0}\triangleq\{\hat{{\mathbf{G}}}_{0},\hat{{\mathbf{c}}}_{0},\hat{{\mathbf{A}}}_{0},\hat{{\mathbf{b}}}_{0}\}$, with $\hat{{\mathbf{G}}}_{0}=[{\mathbf{G}}_{0}\,\;{\mathbf{0}}]$, $\hat{{\mathbf{c}}}_{0}={\mathbf{c}}_{0}$,

$\displaystyle\hat{{\mathbf{A}}}_{0}=\begin{bmatrix}{\mathbf{A}}_{0}&{\mathbf{0}}\\ \check{{\mathbf{A}}}{\mathbf{G}}_{0}&\check{{\mathbf{B}}}_{w}{\mathbf{G}}_{w}\end{bmatrix},\;\hat{{\mathbf{b}}}_{0}=\begin{bmatrix}{\mathbf{b}}_{0}\\ -\check{{\mathbf{A}}}{\mathbf{c}}_{0}-\check{{\mathbf{B}}}_{w}{\mathbf{c}}_{w}-\check{{\mathbf{B}}}{\mathbf{u}}_{0}\end{bmatrix}.$

Note that the extra columns in $\hat{{\mathbf{G}}}_{0}$ and $\hat{{\mathbf{A}}}_{0}$ come from ${\mathbf{w}}_{0}\in W$. Having defined $\hat{Z}_{0}$, for the purpose of state estimation the static relation (13b) can be shifted forward to time $k$ without loss of information. By doing so, from (13), the variables $\tilde{{\mathbf{z}}}_{k}$ are fully determined by (13a), while $\check{{\mathbf{z}}}_{k}$ are obtained a posteriori by $\check{{\mathbf{A}}}{\mathbf{z}}_{k}+\check{{\mathbf{B}}}{\mathbf{u}}_{k}+\check{{\mathbf{B}}}_{w}{\mathbf{w}}_{k}={\mathbf{0}}$.

Consider the set $Z_{\text{a}}={{\mathbf{T}}}{\,\!}^{-1}X_{\text{a}}=\{{{\mathbf{T}}}{\,\!}^{-1}{\mathbf{G}}_{\text{a}},{{\mathbf{T}}}{\,\!}^{-1}{\mathbf{c}}_{\text{a}},{\mathbf{A}}_{\text{a}},{\mathbf{b}}_{\text{a}}\}$, and let ${{\mathbf{T}}}{\,\!}^{-1}{\mathbf{c}}_{\text{a}}=[\tilde{{\mathbf{c}}}_{\text{a}}^{T}\,\;\check{{\mathbf{c}}}_{\text{a}}^{T}]^{T}$, ${{\mathbf{T}}}{\,\!}^{-1}{\mathbf{G}}_{\text{a}}=[\tilde{{\mathbf{G}}}_{\text{a}}^{T}\,\;\check{{\mathbf{G}}}_{\text{a}}^{T}]^{T}$. An effective enclosure of the prediction step for the descriptor system (13) can be obtained in CG-rep as follows.

{pf}

Since by assumption $({\mathbf{z}}_{k-1},{\mathbf{w}}_{k-1},{\mathbf{w}}_{k})\in\hat{Z}_{k-1}\times W\times W$, there exists $(\bm{\xi}_{k-1},\bm{\delta}_{k-1},\bm{\delta}_{k})\in B_{\infty}(\hat{{\mathbf{A}}}_{k-1},\hat{{\mathbf{b}}}_{k-1})\times B_{\infty}({\mathbf{A}}_{w},{\mathbf{b}}_{w})\times B_{\infty}({\mathbf{A}}_{w},{\mathbf{b}}_{w})$ such that ${\mathbf{z}}_{k-1}=\hat{{\mathbf{c}}}_{k-1}+\hat{{\mathbf{G}}}_{k-1}\bm{\xi}_{k-1}$, ${\mathbf{w}}_{k-1}={\mathbf{c}}_{w}+{\mathbf{G}}_{w}\bm{\delta}_{k-1}$, and ${\mathbf{w}}_{k}={\mathbf{c}}_{w}+{\mathbf{G}}_{w}\bm{\delta}_{k}$. Moreover, Assumption 1 implies that ${\mathbf{z}}_{k}\in Z_{\text{a}}$. Thus there must exist $\bm{\xi}_{\text{a}}\in B_{\infty}({\mathbf{A}}_{\text{a}},{\mathbf{b}}_{\text{a}})$ such that $\check{{\mathbf{z}}}_{k}=\check{{\mathbf{c}}}_{\text{a}}+\check{{\mathbf{G}}}_{\text{a}}\bm{\xi}_{\text{a}}$. Therefore, substituting these equalities in (13a) leads to

	$\displaystyle(\tilde{{\mathbf{z}}}_{k},\check{{\mathbf{z}}}_{k})=$	$\displaystyle(\tilde{{\mathbf{A}}}\hat{{\mathbf{c}}}_{k-1}{+}\tilde{{\mathbf{B}}}{\mathbf{u}}_{k-1}{+}\tilde{{\mathbf{B}}}_{w}{\mathbf{c}}_{w}$		(14)
		$\displaystyle{+}\tilde{{\mathbf{A}}}\hat{{\mathbf{G}}}_{k-1}\bm{\xi}_{k-1}{+}\tilde{{\mathbf{B}}}_{w}{\mathbf{G}}_{w}\bm{\delta}_{k-1},\check{{\mathbf{c}}}_{\text{a}}+\check{{\mathbf{G}}}_{\text{a}}\bm{\xi}_{\text{a}}).$

From the constraint (13b) shifted to time $k$, we have that

		$\displaystyle\check{{\mathbf{B}}}_{w}{\mathbf{c}}_{w}+\check{{\mathbf{B}}}_{w}{\mathbf{G}}_{w}\bm{\delta}_{k}+\check{{\mathbf{A}}}\begin{bmatrix}\tilde{{\mathbf{A}}}\hat{{\mathbf{c}}}_{k-1}+\tilde{{\mathbf{B}}}{\mathbf{u}}_{k-1}+\tilde{{\mathbf{B}}}_{w}{\mathbf{c}}_{w}\\ \check{{\mathbf{c}}}_{\text{a}}\end{bmatrix}$		(15)
		$\displaystyle+\check{{\mathbf{A}}}\begin{bmatrix}\tilde{{\mathbf{A}}}\hat{{\mathbf{G}}}_{k-1}&\tilde{{\mathbf{B}}}_{w}{\mathbf{G}}_{w}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{0}}&\check{{\mathbf{G}}}_{\text{a}}\end{bmatrix}\begin{bmatrix}\bm{\xi}_{k-1}\\ \bm{\delta}_{k-1}\\ \bm{\xi}_{\text{a}}\end{bmatrix}+\check{{\mathbf{B}}}{\mathbf{u}}_{k}={\mathbf{0}}.$

Rearranging (14) and (15), grouping $(\bm{\xi}_{k-1},\bm{\delta}_{k-1},\bm{\xi}_{\text{a}},$ $\bm{\delta}_{k})$, and writing in the CG-rep (2) proves the lemma. ∎

Lemma 1 provides a predicted enclosure of the state ${\mathbf{z}}_{k}$ in which the equality constraints (13b), shifted to time $k$, are directly taken into account. This is possible thanks to the fact that CZs incorporate equality constraints (see (2)). Finally, the prediction-update algorithm proposed for descriptor systems consists in the computation of CZs $\bar{Z}_{k}$, $\hat{Z}_{k}$, and $\hat{X}_{k}$, such that

	$\displaystyle\bar{Z}_{k}$	$\displaystyle=\{\bar{{\mathbf{G}}}_{k},\bar{{\mathbf{c}}}_{k},\bar{{\mathbf{A}}}_{k},\bar{{\mathbf{b}}}_{k}\},$		(16)
	$\displaystyle\hat{Z}_{k}$	$\displaystyle=\bar{Z}_{k}\cap_{{\mathbf{C}}{\mathbf{T}}}(({\mathbf{y}}_{k}-{\mathbf{D}}_{u}{\mathbf{u}}_{k})\oplus(-{\mathbf{D}}_{v}V_{k})),$		(17)
	$\displaystyle\hat{X}_{k}$	$\displaystyle={\mathbf{T}}\hat{Z}_{k}.$		(18)

For this algorithm, the initial set is $\hat{Z}_{0}$. The algorithm (16)–(18) operates recursively with $\bar{Z}_{k}$ and $\hat{Z}_{k}$ for $k\geq 1$ in the transformed state-space (13), while the estimated enclosure in the original state-space (1) is given by $\hat{X}_{k}$.

The previous section presented a method to address the problem of the set-based estimation of descriptor systems. In the following, this tool is used in the design of an AFD method accounting for a finite number of possible abrupt faults. Consider a linear discrete-time descriptor system whose dynamics obeys one of possible $n_{m}$ known models

	$\displaystyle{\mathbf{E}}^{[i]}{\mathbf{x}}_{k}^{[i]}$	$\displaystyle={\mathbf{A}}^{[i]}{\mathbf{x}}_{k-1}^{[i]}+{\mathbf{B}}^{[i]}{\mathbf{u}}_{k-1}+{\mathbf{B}}_{w}^{[i]}{\mathbf{w}}_{k-1},$		(19)
	$\displaystyle{\mathbf{y}}_{k}^{[i]}$	$\displaystyle={\mathbf{C}}^{[i]}{\mathbf{x}}_{k}^{[i]}+{\mathbf{D}}^{[i]}{\mathbf{u}}_{k}+{\mathbf{D}}_{v}^{[i]}{\mathbf{v}}_{k},$

for $k\geq 1$, with ${\mathbf{E}}^{[i]}\in\mathbb{R}^{n\times n}$, ${\mathbf{A}}^{[i]}\in\mathbb{R}^{n\times n}$, ${\mathbf{B}}^{[i]}\in\mathbb{R}^{n\times n_{u}}$, ${\mathbf{B}}_{w}^{[i]}\in\mathbb{R}^{n\times n_{w}}$, ${\mathbf{C}}^{[i]}\in\mathbb{R}^{n_{y}\times n}$, ${\mathbf{D}}^{[i]}\in\mathbb{R}^{n_{y}\times n_{u}}$, and ${\mathbf{D}}_{v}^{[i]}\in\mathbb{R}^{n_{y}\times n_{v}}$, $i\in\mathbb{I}=\{1,2,\ldots,n_{m}\}$. Moreover, $\text{rank}({\mathbf{E}}^{[i]})\leq n$, and let ${\mathbf{x}}_{0}^{[i]}\in X_{0}$, $({\mathbf{w}}_{k},{\mathbf{v}}_{k})\in W\times V$, and ${\mathbf{u}}_{k}\in U$, with $X_{0}$, $W$, $V$ and $U$ being known polytopic sets.

The goal of AFD is to find which model describes the process behaviour. In the following, the dynamics are assumed to not change during the diagnosis procedure, i.e. the AFD is fast enough to avoid the switching between models. In this sense, a sequence $({\mathbf{u}}_{0},{\mathbf{u}}_{1},...,{\mathbf{u}}_{N})$ of minimal length $N$ is designed such that any output ${\mathbf{y}}_{N}^{[i]}$ is consistent with only one $i\in\mathbb{I}$. If feasible, this problem may admit multiple solutions. For this reason, we introduce a cost function and select among the feasible input sequences the optimal one.

Let ${\overset{\to}{{\mathbf{u}}}}=({\mathbf{u}}_{0},...,{\mathbf{u}}_{N})\in\mathbb{R}^{(N+1)n_{u}}$, ${\overset{\to}{{\mathbf{w}}}}=({\mathbf{w}}_{0},\ldots$, ${\mathbf{w}}_{N})\in\mathbb{R}^{(N+1)n_{w}}$, and ${\overset{\to}{W}}=W\times\ldots\times W$. Consider a variable transformation similar to the one used in the previous section. With a slight abuse of notation, let ${\mathbf{z}}_{k}=(\tilde{{\mathbf{z}}}_{k},\check{{\mathbf{z}}}_{k})=({{\mathbf{T}}}{\,\!}^{-1}{\mathbf{x}}_{k},{\mathbf{w}}_{k}),$ $\tilde{{\mathbf{z}}}_{k}\in\mathbb{R}^{n_{z}}$, $\check{{\mathbf{z}}}_{k}\in\mathbb{R}^{n+n_{w}-n_{z}}$, with ${\mathbf{T}}^{[i]}={(({\mathbf{V}}^{[i]})^{T})}{\,\!}^{-1}$, ${\mathbf{V}}^{[i]}$ being obtained from the SVD ${\mathbf{E}}^{[i]}={\mathbf{U}}^{[i]}{\mathbf{\Sigma}}^{[i]}({\mathbf{V}}^{[i]})^{T}$. Then, (19) can be rewritten as

	$\displaystyle\tilde{{\mathbf{z}}}_{k}^{[i]}$	$\displaystyle=\tilde{{\mathbf{A}}}_{z}^{[i]}{\mathbf{z}}_{k-1}^{[i]}+\tilde{{\mathbf{B}}}^{[i]}{\mathbf{u}}_{k-1},$		(20a)
	$\displaystyle{\mathbf{0}}$	$\displaystyle=\check{{\mathbf{A}}}_{z}^{[i]}{\mathbf{z}}_{k}^{[i]}+\check{{\mathbf{B}}}^{[i]}{\mathbf{u}}_{k},$		(20b)
	$\displaystyle{\mathbf{y}}_{k}^{[i]}$	$\displaystyle={\mathbf{F}}^{[i]}{\mathbf{z}}_{k}^{[i]}+{\mathbf{D}}^{[i]}{\mathbf{u}}_{k}+{\mathbf{D}}_{v}^{[i]}{\mathbf{v}}_{k},$		(20c)

with ${\mathbf{F}}^{[i]}={\mathbf{C}}^{[i]}{\mathbf{T}}^{[i]}{\mathbf{P}}$, where ${\mathbf{P}}=[{\mathbf{I}}_{n}\,\;\bm{0}_{n\times n_{w}}]$, $\tilde{{\mathbf{A}}}_{z}^{[i]}=[\tilde{{\mathbf{A}}}^{[i]}\,\;\tilde{{\mathbf{B}}}_{w}^{[i]}]$, and $\check{{\mathbf{A}}}_{z}^{[i]}=[\check{{\mathbf{A}}}^{[i]}\,\;\check{{\mathbf{B}}}_{w}^{[i]}]$. Note that the $\tilde{(\cdot)}$ and $\check{(\cdot)}$ variables are defined according to (12) for each $i$, and equation (20b) has been already shifted to time $k$.

For each model $i$, consider the CZ $Z_{\text{a}}^{[i]}={({\mathbf{T}}^{[i]})}{\,\!}^{-1}X_{\text{a}}\times W=\{{\mathbf{G}}_{\text{a}}^{[i]},{\mathbf{c}}_{\text{a}}^{[i]},{\mathbf{A}}_{\text{a}}^{[i]},{\mathbf{b}}_{\text{a}}^{[i]}\}$, where $X_{\text{a}}$ satisfies Assumption 1, the set $\{{\mathbf{G}}_{z}^{[i]},{\mathbf{c}}_{z}^{[i]},{\mathbf{A}}_{z}^{[i]},{\mathbf{b}}_{z}^{[i]}\}={({\mathbf{T}}^{[i]})}{\,\!}^{-1}X_{0}\times W$, and define the initial feasible set $Z_{0}^{[i]}({\mathbf{u}}_{0})=\{{\mathbf{z}}\in{({\mathbf{T}}^{[i]})}{\,\!}^{-1}X_{0}\times W:\eqref{eq:systemSVDfaultconstraints}\text{ holds for }k=0\}$. This set is given by $Z_{0}^{[i]}({\mathbf{u}}_{0})=\{{\mathbf{G}}_{0}^{[i]},{\mathbf{c}}_{0}^{[i]},{\mathbf{A}}_{0}^{[i]},$ ${\mathbf{b}}_{0}^{[i]}({\mathbf{u}}_{0})\}$, where ${\mathbf{G}}_{0}^{[i]}={\mathbf{G}}_{z}^{[i]}$, ${\mathbf{c}}_{0}^{[i]}={\mathbf{c}}_{z}^{[i]}$,

$${\mathbf{A}}_{0}^{[i]}=\begin{bmatrix}{\mathbf{A}}_{z}^{[i]}\\ \check{{\mathbf{A}}}_{z}^{[i]}{\mathbf{G}}_{0}^{[i]}\end{bmatrix},\;{\mathbf{b}}_{0}^{[i]}({\mathbf{u}}_{0})=\begin{bmatrix}{\mathbf{b}}_{z}^{[i]}\\ -\check{{\mathbf{A}}}_{z}^{[i]}{\mathbf{c}}_{0}^{[i]}-\check{{\mathbf{B}}}^{[i]}{\mathbf{u}}_{0}\end{bmatrix}.$$

(21)

In addition, define the solution mappings $(\bm{\phi}_{k}^{[i]},\bm{\psi}_{k}^{[i]}):\mathbb{R}^{(k+1)n_{u}}\times\mathbb{R}^{n}\times\mathbb{R}^{(k+1)n_{w}}\times\mathbb{R}^{n_{v}}\to\mathbb{R}^{n+n_{w}}\times\mathbb{R}^{n_{y}}$ such that $\bm{\phi}_{k}^{[i]}({\overset{\to}{{\mathbf{u}}}},{\mathbf{z}}_{0},{\overset{\to}{{\mathbf{w}}}})$ and $\bm{\psi}_{k}^{[i]}({\overset{\to}{{\mathbf{u}}}},{\mathbf{z}}_{0},{\overset{\to}{{\mathbf{w}}}},{\mathbf{v}}_{k})$ are the state and output of (20) at $k$, respectively. Then, for each $i\in\mathbb{I}$, define state and output reachable sets at time $k$ as

	$\displaystyle Z_{k}^{[i]}({\overset{\to}{{\mathbf{u}}}}){=}$	$\displaystyle\{\bm{\phi}_{k}^{[i]}({\overset{\to}{{\mathbf{u}}}},{\mathbf{z}}_{0},{\overset{\to}{{\mathbf{w}}}}):({\mathbf{z}}_{0}^{[i]},{\overset{\to}{{\mathbf{w}}}})\in Z_{0}^{[i]}({\mathbf{u}}_{0})\times{\overset{\to}{W}}\},$
	$\displaystyle Y_{k}^{[i]}({\overset{\to}{{\mathbf{u}}}}){=}$	$\displaystyle\{\bm{\psi}_{k}^{[i]}({\overset{\to}{{\mathbf{u}}}},{\mathbf{z}}_{0},{\overset{\to}{{\mathbf{w}}}},{\mathbf{v}}_{k})\!:\!({\mathbf{z}}_{0},{\overset{\to}{{\mathbf{w}}}},{\mathbf{v}}_{k})\!\in\!Z_{0}^{[i]}({\mathbf{u}}_{0}){\times}{\overset{\to}{W}}{\times}V\}.$

Using (6)–(7), and taking note that by assumption ${\mathbf{z}}_{k}^{[i]}\in Z_{\text{a}}^{[i]}\subset\mathbb{R}^{n}\times W$ for every $k\geq 0$, the set $Z_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})$ is given by the CZ $\{{\mathbf{G}}_{N}^{[i]},{\mathbf{c}}_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}}),{\mathbf{A}}_{N}^{[i]},{\mathbf{b}}_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})\},$ where ${\mathbf{G}}_{N}^{[i]}$, ${\mathbf{c}}_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})$, ${\mathbf{A}}_{N}^{[i]}$, and ${\mathbf{b}}_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})$ are obtained by the recursive relations

	$\displaystyle{\mathbf{c}}_{k}^{[i]}({\overset{\to}{{\mathbf{u}}}}){=}\!\begin{bmatrix}\tilde{{\mathbf{A}}}_{z}^{[i]}{\mathbf{c}}_{k-1}^{[i]}({\overset{\to}{{\mathbf{u}}}})+\tilde{{\mathbf{B}}}^{[i]}{\mathbf{u}}_{k-1}\\ \check{{\mathbf{c}}}_{\text{a}}^{[i]}\end{bmatrix}\!,{\mathbf{G}}_{k}^{[i]}{=}\!\begin{bmatrix}\tilde{{\mathbf{A}}}_{z}^{[i]}{\mathbf{G}}_{k-1}^{[i]}\!\!&{\mathbf{0}}\\ {\mathbf{0}}&\check{{\mathbf{G}}}_{\text{a}}^{[i]}\end{bmatrix}\!,$
	$\displaystyle{\mathbf{A}}_{k}^{[i]}=\begin{bmatrix}{\mathbf{A}}_{k-1}^{[i]}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{A}}_{\text{a}}^{[i]}\\ \lx@intercol\hfil\check{{\mathbf{A}}}_{z}^{[i]}{\mathbf{G}}_{k}^{[i]}\hfil\lx@intercol\end{bmatrix},{\mathbf{b}}_{k}^{[i]}({\overset{\to}{{\mathbf{u}}}})=\begin{bmatrix}{\mathbf{b}}_{k-1}^{[i]}({\overset{\to}{{\mathbf{u}}}})\\ {\mathbf{b}}_{\text{a}}^{[i]}\\ -\check{{\mathbf{A}}}_{z}^{[i]}{\mathbf{c}}_{k}^{[i]}({\overset{\to}{{\mathbf{u}}}})-\check{{\mathbf{B}}}^{[i]}{\mathbf{u}}_{k}\end{bmatrix},$

for $k=1,2,\ldots,N$. Note that the third constraint in ${\mathbf{A}}_{k}^{[i]}$, ${\mathbf{b}}_{k}^{[i]}({\overset{\to}{{\mathbf{u}}}})$, comes from the fact that (20b) must hold. Using the initial values (21), the variables ${\mathbf{c}}_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})$ and ${\mathbf{b}}_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})$ can be written as explicit functions of the input sequence ${\overset{\to}{{\mathbf{u}}}}$ as

	$\displaystyle{\mathbf{c}}_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})$	$\displaystyle{=}\begin{bmatrix}\tilde{{\mathbf{A}}}_{z}^{[i]}\\ {\mathbf{0}}\end{bmatrix}^{N}\!\!\!\!{\mathbf{c}}_{z}^{[i]}{+}\!\!\sum_{m=1}^{N}\!\!\bigg{(}\begin{bmatrix}\tilde{{\mathbf{A}}}_{z}^{[i]}\\ {\mathbf{0}}\end{bmatrix}^{m\!-\!1}\!\!\begin{bmatrix}{\mathbf{0}}\\ \check{{\mathbf{c}}}_{\text{a}}^{[i]}\end{bmatrix}\bigg{)}{+}{\mathbf{H}}_{N}^{[i]}{\overset{\to}{{\mathbf{u}}}},$		(22)
	$\displaystyle{\mathbf{b}}_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})$	$\displaystyle=\bm{\alpha}_{N}^{[i]}+\bm{\Lambda}_{N}^{[i]}{\mathbf{c}}_{z}^{[i]}+\bm{\Omega}_{N}^{[i]}{\overset{\to}{{\mathbf{u}}}},$		(23)

where $\bm{\alpha}_{N}^{[i]}=\bm{\beta}_{N}^{[i]}+\bm{\Upsilon}_{N}^{[i]}{\mathbf{p}}_{N}^{[i]}$, $\bm{\Lambda}_{N}^{[i]}=\bm{\Upsilon}_{N}^{[i]}{\mathbf{Q}}_{N}^{[i]}$, $\bm{\Omega}_{N}^{[i]}=\bm{\Gamma}_{N}^{[i]}+\bm{\Upsilon}_{N}^{[i]}{\overset{\to}{{\mathbf{H}}}}{\,\!}^{[i]}$, $\bm{\beta}_{N}^{[i]}=\begin{bmatrix}[({\mathbf{b}}_{z}^{[i]})^{T}\;{\mathbf{0}}]&[({\mathbf{b}}_{\text{a}}^{[i]})^{T}\;{\mathbf{0}}]&\cdots&[({\mathbf{b}}_{\text{a}}^{[i]})^{T}\;{\mathbf{0}}]\end{bmatrix}^{T},$

	$\displaystyle{\mathbf{p}}_{N}^{[i]}=\bigg{[}\begin{matrix}\begin{bmatrix}{\mathbf{0}}\\ {\mathbf{0}}\end{bmatrix}^{T}&\begin{bmatrix}{\mathbf{0}}\\ \check{{\mathbf{c}}}_{\text{a}}^{[i]}\end{bmatrix}^{T}&\cdots&\sum_{m=1}^{N}\bigg{(}\begin{bmatrix}\tilde{{\mathbf{A}}}_{z}^{[i]}\\ {\mathbf{0}}\end{bmatrix}^{m-1}\begin{bmatrix}{\mathbf{0}}\\ \check{{\mathbf{c}}}_{\text{a}}^{[i]}\end{bmatrix}\bigg{)}^{T}\end{matrix}\bigg{]}^{T},$
	$\displaystyle\bm{\Upsilon}_{N}^{[i]}=\text{blkdiag}\big{(}\big{[}{\mathbf{0}}\,\;(-\check{{\mathbf{A}}}_{z}^{[i]})^{T}\big{]}^{T},\ldots,\big{[}{\mathbf{0}}\,\;(-\check{{\mathbf{A}}}_{z}^{[i]})^{T}\big{]}^{T}\big{)},$
	$\displaystyle\bm{\Gamma}_{N}^{[i]}=\text{blkdiag}\big{(}\big{[}{\mathbf{0}}\,\;(-\check{{\mathbf{B}}}^{[i]})^{T}\big{]}^{T},\ldots,\big{[}{\mathbf{0}}\,\;(-\check{{\mathbf{B}}}^{[i]})^{T}\big{]}^{T}\big{)},$
	$\displaystyle{\mathbf{Q}}_{N}^{[i]}=\bigg{[}\begin{matrix}\cdots&\bigg{(}\begin{bmatrix}\tilde{{\mathbf{A}}}_{z}^{[i]}\\ {\mathbf{0}}\end{bmatrix}^{\ell}\bigg{)}^{T}&\cdots\end{matrix}\bigg{]}^{T}\!\!,\;{\overset{\to}{{\mathbf{H}}}}{\,\!}^{[i]}=\big{[}\begin{matrix}\cdots&({\mathbf{H}}_{\ell}^{[i]})^{T}&\cdots\end{matrix}\big{]}^{T}\!\!,$
	$\displaystyle{\mathbf{H}}_{h}^{[i]}=\left[\begin{matrix}\cdots&\underbrace{\begin{bmatrix}\tilde{{\mathbf{A}}}_{z}^{[i]}\\ {\mathbf{0}}\end{bmatrix}^{h-m}\begin{bmatrix}\tilde{{\mathbf{B}}}^{[i]}\\ {\mathbf{0}}\end{bmatrix}}_{m=1,2,\ldots,h}&\cdots&\underbrace{\begin{bmatrix}{\mathbf{0}}\\ {\mathbf{0}}\end{bmatrix}}_{N-h+1\text{ terms}}&\cdots\end{matrix}\right],$

with $\ell=0,\ldots,N$. The variables $\bm{\beta}_{N}$, ${\mathbf{p}}_{N}$, $\bm{\Upsilon}_{N}$, and $\bm{\Gamma}_{N}$ have $N+1$ block matrices. In addition, the expression ${\mathbf{H}}_{h}^{[i]}$ holds for $h=1,\ldots,N$, while ${\mathbf{H}}_{0}^{[i]}=\bm{0}_{n\times(N{+}1)n_{u}}$.

Since $Z_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})$ is a CZ, the output reachable set $Y_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})$ is then a CZ obtained in accordance with (20c) as $Y_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})={\mathbf{F}}^{[i]}Z_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})\oplus{\mathbf{D}}^{[i]}{\mathbf{u}}_{N}\oplus{\mathbf{D}}_{v}^{[i]}V.$ Using properties (6) and (7), and letting $V=\{{\mathbf{G}}_{v},{\mathbf{c}}_{v},{\mathbf{A}}_{v},$ ${\mathbf{b}}_{v}\}$, this set is $Y_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})=\{{\mathbf{G}}_{N}^{Y[i]},{\mathbf{c}}_{N}^{Y[i]}({\overset{\to}{{\mathbf{u}}}})$, ${\mathbf{A}}_{N}^{Y[i]},$ ${\mathbf{b}}_{N}^{Y[i]}({\overset{\to}{{\mathbf{u}}}})\}$, with

	$\displaystyle{\mathbf{c}}_{N}^{Y[i]}({\overset{\to}{{\mathbf{u}}}})={\mathbf{F}}^{[i]}{\mathbf{c}}_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}})+{\mathbf{D}}^{[i]}{\mathbf{u}}_{N}+{\mathbf{D}}_{v}^{[i]}{\mathbf{c}}_{v},$		(24a)
	$\displaystyle{\mathbf{G}}_{N}^{Y[i]}=[{\mathbf{F}}^{[i]}{\mathbf{G}}_{N}^{[i]}\,\;{\mathbf{G}}_{v}],\;{\mathbf{A}}_{N}^{Y[i]}=\text{blkdiag}({\mathbf{A}}_{N}^{[i]},{\mathbf{A}}_{v}),$		(24b)
	$\displaystyle{\mathbf{b}}_{N}^{Y[i]}({\overset{\to}{{\mathbf{u}}}})=[({\mathbf{b}}_{N}^{[i]}({\overset{\to}{{\mathbf{u}}}}))^{T}\,\;{\mathbf{b}}_{v}^{T}]^{T}.$		(24c)