Exact recursive updating of uncertainty sets for discrete-time plants with a lag

TL;DR

This paper introduces an exact recursive algorithm for updating the uncertainty sets of linear discrete-time systems with a lag, improving state estimation accuracy under bounded disturbances and noise.

Contribution

It presents a novel recursive method specifically for plants with a lag, enhancing the existing set-based state estimation techniques.

Findings

Efficient recursive updating of uncertainty sets for lag plants.

Improved accuracy in state estimation under bounded disturbances.

Applicable to systems with measurement noise and disturbances.

Abstract

In a recent paper [arXiv:1612.04918] there are new results concerning the polytopic set of possible states of a linear discrete-time SISO system subject to bounded disturbances from measurements corrupted by bounded noise. Using these results we construct an algorithm which, for the special case of a plant with a lag, recursively updates these polytopic sets when new measurements arrive.

Figures

Tables

Table 1. Table 1: The incidence matrix for S k + 1 subscript 𝑆 𝑘 1 S_{k+1} prior to the removal of rows and columns composed entirely of zeros

	$g^T\left(𝐕\right)$	$g^B\left(𝐕\right)$
$A^{\ast }{𝐟}_i$(with $u_k=1$)	${ℐ}^T$	$\mathrm{𝟎}$
$A^{\ast }{𝐟}_i$(with $u_k=-1$)	$\mathrm{𝟎}$	${ℐ}^B$
$A^{\ast }{𝐟}_i$(where ${\left(𝐟\right)}_1=0$)	$ℐO^T$	$ℐO^B$
$A^{\ast }{𝐟}_l^R$(Facets of $S_{k+1}$ from ridges)	$ℐR^T$	$ℐR^B$

Equations

$$\hat{\mathbf{n}}(\lambda )$$

$$\hat{\mathbf{d}}(\lambda )$$

$${\mathbf{x}}_{k+1}$$

$$y_k$$

$$z_k$$

$$\mathbf{A}$$

$$\mathbf{C}$$

$$P^{\ast }(\lambda )=-{\hat{\mathbf{n}}}_{\mathrm{dual}}(\lambda )/{\hat{\mathbf{d}}}_{\mathrm{dual}}(\lambda )$$

$${\mathbf{x}}_{k+1}^{\ast }$$

$$u_k^{\ast }$$

$${\mathbf{A}}^{\ast }$$

$${\mathbf{B}}^{\ast }$$

$${\mathbf{C}}^{\ast }$$

$$h_S(\mathbf{f})=\underset{\mathbf{x}\in S}{\max }⟨\mathbf{f},\mathbf{x}⟩,$$

$$H_S(\mathbf{f}):=\left\{\mathbf{x}\in {\mathbb{R}}^m:⟨\mathbf{f},\mathbf{x}⟩=h_S(\mathbf{f})\right\}$$

$$\left\{\mathbf{f}\in {\mathbb{R}}^m:⟨\mathbf{f},\mathbf{x}⟩=h_S(\mathbf{f})\right\}$$

$${\mathcal{N}}_S\left(\mathbf{x}\right)=\{{\lambda }_1{\mathbf{f}}_1+\dots +{\lambda }_n{\mathbf{f}}_n:{\lambda }_1,\dots ,{\lambda }_n\ge 0\},$$

$$S_{k+1}=\left\{\mathbf{x}^{+}\;:\;\begin{array}[c]{l}\mathbf{x}\in S_{k},\mathbf{x}^{+}=\mathbf{A}\mathbf{x}+\mathbf{B}u_{k},\\ y_{k}=\mathbf{C}\mathbf{x},\\ \left|u_{k}\right|\leq 1,\text{ }\left|y_{k}-z_{k}\right|\leq 1.\end{array}\right\}$$

$$\begin{array}[c]{c}u_{k}^{\ast}>0\implies u_{k}=1\\ u_{k}^{\ast}<0\implies u_{k}=-1\\ \left|u_{k}\right|<1\implies u_{k}^{\ast}=0\end{array}$$

$$\begin{array}[c]{c}y_{k}^{\ast}>0\implies y_{k}=1+z_{k}\\ y_{k}^{\ast}<0\implies y_{k}=-1+z_{k}\\ \left|y_{k}-z_{k}\right|<1\implies y_{k}^{\ast}=0.\end{array}$$

$$\begin{array}[c]{c}g^{\mathrm{T}}\left(\mathbf{x}\right)=\mathbf{A}\mathbf{x}+\mathbf{B}\\ g^{\mathrm{B}}\left(\mathbf{x}\right)=\mathbf{A}\mathbf{x}-\mathbf{B}.\\ \end{array}$$

$$F_i^{+}=\{{\mathbf{x}}_j^{+}:{\mathbf{x}}_j^{+}=g^{\mathrm{T}}\left({\mathbf{v}}_j\right)\text{ or }{\mathbf{x}}_j^{+}=g^{\mathrm{B}}\left({\mathbf{v}}_j\right),{\mathbf{v}}_j\in \mathrm{vert}(F_i)\}.$$

Full text

1 INTRODUCTION

Determining the possible states of a dynamic system from noisy measurements, when the input and measurement errors are unknown except for bounds on their magnitude, is a computationally challenging task. We use the term uncertainty set to describe such a set of states consistent with past measurements and a priori assumptions on the exact plant model and bounds on the input and measurement noise.

The problem of recursively calculating uncertainty sets was formulated, and in principle solved, in the seminal papers [21, 16, 2]. There is an extensive literature on the topic of uncertainty set propagation. See [6] and [12] for background, the survey paper [11] and the book [3]. Some of the many other papers which consider this problem are [4], [14], [19] and [20]. In real-time applications, for example fault detection and isolation [15], existing exact algorithms are of limited use because of their computational complexity. For this reason there has been a lot of research recently on the use of zonotopes and constrained zonotopes to approximate the exact polytopic uncertainty set, for example [5, 1, 17].

In [7] two theorems have been presented from which new algorithms can be derived which update polytopic uncertainty sets exactly for linear, SISO, discrete-time plants. We believe that these algorithms run faster than any current exact schemes. The particular form that the new algorithms take depends on properties of the plant, and in this paper we present details for the special case of plants without a feedthrough term, that is plants with a lag. The theorems quoted from [7] cover both the lag and non-lag cases, but the details of algorithm development for a plant with a lag are a lot simpler than the non-lag case treated in [7].

Exact, recursive methods currently require polytopic projection techniques, for example Fourier-Motzkin elimination, see [8], [13] and [18]. In these papers it is the identification of redundant inequality constraints that is most demanding computationally. The algorithm we present relies on updating the vertex-facet incidence matrix and does not generate any redundant constraints. Another advantage is that, along with the incidence matrix, both the vertices and facets of the uncertainty set are updated. Existing methods update either the facets or vertices, but not both. It is well known that converting from a representation of the polytope in terms of inequality constraints to a vertex representation is not easy, so having both at our disposal is valuable.

We give representative plots which demonstrate that use our algorithm requires a smaller update time than two other commonly used methods.

2 Basic Setup

The material in this section and the next is a summary of results in [7], specialized to the case of a plant with no direct feedthrough. The plant $P$, a linear, time-invariant, causal discrete-time, $m^{\rm th}$ order scalar system, is assumed known. There are two sources of uncertainty, an input noise disturbance $(u_{j})_{j=0}^{\infty}=\mathbf{u},$ and output measurement noise $(w_{j})_{j=0}^{\infty}=\mathbf{w}.$ The plant output is $(y_{j})_{j=0}^{\infty}=\mathbf{y}$, and the measurement at time $k$ is $z_{k}=y_{k}+w_{k}.$ The initial state, at time $k=0,$ is assumed to be known exactly, but nothing is known about the uncertainties except that they satisfy $\left|u_{j}\right|\leq 1$ and $\left|w_{j}\right|\leq 1.$

Given an initial state $\mathbf{x}_{0},$ the measurement history $z_{0},z_{1},\ldots,z_{k},$ and the plant dynamics, we seek the uncertainty set at time $k+1,$ denoted $S_{k+1};$ it is the set of possible states consistent with the measurements up to and including $z_{k}$, and is easily seen to be a closed, convex polytope.

2.1 Notation

The names of incidence vectors and matrices, whose elements are zeros or ones, will start with the letter $\mathcal{I}$. The $\lambda$-transform (generating function) of an arbitrary sequence $\mathbf{y}=(y_{k})_{k=0}^{\infty}$ is defined to be $\hat{\mathbf{y}}(\lambda):=\sum_{k=0}^{\infty}y_{k}\lambda^{k}.$ Real Euclidean space of dimension $m$ is denoted $\mathbb{R}^{m}$, where $m$ is the order of the plant $P$. States of the plant $P$ are represented by vectors, or points, in $\mathbb{R}^{m}$.

2.2 Transfer function description

The plant for the estimation system has the transfer function representation $P(\lambda)=\hat{\mathbf{n}}(\lambda)/\hat{\mathbf{d}}(\lambda)$ where

[TABLE]

$m\geq 1$, $\hat{\mathbf{n}}(\lambda)$ and $\hat{\mathbf{d}}(\lambda)$ are coprime polynomials with real coefficients, and it is assumed that both the plant $P(\lambda)$ and the plant $P^{\ast}(\lambda)$ for the dual system, defined below, are causal, implying $d_{0}\neq 0$ and $d_{m}\neq 0$. It is assumed that the plant has a lag so $n_{0}=0$, and without loss of generality we take $d_{0}=1.$

Assuming zero initial conditions, $\mathbf{y}$ and $\mathbf{u}$ are related by $\hat{\mathbf{d}}(\lambda)\hat{\mathbf{y}}(\lambda)=\hat{\mathbf{n}}(\lambda)\hat{\mathbf{u}}(\lambda).$ The indeterminate $\lambda$ represents the unit delay operator in the time domain.

2.3 State-space representations

A state-space description of the estimation system is

[TABLE]

where

[TABLE]

In (3f) $\mathbf{I}_{m-1}$ denotes the $m-1$ dimensional identity matrix, and $\mathbf{0}$ denotes a column vector of zeros of length $m-1$. The system matrix $\mathbf{A}$ is non-singular because $d_{m}\neq 0$.

There is a system closely related to the estimation system that we refer to as the dual system. Its input and output sequences are $(y_{j}^{\ast})_{j=1}^{\infty}$ and $(u_{j}^{\ast})_{j=1}^{\infty},$ and the dual plant, denoted $P^{\ast},$ has the transfer function representation

[TABLE]

where $\mathbf{n}_{\mathrm{dual}}=\left(n_{m+1},\ldots,n_{1}\right)$ and $\mathbf{d}_{\mathrm{dual}}=\left(d_{m+1},\ldots,d_{2},1\right)$. A minimal state-space realization of the dual system is

[TABLE]

[TABLE]

2.4 Polytopes

The primal and dual states defined in the previous Section will be interpreted in terms of the geometry of the polytopic uncertainty set $S_{k}$, so in this Section we introduce notation and briefly summarise the relevant theory of convex polytopes. For more information and background, including the definition of a polytope, see for example [10] or [22]. Let $S$ denote a closed, convex polytope.

The support function of $S$ in $\mathbb{R}^{m}$ is

[TABLE]

where $\mathbf{f}\in\mathbb{R}^{m}.$

When $\mathbf{f}\neq\mathbf{0}$ the set

[TABLE]

is the supporting hyperplane of $S$ with direction (outer normal vector) $\mathbf{f}.$ If $\mathbf{f}=\mathbf{0}$ then $H_{S}(\mathbf{f})=\mathbb{R}^{m}$.

The intersection of $S$ with a supporting hyperplane is called a face of $S,$ and a face of dimension $m-1$ is a facet of $S$. A face of dimension $m-2$ is called a ridge, and the faces of dimensions [math] and $1$ are termed vertices and edges, respectively.

The normal cone at a boundary point $\mathbf{x}$ is the set

[TABLE]

and is denoted $\mathcal{N}_{S}\left(\mathbf{x}\right)$. It is generated by the outward normals to the facets that form the polytope at $\mathbf{x}$, that is

[TABLE]

where $\mathbf{f}_{1},\ldots,\mathbf{f}_{n}$ are the directions of the facets containing $\mathbf{x}$. Thus $\mathcal{N_{S}}\left(\mathbf{x}\right)$ contains the directions of all hyperplanes which touch $S$ at $\mathbf{x.}$

The sets of vertices and boundary points of $S$ are respectively denoted $\operatorname*{vert}(S)$ and $\partial S$. The interior of $S$ is denoted $\operatorname*{int}\left(S\right)$. If $\mathbf{x}\in\operatorname*{int}\left(S\right)$, then $\mathcal{N_{S}}\left(\mathbf{x}\right):=\left\{\mathbf{0}\right\}$. By definition the directions of facets, and of the hyperplanes that contain facets, point outwards from the polytope. A direction is a non-zero vector.

The dual state $\mathbf{x}^{*}_{k}\left(\mathbf{y}^{*},\mathbf{u}^{*}\right)$ can be interpreted as a direction vector $\mathbf{f}$ in the normal cone of the primal state $\mathbf{x}_{k}\left(\mathbf{y},\mathbf{u}\right)\in S_{k}$. We shall sometimes drop the subscript $k$, and indicate the next time instant with the superscript +. For example, $\mathbf{f}$ will denote $\mathbf{x}^{*}_{k}$, and $\mathbf{f^{+}}$ will denote $\mathbf{x}^{*}_{k+1}$. Likewise, states $\mathbf{x}$ and facets $F$ belong to $S_{k}$, while $\mathbf{x}^{+}$ and $F^{+}$ belong to $S_{k+1}$.

2.4.1 Geometric and combinatorial description of $S_{k}$

Three matrices, $\mathbf{V}$, $\mathbf{F}$ and $\mathcal{I}$ are used to describe the polytope $S_{k}$, and the algorithm updates all of them. Geometric data are contained in $\mathbf{V}$ and $\mathbf{F}$, while combinatorial information on the structure of $S_{k}$ is contained in $\mathcal{I}$.

The $m$ by $n_{v}$ vertex matrix $\mathbf{V}$ has columns containing the coordinates of the vertices of $S_{k}$. 2. 2.

The $m$ by $n_{f}$ facet matrix $\mathbf{F}$ has columns containing the coordinates of the directions of the facets of $S_{k}$. Thus $\mathbf{F}=\left[\mathbf{f}_{1},\dotsc,\mathbf{f}_{n_{f}}\right],$ where $\mathbf{f}_{i}$ is a column vector representing the direction of facet number $i.$ The magnitude of $\mathbf{f}_{i}$ is not important, that is, for any $\alpha>0$, $\alpha\mathbf{f}_{i}$ can be used in place of $\mathbf{f}_{i}.$ 3. 3.

The vertex-facet incidence matrix of $S_{k}$ is the matrix $\mathcal{I}\in\{0,1\}^{n_{f}\times n_{v}}$ which has entry $\mathcal{I}(i,j)=1$ if $\mathbf{v}_{j}\in F_{i}$, and $\mathcal{I}(i,j)=0$ otherwise. So the facet $F_{i}$ is the convex hull of the vertices identified by the ones in the $i^{\rm th}$ row of $\mathcal{I}$.

3 Statement of Required Theorems

Following Witsenhausen, [21], $S_{k+1}$ is given recursively in terms of $S_{k}$ and the new observation $z_{k}$ by

[TABLE]

Special terminology is now introduced to describe states $\mathbf{x}$ and $\mathbf{x}^{+}$ related as in (23).

Definition 1**.**

The state $\mathbf{x}\in S_{k}$ is said to be a precursor of the state $\mathbf{x}^{+}$, $\mathbf{x}$ is propagated to $\mathbf{x}^{+},$ and $\mathbf{x}^{+}$ is a successor to $\mathbf{x},$ if there exists a scalar $u_{k}$ satisfying (1), (2), $\left|u_{k}\right|\leq 1$, and $\left|y_{k}-z_{k}\right|\leq 1$.

So $S_{k+1}$ is the set of all successors to all states in $S_{k}$, and any precursor of any state $\mathbf{x}^{+}\in S_{k+1}$ is in $S_{k}$.

We now define a relation between the inputs and outputs of the primal and dual systems.

Definition 2**.**

The scalar pair $\left(y_{k},u_{k}\right)$ is said to be aligned with $\left(y_{k}^{\ast},u_{k}^{\ast}\right)$ if

[TABLE]

and

[TABLE]

Associated with any $\mathbf{f}\in\mathbb{R}^{m}$ define in the $y^{\ast}_{k}u^{\ast}_{k}$-plane the dual line $L^{\ast}\left(\mathbf{f}\right)$:

Definition 3**.**

$L^{\ast}\left(\mathbf{f}\right)=\{(y^{*}_{k},u^{*}_{k}):n_{m}y^{\ast}_{k}+d_{m}u^{\ast}_{k}=-\left(\mathbf{f}\right)_{1}\}.$**

The scalars $y^{*}_{k}$ and $u^{*}_{k}$, the input and output of the dual plant at time $k$, are constrained to lie on the line $L^{\ast}\left(\mathbf{f}\right)$ by (6) and (22).

The square $Q$ in the $y_{k}u_{k}$-plane contains points $(y_{k},u_{k})$ satisfying $\left|u_{k}\right|\leq 1$ and $\left|y_{k}-z_{k}\right|\leq 1$.

One more definition is required in order to state the Theorems quoted from [7].

Definition 4**.**

Given $\mathbf{x}\in S_{k}$, $\mathbf{f}\in\mathcal{N}_{S_{k}}(\mathbf{x})$ and $z_{k},$ the set $M\left(\mathbf{x},\mathbf{f},z_{k}\right)$ is the set of scalar pairs $\left(u_{k},y_{k}^{\ast}\right)$ which satisfy

$(y_{k},u_{k})\in Q$, where $y_{k}=\mathbf{C}\mathbf{x}$ and 2. 2.

$\left(y_{k},u_{k}\right)$* is aligned with $\left(y_{k}^{\ast},u_{k}^{\ast}\right)$, where $(y_{k}^{\ast},u_{k}^{\ast})\in L^{\ast}\left(\mathbf{f}\right).$*

Finding $M$ is computationally very simple. For example, in Fig. 1 alignment occurs solely between the two circled points, so $M$ contains the singleton pair $(u_{k},y^{*}_{k})=(1,0)$.

We now present the Theorems from which the algorithm is derived. See Fig. 2 for a geometric depiction of the vectors in Theorem 5 for the important special case where $\mathbf{x}$ and $\mathbf{x}^{+}$ are both boundary points. Proofs are in [7].

Theorem 5**.**

Suppose $\mathbf{x}\in S_{k}$ and $\mathbf{f}\in\mathcal{N}_{S_{k}}(\mathbf{x})$. Then $\mathbf{x^{+}}=\mathbf{A}\mathbf{x}+\mathbf{B}u_{k}\in S_{k+1}$ and $\mathbf{f^{+}}=\mathbf{A}^{\ast}\mathbf{f}+\mathbf{B}^{\ast}y_{k}^{\ast}\in\mathcal{N}_{S_{k+1}}(\mathbf{x^{+}})$ if and only if $\left(u_{k},y_{k}^{\ast}\right)\in M\left(\mathbf{x},\mathbf{f},z_{k}\right)$.

In a typical application a point on the boundary of $S_{k}$, and the direction of any one of its supporting hyperplanes, are propagated to a point on the boundary of $S_{k+1}$ along with the direction of one of its supporting hyperplanes.

The companion result Theorem 6, given below, shows that any boundary point $\mathbf{x^{+}}\in S_{k+1}$, and any direction in the normal cone of $\mathbf{x^{+}}$, are attainable from any precursor $\mathbf{x}$ of $\mathbf{x^{+}}$ and some direction in the normal cone of $\mathbf{x}$.

Theorem 6**.**

Select any $\mathbf{x^{+}}\in S_{k+1}$, any $\mathbf{f^{+}}\in\mathcal{N}_{S_{k+1}}\left(\mathbf{x^{+}}\right)$ and any precursor $\mathbf{x}$ of $\mathbf{x^{+}}$. There exists $\mathbf{f}\in\mathcal{N}_{S_{k}}\left(\mathbf{x}\right)$ and $\left(u_{k},y_{k}^{\ast}\right)\in M\left(\mathbf{x},\mathbf{f},z_{k}\right)$ for which $\mathbf{x^{+}}=\mathbf{A}\mathbf{x}+\mathbf{B}u_{k}$ and $\mathbf{f^{+}}=\mathbf{A}^{\ast}\mathbf{f}+\mathbf{B}^{\ast}y_{k}^{\ast}$.

In the rest of the paper we show how to use the propagation of states and their normal cones to update the uncertainty set for plants with a lag. An implementation of our algorithm in Matlab code, uncertaintysetwithlag.m, is available on the link below111Go to https://mathworks.com/matlabcentral/fileexchange. Also given is code tcomplag.m which allows comparisons of computation time for our algorithm, denoted F-V (facets-vertices), Fourier-Motzkin and mplp (multi-parametric linear programming). We used the software package [9] for Fourier-Motzkin projection and mplp. See [18] for an explanation of the use of Fourier-Motzkin projection to compute updated uncertainty sets.

4 The incidence matrix for $S_{k}$

We turn now to the construction of $\mathcal{I}^{+}$, the incidence matrix for $S_{k+1}$, shown in Table 1, where the six incidence matrices in the body of the Table are defined in Section 5. Initially $\mathcal{I}^{+}$ is put equal to the $\left(3n_{f}+n_{R}\right)$ by $2n_{v}$ zero matrix, that is $\mathcal{I}^{+}$ is allocated sufficient space to accommodate all conceivable facets and vertices of $S_{k+1}.$ So there are $3n_{f}+n_{R}$ rows in the Table, where $n_{f}$ and $n_{R}$ are respectively the number of facets and ridges in $S_{k}$. The rows are divided into fours block, each of the first three blocks comprising $n_{f}$ rows, and the last block $n_{R}$ rows. The first two blocks point to facets of $S_{k+1}$ whose directions have the property $\left(\mathbf{f}\right)_{m}\neq 0$, while the last two blocks point to facets whose directions $\mathbf{f}$ have the property $\left(\mathbf{f}\right)_{m}=0$.

Turning now to the columns of Table 1, the first block of $n_{v}$ columns points to states $\mathbf{A}\mathbf{v}_{j}+\mathbf{B}$, \mathbf{v}_{l}\in\operatorname{vert}\bigl{(}S_{k}\bigr{)}, which, depending on $z_{k}$, might be vertices of $S_{k+1}$. The second block of $n_{v}$ columns points to states $\mathbf{A}\mathbf{v}_{j}-\mathbf{B}$, also possible vertices of $S_{k+1}$.

After $\mathcal{I}^{+}$ has been constructed in accordance with Table 1 it must have some rows and columns composed entirely of zeros. The final step is for these rows and columns of zeros to be deleted.

5 Filling in Table 1

We introduce notation to describe two ways a state can be propagated from $S_{k}$ to $S_{k+1}$. The mappings $g^{\mathrm{T}}$ and $g^{\mathrm{B}}$, where the superscripts indicate the top and bottom sides of $Q$, map a state on the boundary of $S_{k}$ to the boundary of $S_{k+1}$.

On the top (bottom) side of $Q$, $u_{k}=1$ $(u_{k}=-1)$, and (1) is used to propagate $\mathbf{x}$.

Definition 7**.**

For any $\mathbf{x}\in S_{k}$ put

[TABLE]

Definition 8**.**

A point $\mathbf{x}\in S_{k}$ is said to propagate up (resp. down) if $g^{\mathrm{T}}\left(\mathbf{x}\right)$ \bigl{(}(g^{\mathrm{B}}\left(\mathbf{x}\right)\bigr{)} lies on the boundary of $S_{k+1}$.

These propagations occur when the maximum or minimum allowable value of the input ($\pm$1) is applied to the plant, and furthermore the resulting state of the plant lies on the boundary of $S_{k+1}$. Theorem 5 tells us, given $\mathbf{x}$ and $\mathcal{N}_{S_{k}}\left(\mathbf{x}\right)$, if $\mathbf{x}$ propagates up or down. For example, in Fig. 1, if $\mathbf{x}\in\partial S_{k}$ then $\mathbf{x}$ propagates up to $g^{\mathrm{T}}\left(\mathbf{x}\right)$ because there exists $\mathbf{0}\neq\mathbf{f}\in\mathcal{N}_{S_{k}}\left(\mathbf{x}\right)$ which implies $(u_{k},y^{*}_{k})=(1,0)\in M\left(\mathbf{x},\mathbf{f},z_{k}\right)$ and so, by Theorem 5, $\mathbf{x}^{+}=\mathbf{A}\mathbf{x}+\mathbf{B}\in S_{k+1}$ and $\mathbf{f^{+}}=\mathbf{A}^{\ast}\mathbf{f}\in\mathcal{N}_{S_{k+1}}(\mathbf{x^{+}})$ and is non-zero, further implying that $\mathbf{x}^{+}\in\partial S_{k+1}$.

Definition 8 can be extended to facets as follows.

Definition 9**.**

Suppose $F$ is a facet $S_{k}$, and $F^{+}$ is a facet $S_{k+1}$. Then $F$ is said to propagate up (down) to $F^{+}$ if $F^{+}=g^{T}(F)$, \bigl{(}F^{+}=g^{B}(F)\bigr{)}.

Here $g^{\mathrm{T}}\left(F\right)$, for example, denotes the set $\left\{g^{\mathrm{T}}\left(\mathbf{x}\right):\mathbf{x}\in F\right\}$.

If $F$ propagates either up or down to $F^{+}$ then $F^{+}$ is an affinely isomorphic image of $F$ because $\mathbf{A}$ is invertible.

Let $F_{i}$ be any facet of $S_{k}$, and denote its direction by $\mathbf{f}_{i}$. There are three possibilities:

The line L^{*}\bigl{(}\mathbf{f}_{i}\bigr{)} intersects the positive $u^{*}_{k}$ axis, that is d_{m}\bigl{(}\mathbf{f}_{i}\bigr{)}_{1}<0. The $n_{f}\times 1$ logical vector $\mathcal{I}F^{\mathrm{T}}$ points to facets with such directions. Explicitly, $\mathcal{I}F^{\mathrm{T}}(i)=1$ if $\mathbf{f}_{i}$ satisfies d_{m}\bigl{(}\mathbf{f}_{i}\bigr{)}_{1}<0, and $\mathcal{I}F^{\mathrm{T}}(i)=0$ otherwise; 2. 2.

L^{*}\bigl{(}\mathbf{f}_{i}\bigr{)} intersects the negative $u^{*}_{k}$ axis, that is d_{m}\bigl{(}\mathbf{f}_{i}\bigr{)}_{1}>0. Define $\mathcal{I}F^{\mathrm{B}}(i)=1$ if $\mathbf{f}_{i}$ satisfies d_{m}\bigl{(}\mathbf{f}_{i}\bigr{)}_{1}>0, and $\mathcal{I}F^{\mathrm{B}}(i)=0$ otherwise; and 3. 3.

L^{*}\bigl{(}\mathbf{f}_{i}\bigr{)} passes through the origin. Define $\mathcal{I}F^{\mathrm{O}}(i)=1$ if $\mathbf{f}_{i}$ satisfies \bigl{(}\mathbf{f}_{i}\bigr{)}_{1}=0, and $\mathcal{I}F^{\mathrm{O}}(i)=0$ otherwise.

5.0.1 Case 1

For every $\mathbf{x}_{j}$ in $F_{i}$ we have the situation depicted in Fig. 1. By Theorem 5, $\mathbf{x}_{j}$ is propagated up, and $F_{i}$ is propagated up to the facet F_{i}^{+}:=g^{T}\bigl{(}F_{i}\bigr{)} of $S_{k+1}$, which has direction $\mathbf{A}^{*}\mathbf{f}_{i}$. The logical matrix $\mathcal{I}^{\mathrm{T}}\in\{0,1\}^{n_{f}\times n_{v}}$ is defined to have entry $\mathcal{I}^{\mathrm{T}}(i,j)=1$ if $\mathcal{I}(i,j)=1$ and $\mathcal{I}F^{\mathrm{T}}(i)=1$, and $\mathcal{I}^{\mathrm{T}}(i,j)=0$ otherwise. A one in the $i^{\rm th}$ row and $j^{\rm th}$ column of $\mathcal{I}^{\mathrm{T}}$ implies that there is a facet, $F_{i}^{+}$, of $S_{k+1}$ with direction $\mathbf{A}^{*}\mathbf{f}_{i}$, and that $g^{\mathrm{T}}\left(\mathbf{v}_{j}\right)\in F_{i}^{+}$.

Case 2) is similar to Case 1), except that $F_{i}$ is propagated down. The matrix $\mathcal{I}^{\mathrm{B}}\in\{0,1\}^{n_{f}\times n_{v}}$ is defined to have entry $\mathcal{I}^{\mathrm{B}}(i,j)=1$ if $\mathcal{I}(i,j)=1$ and $\mathcal{I}F^{\mathrm{B}}(i)=1$, and $\mathcal{I}^{\mathrm{B}}(i,j)=0$ otherwise.

We require for Case 3), covered below, as well as for the case of propagation of facets from ridges, to be treated in the next section, two logical matrices which point to those vertices of $S_{k}$ that are propagated to vertices of $S_{k+1}$. The $n_{v}\times 1$ vector $\mathcal{I}V^{\mathrm{PT}}$ is defined by $\mathcal{I}V^{\mathrm{PT}}(j)=1$ if the $j^{\rm th}$ column of $\mathcal{I}^{\mathrm{T}}$ is not composed entirely of zeros, and $\mathcal{I}V^{\mathrm{PT}}(j)=0$ otherwise. Thus $\mathcal{I}V^{\mathrm{PT}}$ points to those vertices which propagate up to a vertex of $S_{k+1}$. Likewise $\mathcal{I}V^{\mathrm{PB}}$ is defined by $\mathcal{I}V^{\mathrm{PB}}(j)=1$ if the $j^{\rm th}$ column of $\mathcal{I}^{\mathrm{B}}$ is not composed entirely of zeros, and $\mathcal{I}V^{\mathrm{PB}}(j)=0$ otherwise.

Every point of $Q$ is aligned with the origin in the $y^{*}_{k}u^{*}_{k}$ plane so, in Case 3), by Theorem 5 every $\mathbf{x}_{j}$ in $F_{i}$ is propagated to $\mathbf{x}_{j}^{+}=\mathbf{A}\mathbf{x}_{j}+\mathbf{B}u_{k}\in\partial S_{k+1}$ for all $\left|u_{k}\right|\leq 1$. We show the method for constructing all vertices of the facet $F_{i}^{+}$, and provide a sketch of the proof that it works. Put $F_{i}^{+}=\{\mathbf{x}_{j}^{+}:\mathbf{x}_{j}^{+}=\mathbf{A}\mathbf{x}_{j}+\mathbf{B}u_{k}\text{ where }\mathbf{x}_{j}\in F_{i},u_{k}\in[-1,1]\}.$ A convexity argument shows

[TABLE]

We seek the vertices of the facet $F_{i}^{+}$. Now every precursor $\mathbf{v}_{j}$ of any $\mathbf{v}_{j}^{+}$ in the vertex set of $F_{i}^{+}$ is a vertex of $F_{i}$. However, it is not necessarily the case that every successor of every vertex of $F_{i}$ is a vertex of $F_{i}^{+}$. We must remove from $F_{i}^{+}$ those $\mathbf{x}_{j}^{+}$ which are not vertices of $F_{i}^{+}$, and this can be done by removing those $\mathbf{x}_{j}^{+}$ which are not vertices of $S_{k}$.

Define the incidence matrix $\mathcal{I}O^{\mathrm{T}}(i,j)=1$ if $\mathcal{I}F^{\mathrm{O}}(i)=1$, $\mathcal{I}V^{\mathrm{PT}}(j)=1$ and $\mathcal{I}(i,j)=1$, and $\mathcal{I}O^{\mathrm{T}}(i,j)=0$ otherwise. In words, $\mathcal{I}O^{\mathrm{T}}(i,j)=1$ if $\mathbf{v}_{j}$ propagates up to a vertex of $S_{k}$, and $\mathbf{v}_{j}$ is in $F_{i}$, with $(\mathbf{f}_{i})_{1}=0$.

The matrix $\mathcal{I}O^{\mathrm{B}}$, pointing to those vertices of facets with direction vector whose first component is zero that propagate down, is defined similarly. Thus $\mathcal{I}O^{\mathrm{B}}(i,j)=1$ if $\mathcal{I}F^{\mathrm{O}}(i)=1$, $\mathcal{I}V^{\mathrm{PB}}(j)=1$ and $\mathcal{I}(i,j)=1$, and $\mathcal{I}O^{\mathrm{B}}(i,j)=0$ otherwise.

This completes the construction of all of Table 1 except the last block.

5.1 Propagation from a ridge in $S_{k}$ to a facet in $S_{k+1}$

In the previous Section Case 3) was more intricate than Cases 1) and 2) because it treated facets with directions whose dual lines passed through the origin, opening up an infinite number of alignment possibilities. There is one, and only one, other situation in which a boundary point $\mathbf{x}$ of $S_{k}$ can have a direction $\mathbf{f}$ in its normal cone which satisfies \bigl{(}\mathbf{f}\bigr{)}_{1}=0. If such $\mathbf{x}$ does not lie on a facet covered by case 3) above, $\mathbf{x}$ must belong to a special type of ridge, to be defined next.

See Fig. 3 where $S_{k}$ possesses two intersecting facets, with directions $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$, whose first components have opposite signs. These facets have the state $\mathbf{x}$ in their intersection, so necessarily there exists a direction in the normal cone of $\mathbf{x}$ whose first component is zero. We define $R$ to be the set of all such ridges of $S_{k}$, and denote its cardinality by $n_{R}$. Clearly $R$ is non-empty if $S_{k}$ is non-empty and full-dimensional.

Definition 10**.**

Define $R$ to be the set of ridges, $R_{l}$, $l=1,\ldots,n_{R}$, of $S_{k}$ with the property that $R_{l}$ is the intersection of two facets $F_{1}$ and $F_{2}$ say, for which $\left(\mathbf{f}_{1}\right)_{1}<0$ and $\left(\mathbf{f}_{2}\right)_{1}>0$, where $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ are the directions of $F_{1}$ and $F_{2}$.

Associated with any ridge in $R$ there is by definition a direction in the cone generated by $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ having the property that its first component is zero. This direction is unique.

Definition 11**.**

For any $R_{l}\in R$ define $\mathbf{f}_{l}^{R}$ to be the direction in the convex hull of $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ for which $\left(\mathbf{f}_{l}^{R}\right)_{1}=0$.

Clearly $\mathbf{f}_{l}^{R}\in\mathcal{N}_{S_{k}}\left(\mathbf{x}\right)$ for all $\mathbf{x}\in R_{l}$. The proof of the following Proposition is omitted.

Proposition 12**.**

For all $R_{l}\in R$ there is a facet of $S_{k}$ with direction $\mathbf{A}^{*}\mathbf{f}_{l}^{R}$.

For any $S_{k}$ and any $R_{l}\in R$, the direction $\mathbf{f}_{l}^{R}$ can be found using linear interpolation. More precisely $\mathbf{f}_{l}^{R}=t\mathbf{f}_{1}+(1-t)\mathbf{f}_{2}$ where $t=\left(\mathbf{f}_{1}\right)_{1}/[\left(\mathbf{f}_{2}\right)_{1}-\left(\mathbf{f}_{1}\right)_{1}]$, and $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ are the directions of the facets which intersect to form $R_{l}$.

The incidence matrices in the bottom row of Table 1 can now be defined.

The vertex set of $R_{l}$ is the intersection of the vertex sets of the facets which intersect to form $R_{l}$. The matrix $\mathcal{I}RV\in\{0,1\}^{{n_{R}}\times n_{v}}$ has entry $\mathcal{I}RV(l,j)=1$ if the ridge $R_{l}$ contains $\mathbf{v}_{j}$, and $\mathcal{I}RV(l,j)=0$ otherwise. So the $l^{\rm th}$ row of $\mathcal{I}RV$ points to the vertex set of $R_{l}$.

The matrix $\mathcal{I}R^{\mathrm{T}}\in\{0,1\}^{{n_{R}}\times n_{v}}$ has entry $\mathcal{I}R^{\mathrm{T}}(l,j)=1$ if $\mathcal{I}V^{\mathrm{PT}}(j)=1$ and $\mathcal{I}RV(l,j)=1$, and $\mathcal{I}R^{\mathrm{T}}(l,j)=0$ otherwise. So the $l^{\rm th}$ row of $\mathcal{I}R^{\mathrm{T}}$ points to vertices $\mathbf{v}_{j}$ of $R_{l}$ which propagate up to vertices of the ridge $R_{l}$ in $S_{k}$.

The matrix $\mathcal{I}R^{\mathrm{B}}\in\{0,1\}^{{n_{R}}\times n_{v}}$ has entry $\mathcal{I}R^{\mathrm{B}}(l,j)=1$ if $\mathcal{I}V^{\mathrm{PB}}(j)=1$ and $\mathcal{I}RV(l,j)=1$, and $\mathcal{I}R^{\mathrm{T}}(l,j)=0$ otherwise. So the $l^{\rm th}$ row of $\mathcal{I}R^{\mathrm{B}}$ points to vertices $\mathbf{v}_{j}$ of $R_{l}$ which propagate down to vertices of the ridge $R_{l}$ in $S_{k}$.

This completes the definition of all of the incidence matrices occurring in Table 1. We claim that $\mathcal{I}^{+}$ constructed according to Table 1 is, after removal of rows and columns containing all zeros, the incidence matrix of $S_{k+1}$. We have sketched the proof, which ultimately relies on Theorems 5 and 6.

6 Examples

Theorems 5 and 6 are most naturally interpreted as yielding $S_{k+1}$, consistent with measurements up to $z_{k}$. The input is $S_{k}$ consistent with measurements up to $z_{k-1}$. In common terminology we “update then propagate”. It is more common to determine $S_{k+1}$ consistent with measurements up to $z_{k+1}$, with input $S_{k}$ consistent with measurements up to $z_{k}$, essentially “propagate then update”. The accompanying code implements this latter case. For plants with a lag changing the order of the operations update and propagate is very easy.

We give representative plots comparing our algorithm, denoted F-V, projection using Fourier-Motzkin elimination, and projection using multi-parametric linear programming. They were generated using the Matlab code accompanying the paper. The reader can generate plots for her/himself with differing plant sizes and lengths of measurement sequences. Figure 4 compares the computation times of the three methods for a third order randomly selected stable plant, with randomly selected measurements. Figure 5 compares the computation times for a fifth order randomly selected stable plant, again with randomly chosen measurements. It can be seen that the algorithm presented in this paper is significantly faster. For these examples the F-V algorithm was between twenty and forty times faster, on average, than the better of the two other methods.

7 CONCLUSIONS

For the case of plants with a lag we have described a new method of uncertainty set propagation that appears to be more efficient than those in current use. There are in addition other advantages. It provides more information than most current methods, because the vertices and facets are both updated at every time step. Finally, because the computationally intensive steps are done using Boolean matrices, for moderately sized problems all calculations can be performed exactly over the rationals.