Navigation

Reachability

Basics of Reachability Analysis

Systems without disturbances

Consider a general continuous-time

(1)\dot{x}(t) = f(t, x, u),

or discrete-time dynamical system

(2)x(t+1) = f(t, x, u),

wherein t is time [1], x\in{\bf R}^n is the state, u\in{\bf R}^m is the control, and f is a measurable vector function taking values in {\bf R}^n. [2] The control values u(t, x(t)) are restricted to a closed compact control set {\mathcal U}(t)\subset{\bf R}^m. An open-loop control does not depend on the state, u=u(t); for a closed-loop control, u=u(t, x(t)).

Definition. The (forward) reach set {\mathcal X}(t, t_0, x_0) at time t>t_0 from the initial position (t_0, x_0) is the set of all states x(t) reachable at time t by system (1), or (2), with x(t_0)=x_0 through all possible controls u(\tau, x(\tau))\in{\mathcal U}(\tau), t_0\leqslant\tau< t. For a given set of initial states {\mathcal X}_0, the reach set {\mathcal X}(t, t_0, {\mathcal X}_0) is

(3){\mathcal X}(t, t_0, {\mathcal X}_0) = \bigcup_{x_0\in{\mathcal X}_0}{\mathcal X}(t, t_0, x_0).

Here are two facts about forward reach sets.

  1. {\mathcal X}(t, t_0, {\mathcal X}_0) is the same for open-loop and closed-loop control.

  2. {\mathcal X}(t, t_0, {\mathcal X}_0) satisfies the semigroup property,

    (4){\mathcal X}(t, t_0, {\mathcal X}_0) = {\mathcal X}(t, \tau, {\mathcal X}(\tau, t_0, {\mathcal X}_0)), \;\;\; t_0\leqslant\tau< t.

For linear systems

(5)f(t, x, u) = A(t)x(t) + B(t)u,

with matrices A(t) in {\bf R}^{n\times n} and B(t) in {\bf R}^{m\times n}. For continuous-time linear system the state transition matrix is

\dot{\Phi}(t, t_0) = A(t)\Phi(t, t_0), \Phi(t, t) = I,

which for constant A(t)\equiv A simplifies as

\Phi(t, t_0) = e^{A(t-t_0)} .

For discrete-time linear system the state transition matrix is

\Phi(t+1, t_0) = A(t)\Phi(t, t_0), \Phi(t, t) = I,

which for constant A(t)\equiv A simplifies as

\Phi(t, t_0) = A^{t-t_0} .

If the state transition matrix is invertible, \Phi^{-1}(t, t_0) = \Phi(t_0, t). The transition matrix is always invertible for continuous-time and for sampled discrete-time systems. However, if for some \tau, t_0\leqslant\tau<t, A(\tau) is degenerate (singular), \Phi(t, t_0)=\prod_{\tau=t_0}^{t-1}A(\tau), is also degenerate and cannot be inverted.

Following Cauchy’s formula, the reach set {\mathcal X}(t, t_0, {\mathcal X}_0) for a linear system can be expressed as

(6){\mathcal X}(t, t_0, {\mathcal X}_0) = \Phi(t, t_0){\mathcal X}_0 \oplus \int_{t_0}^t\Phi(t, \tau)B(\tau){\mathcal U}(\tau)d\tau

in continuous-time, and as

(7){\mathcal X}(t, t_0, {\mathcal X}_0) = \Phi(t, t_0){\mathcal X}_0 \oplus \sum_{\tau=t_0}^{t-1}\Phi(t, \tau+1)B(\tau){\mathcal U}(\tau)

in discrete-time case.

The operation ‘\oplus’ is the geometric sum, also known as Minkowski sum. [3] The geometric sum and linear (or affine) transformations preserve compactness and convexity. Hence, if the initial set {\mathcal X}_0 and the control sets {\mathcal U}(\tau), t_0\leqslant\tau<t, are compact and convex, so is the reach set {\mathcal X}(t, t_0, {\mathcal X}_0).

Definition. The backward reach set {\mathcal Y}(t_1, t, y_1) for the target position (t_1, y_1) is the set of all states y(t) for which there exists some control u(\tau, x(\tau))\in{\mathcal U}(\tau), t\leqslant\tau<t_1, that steers system (1), or (2) to the state y_1 at time t_1. For the target set {\mathcal Y}_1 at time t_1, the backward reach set {\mathcal Y}(t_1, t, {\mathcal Y}_1) is

(8){\mathcal Y}(t_1, t, {\mathcal Y}_1) = \bigcup_{y_1\in{\mathcal Y}_1}{\mathcal Y}(t_1, t, y_1).

The backward reach set {\mathcal Y}(t_1, t, {\mathcal Y}_1) is the largest weakly invariant set with respect to the target set {\mathcal Y}_1 and time values t and t_1. [4]

Remark. Backward reach set can be computed for continuous-time system only if the solution of (1) exists for t<t_1; and for discrete-time system only if the right hand side of (2) is invertible [5].

These two facts about the backward reach set {\mathcal Y} are similar to those for forward reach sets.

  1. {\mathcal Y}(t_1, t, {\mathcal Y}_1) is the same for open-loop and closed-loop control.

  2. {\mathcal Y}(t_1, t, {\mathcal Y}_1) satisfies the semigroup property,

    (9){\mathcal Y}(t_1, t, {\mathcal Y}_1) = {\mathcal Y}(\tau, t, {\mathcal Y}(t_1, \tau, {\mathcal Y}_1)), \;\;\; t\leqslant\tau< t_1.

For the linear system (5) the backward reach set can be expressed as

(10){\mathcal Y}(t_1, t, {\mathcal Y}_1) = \Phi(t, t_1){\mathcal Y}_1 \oplus \int_{t_1}^t\Phi(t, \tau)B(\tau){\mathcal U}(\tau)d\tau

in the continuous-time case, and as

(11){\mathcal Y}(t_1, t, {\mathcal Y}_1) = \Phi(t, t_1){\mathcal Y}_1 \oplus \sum_{\tau =t}^{t_1-1}-\Phi(t, \tau)B(\tau){\mathcal U}(\tau)

in discrete-time case. The last formula makes sense only for discrete-time linear systems with invertible state transition matrix. Degenerate discrete-time linear systems have unbounded backward reach sets and such sets cannot be computed with available software tools.

Just as in the case of forward reach set, the backward reach set of a linear system {\mathcal Y}(t_1, t, {\mathcal Y}_1) is compact and convex if the target set {\mathcal Y}_1 and the control sets {\mathcal U}(\tau), t\leqslant\tau<t_1, are compact and convex.

Remark. In the computer science literature the reach set is said to be the result of operator post, and the backward reach set is the result of operator pre. In the control literature the backward reach set is also called the solvability set.

Systems with disturbances

Consider the continuous-time dynamical system with disturbance

(12)\dot{x}(t) = f(t, x, u, v),

or the discrete-time dynamical system with disturbance

(13)x(t+1) = f(t, x, u, v),

in which we also have the disturbance input v\in{\bf R}^d with values v(t) restricted to a closed compact set {\mathcal V}(t)\subset{\bf R}^d.

In the presence of disturbances the open-loop reach set (OLRS) is different from the closed-loop reach set (CLRS).

Given the initial time t_0, the set of initial states {\mathcal X}_0, and terminal time t, there are two types of OLRS.

Definition. The maxmin open-loop reach set \overline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) is the set of all states x, such that for any disturbance v(\tau)\in{\mathcal V}(\tau), there exist an initial state x_0\in{\mathcal X}_0 and a control u(\tau)\in{\mathcal U}(\tau), t_0\leqslant\tau<t, that steers system (12) or (13) from x(t_0)=x_0 to x(t)=x.

Definition. The minmax open-loop reach set \underline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) is the set of all states x, such that there exists a control u(\tau)\in{\mathcal U}(\tau) that for all disturbances v(\tau)\in{\mathcal V}(\tau), t_0\leqslant\tau<t, assigns an initial state x_0\in{\mathcal X}_0 and steers system (12), or (13), from x(t_0)=x_0 to x(t)=x.

In the maxmin case the control is chosen after knowing the disturbance over the entire time interval [t_0, t], whereas in the minmax case the control is chosen before any knowledge of the disturbance. Consequently, the OLRS do not satisfy the semigroup property.

The terms ‘maxmin’ and ‘minmax’ come from the fact that \overline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) is the subzero level set of the value function

(14)\underline{V}(t, x) = \max_v\min_u\{{\bf dist}(x(t_0), {\mathcal X}_0) ~|~ x(t)=x, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; t_0\leqslant\tau<t\},

i.e., \overline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) = \{ x~|~\underline{V}(t, x) \leqslant0\}, and \underline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) is the subzero level set of the value function

(15)\overline{V}(t, x) = \min_u\max_v\{{\bf dist}(x(t_0), {\mathcal X}_0) ~|~ x(t)=x, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; t_0\leqslant\tau<t\},

in which {\bf dist}(\cdot, \cdot) denotes Hausdorff semidistance. [6] Since \underline{V}(t, x)\leqslant\overline{V}(t, x), \underline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0)\subseteq\overline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0).

Note that maxmin and minmax OLRS imply guarantees: these are states that can be reached no matter what the disturbance is, whether it is known in advance (maxmin case) or not (minmax case). The OLRS may be empty.

Fixing time instant \tau_1, t_0<\tau_1<t, define the piecewise maxmin open-loop reach set with one correction,

(16)\overline{{\mathcal X}}_{OL}^1(t, t_0, {\mathcal X}_0) = \overline{{\mathcal X}}_{OL}(t, \tau_1, \overline{{\mathcal X}}_{OL}(\tau_1, t_0, {\mathcal X}_0)),

and the piecewise minmax open-loop reach set with one correction,

(17)\underline{{\mathcal X}}_{OL}^1(t, t_0, {\mathcal X}_0) = \underline{{\mathcal X}}_{OL}(t, \tau_1, \underline{{\mathcal X}}_{OL}(\tau_1, t_0, {\mathcal X}_0)).

The piecewise maxmin OLRS \overline{{\mathcal X}}_{OL}^1(t, t_0, {\mathcal X}_0) is the subzero level set of the value function

(18)\underline{V}^1(t, x) = \max_v\min_u\{\underline{V}(\tau_1, x(\tau_1)) ~|~ x(t)=x, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; \tau_1\leqslant\tau<t\},

with V(\tau_1, x(\tau_1)) given by (14), which yields

\underline{V}^1(t, x) \geqslant\underline{V}(t, x),

and thus,

\overline{{\mathcal X}}_{OL}^1(t, t_0 {\mathcal X}_0) \subseteq \overline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) .

On the other hand, the piecewise minmax OLRS \underline{{\mathcal X}}_{OL}^1(t, t_0, {\mathcal X}_0) is the subzero level set of the value function

(19)\overline{V}^1(t, x) = \min_u\max_v\{\overline{V}(\tau_1, x(\tau_1)) ~|~ x(t)=x, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; \tau_1\leqslant\tau<t\},

with V(\tau_1, x(\tau_1)) given by (15), which yields

\overline{V}(t, x) \geqslant\overline{V}^1(t, x),

and thus,

\underline{{\mathcal X}}_{OL}(t, t_0 {\mathcal X}_0) \subseteq \underline{{\mathcal X}}_{OL}^1(t, t_0, {\mathcal X}_0) .

We can now recursively define piecewise maxmin and minmax OLRS with k corrections for t_0<\tau_1<\cdots<\tau_k<t. The maxmin piecewise OLRS with k corrections is

(20)\overline{{\mathcal X}}_{OL}^k(t, t_0, {\mathcal X}_0) = \overline{{\mathcal X}}_{OL}(t, \tau_k, \overline{{\mathcal X}}_{OL}^{k-1}(\tau_k, t_0, {\mathcal X}_0)),

which is the subzero level set of the corresponding value function

(21)\underline{V}^k(t, x) = \max_v\min_u\{\underline{V}^{k-1}(\tau_k, x(\tau_k)) ~|~ x(t)=x, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; \tau_k\leqslant\tau<t\}.

The minmax piecewise OLRS with k corrections is

(22)\underline{{\mathcal X}}_{OL}^k(t, t_0, {\mathcal X}_0) = \underline{{\mathcal X}}_{OL}(t, \tau_k, \underline{{\mathcal X}}_{OL}^{k-1}(\tau_k, t_0, {\mathcal X}_0)),

which is the subzero level set of the corresponding value function

(23)\begin{aligned} \overline{V}^k(t, x) = \min_u\max_v\{\overline{V}^{k-1}(\tau_k, x(\tau_k)) ~|~ x(t)=x, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; \tau_k\leqslant\tau<t\}. \end{aligned}

From (18), (19), (21) and (23) it follows that

\underline{V}(t, x) \leqslant\underline{V}^1(t, x)\leqslant\cdots \leqslant\underline{V}^k(t, x) \leqslant\overline{V}^k(t, x) \leqslant\cdots \leqslant\overline{V}^1(t, x) \leqslant\overline{V}(t, x) .

Hence,

(24)\begin{aligned} &&\underline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) \subseteq \underline{{\mathcal X}}_{OL}^1(t, t_0, {\mathcal X}_0) \subseteq \cdots \subseteq \underline{{\mathcal X}}_{OL}^k(t, t_0, {\mathcal X}_0) \subseteq \nonumber \\ &&\overline{{\mathcal X}}_{OL}^k(t, t_0, {\mathcal X}_0) \subseteq \cdots \subseteq \overline{{\mathcal X}}_{OL}^1(t, t_0, {\mathcal X}_0) \subseteq \overline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) . \end{aligned}

We call

(25)\overline{{\mathcal X}}_{CL}(t, t_0, {\mathcal X}_0) = \overline{{\mathcal X}}_{OL}^k(t, t_0, {\mathcal X}_0), \;\; k = \left\{\begin{array}{ll} \infty & \mbox{ for continuous-time system}\\ t-t_0-1 & \mbox{ for discrete-time system}\end{array}\right.

the maxmin closed-loop reach set of system (12) or (13) at time t, and we call

(26)\underline{{\mathcal X}}_{CL}(t, t_0, {\mathcal X}_0) = \underline{{\mathcal X}}_{OL}^k(t, t_0, {\mathcal X}_0), \;\; k = \left\{\begin{array}{ll} \infty & \mbox{ for continuous-time system}\\ t-t_0-1 & \mbox{ for discrete-time system}\end{array}\right.

the minmax closed-loop reach set of system (12) or (13) at time t.

Definition. Given initial time t_0 and the set of initial states {\mathcal X}_0, the maxmin CLRS \overline{{\mathcal X}}_{CL}(t, t_0, {\mathcal X}_0) of system (12) or (13) at time t>t_0, is the set of all states x, for each of which and for every disturbance v(\tau)\in{\mathcal V}(\tau), there exist an initial state x_0\in{\mathcal X}_0 and a control u(\tau, x(\tau))\in{\mathcal U}(\tau), such that the trajectory x(\tau | v(\tau), u(\tau, x(\tau))) satisfying x(t_0) = x_0 and

\dot{x}(\tau | v(\tau), u(\tau, x(\tau))) \in f(\tau, x(\tau), u(\tau, x(\tau)), v(\tau))

in the continuous-time case, or

x(\tau+1 | v(\tau), u(\tau, x(\tau))) \in f(\tau, x(\tau), u(\tau, x(\tau)), v(\tau))

in the discrete-time case, with t_0\leqslant\tau<t, is such that x(t)=x.

Definition. Given initial time t_0 and the set of initial states {\mathcal X}_0, the maxmin CLRS \underline{{\mathcal X}}_{CL}(t, t_0, {\mathcal X}_0) of system (12) or (13), at time t>t_0, is the set of all states x, for each of which there exists a control u(\tau, x(\tau))\in{\mathcal U}(\tau), and for every disturbance v(\tau)\in{\mathcal V}(\tau) there exists an initial state x_0\in{\mathcal X}_0, such that the trajectory x(\tau, v(\tau) | u(\tau, x(\tau))) satisfying x(t_0) = x_0 and

\dot{x}(\tau, v(\tau) | u(\tau, x(\tau))) \in f(\tau, x(\tau), u(\tau, x(\tau)), v(\tau))

in the continuous-time case, or

x(\tau+1, v(\tau) | u(\tau, x(\tau))) \in f(\tau, x(\tau), u(\tau, x(\tau)), v(\tau))

in the discrete-time case, with t_0\leqslant\tau<t, is such that x(t)=x. By construction, both maxmin and minmax CLRS satisfy the semigroup property (4).

For some classes of dynamical systems and some types of constraints on initial conditions, controls and disturbances, the maxmin and minmax CLRS may coincide. This is the case for continuous-time linear systems with convex compact bounds on the initial set, controls and disturbances under the condition that the initial set {\mathcal X}_0 is large enough to ensure that {\mathcal X}(t_0+\epsilon, t_0, {\mathcal X}_0) is nonempty for some small \epsilon>0.

Consider the linear system case,

(27)f(t, x, u) = A(t)x(t) + B(t)u + G(t)v,

where A(t) and B(t) are as in (5), and G(t) takes its values in {\bf R}^d.

The maxmin OLRS for the continuous-time linear system can be expressed through set valued integrals,

(28)\begin{array}{l} \overline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) = \\ \left(\Phi(t, t_0){\mathcal X}_0 \oplus \int_{t_0}^t\Phi(t, \tau)B(\tau){\mathcal U}(\tau)d\tau\right) \dot{-} \\ \int_{t_0}^t\Phi(t, \tau)(-G(\tau)){\mathcal V}(\tau)d\tau, \end{array}

and for discrete-time linear system through set-valued sums,

(29)\begin{array}{l} \overline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) = \\ \left(\Phi(t, t_0){\mathcal X}_0 \oplus \sum_{\tau=t_0}^{t-1}\Phi(t, \tau+1)B(\tau){\mathcal U}(\tau)\right) \dot{-} \\ \sum_{\tau=t_0}^{t-1}\Phi(t, \tau+1)(-G(\tau)){\mathcal V}(\tau). \end{array}

Similarly, the minmax OLRS for the continuous-time linear system is

(30)\begin{array}{l} \underline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) = \\ \left(\Phi(t, t_0){\mathcal X}_0 \dot{-} \int_{t_0}^t\Phi(t, \tau)(-G(\tau)){\mathcal V}(\tau)d\tau\right) \oplus \\ \int_{t_0}^t\Phi(t, \tau)B(\tau){\mathcal U}(\tau)d\tau, \end{array}

and for the discrete-time linear system it is

(31)\begin{array}{l} \underline{{\mathcal X}}_{OL}(t, t_0, {\mathcal X}_0) = \\ \left(\Phi(t, t_0){\mathcal X}_0 \dot{-} \sum_{\tau=t_0}^{t-1}\Phi(t, \tau+1)(-G(\tau)){\mathcal V}(\tau)\right) \oplus \\ \sum_{\tau=t_0}^{t-1}\Phi(t, \tau+1)B(\tau){\mathcal U}(\tau). \end{array}

The operation ‘\dot{-}’ is geometric difference, also known as Minkowski difference. [7]

Now consider the piecewise OLRS with k corrections. Expression (20) translates into

(32)\begin{array}{l} \overline{{\mathcal X}}_{OL}^k(t, t_0, {\mathcal X}_0) = \\ \left(\Phi(t, \tau_k)\overline{{\mathcal X}}_{OL}^{k-1}(\tau_k, t_0, {\mathcal X}_0) \oplus \int_{\tau_k}^t\Phi(t, \tau)B(\tau){\mathcal U}(\tau)d\tau\right) \dot{-} \\ \int_{\tau_k}^t\Phi(t, \tau)(-G(\tau)){\mathcal V}(\tau)d\tau, \end{array}

in the continuous-time case, and for the discrete-time case into

(33)\begin{array}{l} \overline{{\mathcal X}}_{OL}^k(t, t_0, {\mathcal X}_0) = \\ \left(\Phi(t, \tau_k)\overline{{\mathcal X}}_{OL}^{k-1}(\tau_k, t_0, {\mathcal X}_0) \oplus \sum_{\tau=\tau_k}^{t-1}\Phi(t, \tau+1)B(\tau){\mathcal U}(\tau)\right) \dot{-} \\ \sum_{\tau=\tau_k}^{t-1}\Phi(t, \tau+1)(-G(\tau)){\mathcal V}(\tau). \end{array}

Expression (22) translates into

(34)\begin{array}{l} \underline{{\mathcal X}}_{OL}^k(t, t_0, {\mathcal X}_0) = \\ \left(\Phi(t, \tau_k)\underline{{\mathcal X}}_{OL}^{k-1}(t, t_0, {\mathcal X}_0) \dot{-} \int_{\tau_k}^t\Phi(t, \tau)(-G(\tau)){\mathcal V}(\tau)d\tau\right) \oplus \\ \int_{\tau_k}^t\Phi(t, \tau)B(\tau){\mathcal U}(\tau)d\tau, \end{array}

in the continuous-time case, and for the discrete-time case into

(35)\begin{array}{l} \underline{{\mathcal X}}_{OL}^k(t, t_0, {\mathcal X}_0) = \\ \left(\Phi(t, \tau_k)\underline{{\mathcal X}}_{OL}^{k-1}(\tau_k, t_0, {\mathcal X}_0) \dot{-} \sum_{\tau=\tau_k}^{t-1}\Phi(t, \tau+1)(-G(\tau)){\mathcal V}(\tau)\right) \oplus \\ \sum_{\tau=\tau_k}^{t-1}\Phi(t, \tau+1)B(\tau){\mathcal U}(\tau). \end{array}

Since for any {\mathcal W}_1, {\mathcal W}_2, {\mathcal W}_3 \subseteq {\bf R}^n it is true that

({\mathcal W}_1 \dot{-} {\mathcal W}_2) \oplus {\mathcal W}_3 = ({\mathcal W}_1 \oplus {\mathcal W}_3) \dot{-} ({\mathcal W}_2 \oplus {\mathcal W}_3) \subseteq ({\mathcal W}_1 \oplus {\mathcal W}_3) \dot{-} {\mathcal W}_2,

from (32), (34) and from (33), (35), it is clear that (24) is true. For linear systems, if the initial set {\mathcal X}_0, control bounds {\mathcal U}(\tau) and disturbance bounds {\mathcal V}(\tau), t_0\leqslant\tau<t, are compact and convex, the CLRS \overline{{\mathcal X}}_{CL}(t, t_0, {\mathcal X}_0) and \underline{{\mathcal X}}_{CL}(t, t_0, {\mathcal X}_0) are compact and convex, provided they are nonempty. For continuous-time linear systems, \overline{{\mathcal X}}_{CL}(t, t_0, {\mathcal X}_0) = \underline{{\mathcal X}}_{CL}(t, t_0, {\mathcal X}_0) = {\mathcal X}_{CL}(t, t_0, {\mathcal X}_0).

Just as for forward reach sets, the backward reach sets can be open-loop (OLBRS) or closed-loop (CLBRS).

Definition. Given the terminal time t_1 and target set {\mathcal Y}_1, the maxmin open-loop backward reach set \overline{{\mathcal Y}}_{OL}(t_1, t, {\mathcal Y}_1) of system (12) or (13) at time t<t_1, is the set of all y, such that for any disturbance v(\tau)\in{\mathcal V}(\tau) there exists a terminal state y_1\in{\mathcal Y}_1 and control u(\tau)\in{\mathcal U}(\tau), t\leqslant\tau<t_1, which steers the system from y(t)=y to y(t_1)=y_1.

\overline{{\mathcal Y}}_{OL}(t_1, t, {\mathcal Y}_1) is the subzero level set of the value function

(36)\begin{aligned} \underline{V}_b(t, y) = \max_v\min_u\{{\bf dist}(y(t_1), {\mathcal Y}_1) ~|~ y(t)=y, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; t\leqslant\tau<t_1\}, \end{aligned}

Definition. Given the terminal time t_1 and target set {\mathcal Y}_1, the minmax open-loop backward reach set \underline{{\mathcal Y}}_{OL}(t_1, t, {\mathcal Y}_1) of system (12) or (13) at time t<t_1, is the set of all y, such that there exists a control u(\tau)\in{\mathcal U}(\tau) that for all disturbances v(\tau\in{\mathcal V}(\tau), t\leqslant\tau<t_1, assigns a terminal state y_1\in{\mathcal Y}_1 and steers the system from y(t)=y to y(t_1)=y_1. \underline{{\mathcal Y}}_{OL}(t_1, t, {\mathcal Y}_1) is the subzero level set of the value function

(37)\begin{aligned} \overline{V}_b(t, y) = \min_u\max_v\{{\bf dist}(y(t_1), {\mathcal Y}_1) ~|~ y(t)=y, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; t\leqslant\tau<t_1\}, \end{aligned}

Remark. The backward reach set can be computed for a continuous-time system only if the solution of (12) exists for t<t_1, and for a discrete-time system only if the right hand side of (13) is invertible.

Similarly to the forward reachability case, we construct piecewise OLBRS with one correction at time \tau_1, t<\tau_1<t_1. The piecewise maxmin OLBRS with one correction is

(38)\overline{{\mathcal Y}}_{OL}^1(t_1, t, {\mathcal Y}_1) = \overline{{\mathcal Y}}_{OL}(\tau_1, t, \overline{{\mathcal Y}}_{OL}(t_1, \tau_1, {\mathcal Y}_1)),

and it is the subzero level set of the function

(39)\begin{aligned} \underline{V}^1_b(t, y) = \max_v\min_u\{\underline{V}_b(\tau_1, y(\tau_1)) ~|~ y(t)=y, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; t\leqslant\tau<\tau_1\}. \end{aligned}

The piecewise minmax OLBRS with one correction is

(40)\underline{{\mathcal Y}}_{OL}^1(t_1, t, {\mathcal Y}_1) = \underline{{\mathcal Y}}_{OL}(\tau_1, t, \underline{{\mathcal Y}}_{OL}(t_1, \tau_1, {\mathcal Y}_1)),

and it is the subzero level set of the function

(41)\begin{aligned} \overline{V}^1_b(t, y) = \min_u\max_v\{\overline{V}_b(\tau_1, y(\tau_1)) ~|~ y(t)=y, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; t\leqslant\tau<\tau_1\}, \end{aligned}

Recursively define maxmin and minmax OLBRS with k corrections for t<\tau_k<\cdots<\tau_1<t_1. The maxmin OLBRS with k corrections is

(42)\overline{{\mathcal Y}}_{OL}^k(t_1, t, {\mathcal Y}_1) = \overline{{\mathcal Y}}_{OL}(\tau_k, t, \overline{{\mathcal Y}}_{OL}^{k-1}(t_1, \tau_k, {\mathcal Y}_1)),

which is the subzero level set of function

(43)\begin{aligned} \underline{V}^k_b(t, y) = \max_v\min_u\{\underline{V}^{k-1}_b(\tau_k, y(\tau_k)) ~|~ y(t)=y, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; t\leqslant\tau<\tau_k\}. \end{aligned}

The minmax OLBRS with k corrections is

(44)\underline{{\mathcal Y}}_{OL}^k(t_1, t, {\mathcal Y}_1) = \underline{{\mathcal Y}}_{OL}(\tau_k, t, \underline{{\mathcal Y}}_{OL}^{k-1}(t_1, \tau_k, {\mathcal Y}_1)),

which is the subzero level set of the function

(45)\begin{aligned} \overline{V}^k_b(t, y) = \min_u\max_v\{\overline{V}^{k-1}_b(\tau_k, y(\tau_k)) ~|~ y(t)=y, \; u(\tau)\in{\mathcal U}(\tau), \; v(\tau)\in{\mathcal V}(\tau), \; t\leqslant\tau<\tau_k\}, \end{aligned}

From (39), (41), (43) and (45) it follows that

\underline{V}_b(t, y) \leqslant\underline{V}^1_b(t, y)\leqslant\cdots \leqslant\underline{V}^k_b(t, y) \leqslant\overline{V}^k_b(t, y) \leqslant\cdots \leqslant\overline{V}^1_b(t, y) \leqslant\overline{V}_b(t, y) .

Hence,

(46)\begin{aligned} &&\underline{{\mathcal Y}}_{OL}(t_1, t, {\mathcal Y}_1) \subseteq \underline{{\mathcal Y}}_{OL}^1(t_1, t, {\mathcal Y}_1) \subseteq \cdots \subseteq \underline{{\mathcal Y}}_{OL}^k(t_1, t, {\mathcal Y}_1) \subseteq \nonumber \\ &&\overline{{\mathcal Y}}_{OL}^k(t_1, t, {\mathcal Y}_1) \subseteq \cdots \subseteq \overline{{\mathcal Y}}_{OL}^1(t_1, t, {\mathcal Y}_1) \subseteq \overline{{\mathcal Y}}_{OL}(t_1, t, {\mathcal Y}_1) . \end{aligned}

We say that

(47)\overline{{\mathcal Y}}_{CL}(t_1, t, {\mathcal Y}_1) = \overline{{\mathcal Y}}_{OL}^k(t_1, t, {\mathcal Y}_1), \;\; k = \left\{\begin{array}{ll} \infty & \mbox{ for continuous-time system}\\ t_1-t-1 & \mbox{ for discrete-time system}\end{array}\right.

is the maxmin closed-loop backward reach set of system (12) or (13) at time t.

We say that

(48)\underline{{\mathcal Y}}_{CL}(t_1, t, {\mathcal Y}_1) = \underline{{\mathcal Y}}_{OL}^k(t_1, t, {\mathcal Y}_1), \;\; k = \left\{\begin{array}{ll} \infty & \mbox{ for continuous-time system}\\ t_1-t-1 & \mbox{ for discrete-time system}\end{array}\right.

is the minmax closed-loop backward reach set of system (12) or (13) at time t.

Definition. Given the terminal time t_1 and target set {\mathcal Y}_1, the maxmin CLBRS \overline{{\mathcal Y}}_{CL}(t_1, t, {\mathcal Y}_1) of system (12) or (13) at time t<t_1, is the set of all states y, for each of which for every disturbance v(\tau)\in{\mathcal V}(\tau) there exists terminal state y_1\in{\mathcal Y}_1 and control u(\tau, y(\tau))\in{\mathcal U}(\tau) that assigns trajectory y(\tau, | v(\tau), u(\tau, y(\tau))) satisfying

\dot{y}(\tau | v(\tau), u(\tau, y(\tau))) \in f(\tau, y(\tau), u(\tau, y(\tau)), v(\tau))

in continuous-time case, or

y(\tau+1 | v(\tau), u(\tau, y(\tau))) \in f(\tau, y(\tau), u(\tau, y(\tau)), v(\tau))

in discrete-time case, with t\leqslant\tau<t_1, such that y(t) = y and y(t_1)=y_1.

Definition. Given the terminal time t_1 and target set {\mathcal Y}_1, the minmax CLBRS \underline{{\mathcal Y}}_{CL}(t_1, t, {\mathcal Y}_1) of system (12) or (13) at time t<t_1, is the set of all states y, for each of which there exists control u(\tau, y(\tau))\in{\mathcal U}(\tau) that for every disturbance v(\tau)\in{\mathcal V}(\tau) assigns terminal state y_1\in{\mathcal Y}_1 and trajectory y(\tau, v(\tau) | u(\tau, y(\tau))) satisfying

\dot{y}(\tau, v(\tau) | u(\tau, y(\tau))) \in f(\tau, y(\tau), u(\tau, y(\tau)), v(\tau))

in the continuous-time case, or

y(\tau+1, v(\tau) | u(\tau, y(\tau))) \in f(\tau, y(\tau), u(\tau, y(\tau)), v(\tau))

in the discrete-time case, with t\leqslant\tau<t_1, such that y(t) = y and y(t_1)=y_1.

Both maxmin and minmax CLBRS satisfy the semigroup property (9).

The maxmin OLBRS for the continuous-time linear system can be expressed through set valued integrals,

(49)\begin{array}{l} \overline{{\mathcal Y}}_{OL}(t_1, t, {\mathcal Y}_1) = \\ \left(\Phi(t, t_1){\mathcal Y}_1 \oplus \int_{t_1}^t\Phi(t, \tau)B(\tau){\mathcal U}(\tau)d\tau\right) \dot{-} \\ \int_{t}^{t_1}\Phi(t, \tau)G(\tau){\mathcal V}(\tau)d\tau, \end{array}

and for the discrete-time linear system through set-valued sums,

(50)\begin{array}{l} \overline{{\mathcal Y}}_{OL}(t_1, t, {\mathcal Y}_1) = \\ \left(\Phi(t, t_1){\mathcal Y}_1 \oplus \sum_{\tau=t}^{t_1-1}-\Phi(t, \tau+1)B(\tau){\mathcal U}(\tau)\right) \dot{-} \\ \sum_{\tau=t}^{t_1-1}\Phi(t, \tau+1)G(\tau){\mathcal V}(\tau). \end{array}

Similarly, the minmax OLBRS for the continuous-time linear system is

(51)\begin{array}{l} \underline{{\mathcal Y}}_{OL}(t_1, t, {\mathcal Y}_1) = \\ \left(\Phi(t, t_1){\mathcal Y}_1 \dot{-} \int_{t}^{t_1}\Phi(t, \tau)G(\tau){\mathcal V}(\tau)d\tau\right) \oplus \\ \int_{t_1}^{t}\Phi(t, \tau)B(\tau){\mathcal U}(\tau)d\tau, \end{array}

and for the discrete-time linear system it is

(52)\begin{array}{l} \underline{{\mathcal Y}}_{OL}(t_1, t, {\mathcal Y}_1) = \\ \left(\Phi(t, t_1){\mathcal Y}_1 \dot{-} \sum_{\tau=t}^{t_1-1}\Phi(t, \tau+1)G(\tau){\mathcal V}(\tau)\right) \oplus \\ \sum_{\tau=t}^{t_1-1}-\Phi(t, \tau+1)B(\tau){\mathcal U}(\tau). \end{array}

Now consider piecewise OLBRS with k corrections. Expression (42) translates into

(53)\begin{array}{l} \overline{{\mathcal Y}}_{OL}^k(t_1, t, {\mathcal Y}_1) = \\ \left(\Phi(t, \tau_k)\overline{{\mathcal Y}}_{OL}^{k-1}(t_1, \tau_k, {\mathcal Y}_1) \oplus \int_{\tau_k}^t\Phi(t, \tau)B(\tau){\mathcal U}(\tau)d\tau\right) \dot{-} \\ \int^{\tau_k}_t\Phi(t, \tau)G(\tau){\mathcal V}(\tau)d\tau, \end{array}

in the continuous-time case, and for the discrete-time case into

(54)\begin{array}{l} \overline{{\mathcal Y}}_{OL}^k(t_1, t, {\mathcal Y}_1) = \\ \left(\Phi(t, \tau_k)\overline{{\mathcal Y}}_{OL}^{k-1}(t_1, \tau_k, {\mathcal Y}_1) \oplus \sum_{\tau=t}^{\tau_k-1}-\Phi(t, \tau+1)B(\tau){\mathcal U}(\tau)\right) \dot{-} \\ \sum_{\tau=t}^{\tau_k-1}\Phi(t, \tau+1)G(\tau){\mathcal V}(\tau). \end{array}

Expression (44) translates into

(55)\begin{array}{l} \underline{{\mathcal Y}}_{OL}^k(t_1, t, {\mathcal Y}_1) = \\ \left(\Phi(t, \tau_k)\overline{{\mathcal Y}}_{OL}^{k-1}(t_1, \tau_k, {\mathcal Y}_1) \dot{-} \int^{\tau_k}_t\Phi(t, \tau)G(\tau){\mathcal V}(\tau)d\tau\right) \oplus \\ \int_{\tau_k}^t\Phi(t, \tau)B(\tau){\mathcal U}(\tau)d\tau, \end{array}

in the continuous-time case, and for the discrete-time case into

(56)\begin{array}{l} \underline{{\mathcal Y}}_{OL}^k(t_1, t, {\mathcal Y}_1) = \\ (\Phi(t, \tau_k)\overline{{\mathcal Y}}_{OL}^{k-1}(t_1, \tau_k, {\mathcal Y}_1) \dot{-} \sum_{\tau=t}^{\tau_k-1}\Phi(t, \tau+1)G(\tau){\mathcal V}(\tau)) \oplus \\ \sum_{\tau=t}^{\tau_k-1}-\Phi(t, \tau+1)B(\tau){\mathcal U}(\tau). \end{array}

For continuous-time linear systems \overline{{\mathcal Y}}_{CL}(t_1, t, {\mathcal Y}_1) = \underline{{\mathcal Y}}_{CL}(t_1, t, {\mathcal Y}_1) = {\mathcal Y}_{CL}(t_1, t, {\mathcal Y}_1) under the condition that the target set {\mathcal Y}_1 is large enough to ensure that \underline{{\mathcal Y}}_{CL}(t_1, t_1-\epsilon, {\mathcal Y}_1) is nonempty for some small \epsilon>0.

Computation of backward reach sets for discrete-time linear systems makes sense only if the state transition matrix \Phi(t_1, t) is invertible.

If the target set {\mathcal Y}_1, control sets {\mathcal U}(\tau) and disturbance sets {\mathcal V}(\tau), t\leqslant\tau<t_1, are compact and convex, then CLBRS \overline{{\mathcal Y}}_{CL}(t_1, t, {\mathcal Y}_1) and \underline{{\mathcal Y}}_{CL}(t_1, t, {\mathcal Y}_1) are compact and convex, if they are nonempty.

Reachability problem

Reachability analysis is concerned with the computation of the forward {\mathcal X}(t, t_0, {\mathcal X}_0) and backward {\mathcal Y}(t_1, t, {\mathcal Y}_1) reach sets (the reach sets may be maxmin or minmax) in a way that can effectively meet requests like the following:

  1. For the given time interval [t_0, t], determine whether the system can be steered into the given target set {\mathcal Y}_1. In other words, is the set {\mathcal Y}_1\cap\bigcup_{t_0 \leqslant\tau\leqslant t}{\mathcal X}(\tau, t_0, {\mathcal X}_0) nonempty? And if the answer is ‘yes’, find a control that steers the system to the target set (or avoids the target set). [8]

  2. If the target set {\mathcal Y}_1 is reachable from the given initial condition \{t_0, {\mathcal X}_0\} in the time interval [t_0, t], find the shortest time to reach {\mathcal Y}_1,

    \arg\min_{\tau} \{{\mathcal X}(\tau,t_0,{\mathcal X}_0)\cap{\mathcal Y}_1\neq\emptyset ~|~ t_0\leqslant\tau\leqslant t\}.

  3. Given the terminal time t_1, target set {\mathcal Y}_1 and time t<t_1 find the set of states starting at time t from which the system can reach {\mathcal Y}_1 within time interval [t, t_1]. In other words, find \bigcup_{t\leqslant\tau<t_1}{\mathcal Y}(t_1, \tau, {\mathcal Y}_1).

  4. Find a closed-loop control that steers a system with disturbances to the given target set in given time.

  5. Graphically display the projection of the reach set along any specified two- or three-dimensional subspace.

For linear systems, if the initial set {\mathcal X}_0, target set {\mathcal Y}_1, control bounds {\mathcal U}(\cdot) and disturbance bounds {\mathcal V}(\cdot) are compact and convex, so are the forward {\mathcal X}(t, t_0, {\mathcal X}_0) and backward {\mathcal Y}(t_1, t, {\mathcal Y}_1) reach sets. Hence reachability analysis requires the computationally effective manipulation of convex sets, and performing the set-valued operations of unions, intersections, geometric sums and differences.

Existing reach set computation tools can deal reliably only with linear systems with convex constraints. A claim that certain tool or method can be used effectively for nonlinear systems must be treated with caution, and the first question to ask is for what class of nonlinear systems and with what limit on the state space dimension does this tool work? Some “reachability methods for nonlinear systems” reduce to the local linearization of a system followed by the use of well-tested techniques for linear system reach set computation. Thus these approaches in fact use reachability methods for linear systems.

Ellipsoidal Method

Continuous-time systems

Consider the system

(57)\dot{x}(t) = A(t)x(t) + B(t)u + G(t)v,

in which x\in{\bf R}^n is the state, u\in{\bf R}^m is the control and v\in{\bf R}^d is the disturbance. A(t), B(t) and G(t) are continuous and take their values in {\bf R}^{n\times n}, {\bf R}^{n\times m} and {\bf R}^{n\times d} respectively. Control u(t,x(t)) and disturbance v(t) are measurable functions restricted by ellipsoidal constraints: u(t,x(t)) \in {\mathcal E}(p(t), P(t)) and v(t) \in {\mathcal E}(q(t), Q(t)). The set of initial states at initial time t_0 is assumed to be the ellipsoid {\mathcal E}(x_0,X_0).

The reach sets for systems with disturbances computed by the Ellipsoidal Toolbox are CLRS. Henceforth, when describing backward reachability, reach sets refer to CLRS or CLBRS. Recall that for continuous-time linear systems maxmin and minmax CLRS coincide, and the same is true for maxmin and minmax CLBRS.

If the matrix Q(\cdot)=0, the system (57) becomes an ordinary affine system with known v(\cdot)=q(\cdot). If G(\cdot) = 0, the system becomes linear. For these two cases (Q(\cdot)=0 or G(\cdot)=0) the reach set is as given in definition (3), and so the reach set will be denoted as {\mathcal X}_{CL}(t, t_0, {\mathcal E}(x_0, X_0)) = {\mathcal X}(t, t_0, {\mathcal E}(x_0,X_0)).

The reach set {\mathcal X}(t,t_0,{\mathcal E}(x_0,X_0)) is a symmetric compact convex set, whose center evolves in time according to

(58)\dot{x}_c(t) = A(t)x_c(t) + B(t)p(t) + G(t)q(t), \;\;\; x_c(t_0)=x_0.

Fix a vector l_0\in{\bf R}^n, and consider the solution l(t) of the adjoint equation

(59)\dot{l}(t) = -A^T(t)l(t), \;\;\; l(t_0) = l_0,

which is equivalent to

l(t) = \Phi^T(t_0, t)l_0.

If the reach set {\mathcal X}(t, t_0, {\mathcal E}(x_0,X_0)) is nonempty, there exist tight external and tight internal approximating ellipsoids {\mathcal E}(x_c(t), X^+_l(t)) and {\mathcal E}(x_c(t), X^-_l(t)), respectively, such that

(60){\mathcal E}(x_c(t), X^-_l(t))\subseteq{\mathcal X}(t,t_0,{\mathcal E}(x_0,X_0)) \subseteq {\mathcal E}(x_c(t), X^+_l(t)),

and

(61)\rho(l(t) ~|~ {\mathcal E}(x_c(t), X^-_l(t))) = \rho(l(t) ~|~ {\mathcal X}(t, t_0, {\mathcal E}(x_0,X_0))) = \rho(l(t) ~|~ {\mathcal E}(x_c(t), X^+_l(t))) .

The equation for the shape matrix of the external ellipsoid is

(62)\dot{X}^+_l(t) & = A(t)X^+_l(t) + X^+_l(t)A^T(t) +\nonumber \\ &\pi_l(t)X^+_l(t) + \frac{1}{\pi_l(t)}B(t)P(t)B^T(t) -\nonumber \\ & (X_l^{+}(t))^{1/2}S_l(t)(G(t)Q(t)G^T(t))^{1/2} \nonumber -\\ & (G(t)Q(t)G^T(t))^{1/2}S_l^T(t)(X_l^{+}(t))^{1/2}, \\

(63)X^+_l(t_0) =X_0,

in which

\pi_l(t) = \frac{\langle l(t), B(t)P(t)B^T(t)l(t)\rangle^{1/2}}{\langle l(t), X^+_l(t)l(t)\rangle^{1/2}},

and the orthogonal matrix S_l(t) (S_l(t)S_l^T(t) = I) is determined by the equation

S_l(t)(G(t)Q(t)G^T(t))^{1/2}l(t) = \frac{\langle l(t), G(t)Q(t)G^T(t)l(t)\rangle^{1/2}}{\langle l(t), X_l^+(t)l(t)\rangle^{1/2}}(X_l^{+}(t))^{1/2}l(t).

In the presence of disturbance, if the reach set is empty, the matrix X^+_l(t) becomes sign indefinite. For a system without disturbance, the terms containing G(t) and Q(t) vanish from the equation (62).

The equation for the shape matrix of the internal ellipsoid is

(64)\dot{X}^-_l(t) & = A(t)X^-_l(t) + X^-_l(t)A^T(t) +\nonumber \\ & (X_l^{-}(t))^{1/2}T_l(t)(B(t)P(t)B^T(t))^{1/2} +\nonumber \\ & (B(t)P(t)B^T(t))^{1/2}T_l^T(t)(X_l^{-}(t))^{1/2} -\nonumber \\ & \eta_l(t)X^-_l(t) - \frac{1}{\eta_l(t)}G(t)Q(t)G^T(t), \\

(65)X^-_l(t_0) = X_0,

in which

\eta_l(t) = \frac{\langle l(t), G(t)Q(t)G^T(t)l(t)\rangle^{1/2}}{\langle l(t), X^+_l(t)l(t)\rangle^{1/2}},

and the orthogonal matrix T_l(t) is determined by the equation

T_l(t)(B(t)P(t)B^T(t))^{1/2}l(t) = \frac{\langle l(t), B(t)P(t)B^T(t)l(t)\rangle^{1/2}}{\langle l(t), X_l^-(t)l(t)\rangle^{1/2}}(X_l^{-}(t))^{1/2}l(t).

Similarly to the external case, the terms containing G(t) and Q(t) vanish from the equation (64) for a system without disturbance.

The point where the external and internal ellipsoids touch the boundary of the reach set is given by

x_l^*(t) = x_c(t) + \frac{X^+_l(t)l(t)}{\langle l(t), X^+_l(t)l(t)\rangle^{1/2}} .

The boundary points x^*_l(t) form trajectories, which we call extremal trajectories. Due to the nonsingular nature of the state transition matrix \Phi(t,t_0), every boundary point of the reach set belongs to an extremal trajectory. To follow an extremal trajectory specified by parameter l_0, the system has to start at time t_0 at initial state

(66)x^0_l = x_0 + \frac{X_0l_0}{\langle l_0,X_0l_0\rangle^{1/2}}.

In the absence of disturbances, the open-loop control

(67)u_l(t) = p(t) + \frac{P(t)B^T(t)l(t)}{\langle l(t), B(t)P(t)B^T(t)l(t)\rangle^{1/2}}.

steers the system along the extremal trajectory defined by the vector l_0. When a disturbance is present, this control keeps the system on an extremal trajectory if and only if the disturbance plays against the control always taking its extreme values.

Expressions (60) and (61) lead to the following fact,

\bigcup_{\langle l_0,l_0\rangle=1}{\mathcal E}(x_c(t),X^-_l(t)) = {\mathcal X}(t,t_0,{\mathcal E}(x_0,X_0)) = \bigcap_{\langle l_0,l_0\rangle=1}{\mathcal E}(x_c(t),X^+_l(t)).

In practice this means that the more values of l_0 we use to compute X^+_l(t) and X^-_l(t), the better will be our approximation.

Analogous results hold for the backward reach set.

Given the terminal time t_1 and ellipsoidal target set {\mathcal E}_(y_1,Y_1), the CLBRS {\mathcal Y}_{CL}(t_1, t, {\mathcal Y}_1)={\mathcal Y}(t_1, t, {\mathcal Y}_1), t<t_1, if it is nonempty, is a symmetric compact convex set whose center is governed by

(68)y_c(t) = Ay_c(t) + B(t)p(t) + G(t)q(t), \;\;\; y_c(t_1) = y_1.

Fix a vector l_1\in{\bf R}^n, and consider

(69)l(t) = \Phi(t_1, t)^Tl_1 .

If the backward reach set {\mathcal Y}(t_1, t, {\mathcal E}(y_1,Y_1)) is nonempty, there exist tight external and tight internal approximating ellipsoids {\mathcal E}(y_c(t), Y^+_l(t)) and {\mathcal E}(y_c(t), Y^-_l(t)) respectively, such that

(70){\mathcal E}(y_c(t), Y^-_l(t))\subseteq{\mathcal Y}(t_1,t,{\mathcal E}(y_1,Y_1)) \subseteq {\mathcal E}(y_c(t), Y^+_l(t)),

and

(71)\rho(l(t) ~|~ {\mathcal E}(y_c(t), Y^-_l(t))) = \rho(l(t) ~|~ {\mathcal Y}(t_1, t, {\mathcal E}(y_0,Y_0))) = \rho(l(t) ~|~ {\mathcal E}(y_c(t), Y^+_l(t))) .

The equation for the shape matrix of the external ellipsoid is

(72)\dot{Y}^+_l(t) & = A(t)Y^+_l(t) + Y^+_l(t)A^T(t) -\nonumber \\ & \pi_l(t)Y^+_l(t) - \frac{1}{\pi_l(t)}B(t)P(t)B^T(t) +\nonumber \\ & (Y_l^{+}(t))^{1/2}S_l(t)(G(t)Q(t)G^T(t))^{1/2} +\nonumber \\ & (G(t)Q(t)G^T(t))^{1/2}S_l^T(t)(Y_l^{+}(t))^{1/2},\\

(73)Y^+_l(t_1) = Y_1,

in which

\pi_l(t) = \frac{\langle l(t), B(t)P(t)B^T(t)l(t)\rangle^{1/2}}{\langle l(t), Y^+_l(t)l(t)\rangle^{1/2}},

and the orthogonal matrix S_l(t) satisfies the equation

S_l(t)(G(t)Q(t)G^T(t))^{1/2}l(t) = \frac{\langle l(t), G(t)Q(t)G^T(t)l(t)\rangle^{1/2}}{\langle l(t), Y_l^+(t)l(t)\rangle^{1/2}}(Y_l^{+}(t))^{1/2}l(t).

The equation for the shape matrix of the internal ellipsoid is

(74)\dot{Y}^-_l(t) & = A(t)Y^-_l(t) + Y^-_l(t)A^T(t) -\nonumber \\ & (Y_l^{-}(t))^{1/2}T_l(t)(B(t)P(t)B^T(t))^{1/2} -\nonumber \\ & (B(t)P(t)B^T(t))^{1/2}T_l^T(t)(Y_l^{-}(t))^{1/2} +\nonumber \\ & \eta_l(t)Y^-_l(t) + \frac{1}{\eta_l(t)}G(t)Q(t)G^T(t),\\

(75)Y^-_l(t_1) & = Y_1,

in which

\eta_l(t) = \frac{\langle l(t), G(t)Q(t)G^T(t)l(t)\rangle^{1/2}}{\langle l(t), Y^+_l(t)l(t)\rangle^{1/2}},

and the orthogonal matrix T_l(t) is determined by the equation

T_l(t)(B(t)P(t)B^T(t))^{1/2}l(t) = \frac{\langle l(t), B(t)P(t)B^T(t)l(t)\rangle^{1/2}}{\langle l(t), Y_l^-(t)l(t)\rangle^{1/2}}(Y_l^{-}(t))^{1/2}l(t).

Just as in the forward reachability case, the terms containing G(t) and Q(t) vanish from equations (72) and (74) in the absence of disturbances. The boundary value problems (68), (72) and (74) are converted to the initial value problems by the change of variables s = -t.

Due to (70) and (71),

\bigcup_{\langle l_1,l_1\rangle=1}{\mathcal E}(y_c(t),Y^-_l(t)) = {\mathcal Y}(t_1,t,{\mathcal E}(y_1,Y_1)) = \bigcap_{\langle l_1,l_1\rangle=1}{\mathcal E}(y_c(t),Y^+_l(t)).

Remark. In expressions (62), (64), (72) and (74) the terms \frac{1}{\pi_l(t)} and \frac{1}{\eta_l(t)} may not be well defined for some vectors l, because matrices B(t)P(t)B^T(t) and G(t)Q(t)G^T(t) may be singular. In such cases, we set these entire expressions to zero.

Discrete-time systems

Consider the discrete-time linear system,

(76)x(t+1) = A(t)x(t) + B(t)u(t,x(t)) + G(t)v(t),

in which x(t)\in{\bf R}^n is the state, u(t, x(t))\in{\bf R}^m is the control bounded by the ellipsoid {\mathcal E}(p(t),P(t)), v(t)\in{\bf R}^d is disturbance bounded by ellipsoid {\mathcal E}(q(t),Q(t)), and matrices A(t), B(t), G(t) are in {\bf R}^{n\times n}, {\bf R}^{n\times m}, {\bf R}^{n\times d} respectively. Here we shall assume A(t) to be nonsingular. [9] The set of initial conditions at initial time t_0 is ellipsoid {\mathcal E}(x_0,X_0).

Ellipsoidal Toolbox computes maxmin and minmax CLRS \overline{{\mathcal X}}_{CL}(t, t_0, {\mathcal E}(x_0, X_0) and \underline{{\mathcal X}}_{CL}(t, t_0, {\mathcal E}(x_0, X_0) for discrete-time systems.

If matrix Q(\cdot)=0, the system (76) becomes an ordinary affine system with known v(\cdot)=q(\cdot). If matrix G(\cdot)=0, the system reduces to a linear controlled system. In the absence of disturbance (Q(\cdot)=0 or G(\cdot)=0), \overline{{\mathcal X}}_{CL}(t,t_0,{\mathcal E}(x_0,X_0))=\underline{{\mathcal X}}_{CL}(t,t_0,{\mathcal E}(x_0,X_0))={\mathcal X}(t,t_0,{\mathcal E}(x_0,X_0)), the reach set is as in definition (3).

Maxmin and minmax CLRS \overline{{\mathcal X}}_{CL}(t, t_0, {\mathcal E}(x_0, X_0) and \underline{{\mathcal X}}_{CL}(t, t_0, {\mathcal E}(x_0, X_0), if nonempty, are symmetric convex and compact, with the center evolving in time according to

(77)x_c(t+1) = A(t)x_c(t) + B(t)p(t) + G(t)v(t), \\ x_c(t_0)= x_0.

Fix some vector l_0\in{\bf R}^n and consider l(t) that satisfies the discrete-time adjoint equation, [10]

(78)l(t+1) = \left(A^T\right)^{-1}(t)l(t), \\ l(t_0) = l_0,

or, equivalently

l(t) = \Phi^T(t_0, t)l_0 .

There exist tight external ellipsoids {\mathcal E}(x_c(t), \overline{X}^+_l(t)), {\mathcal E}(x_c(t), \underline{X}^+_l(t)) and tight internal ellipsoids {\mathcal E}(x_c(t), \overline{X}^-_l(t)), {\mathcal E}(x_c(t), \underline{X}^-_l(t)) such that

(79){\mathcal E}(x_c(t), \overline{X}^-_l(t))\subseteq\overline{{\mathcal X}}_{CL}(t,t_0,{\mathcal E}(x_0,X_0)) \subseteq {\mathcal E}(x_c(t), \overline{X}^+_l(t)),

(80)\rho(l(t) ~|~ {\mathcal E}(x_c(t), \overline{X}^-_l(t))) = \rho(l(t) ~|~ \overline{{\mathcal X}}_{CL}(t, t_0, {\mathcal E}(x_0,X_0))) = \rho(l(t) ~|~ {\mathcal E}(x_c(t), \overline{X}^+_l(t))) .

and

(81){\mathcal E}(x_c(t), \underline{X}^-_l(t))\subseteq\underline{{\mathcal X}}_{CL}(t,t_0,{\mathcal E}(x_0,X_0)) \subseteq {\mathcal E}(x_c(t), \underline{X}^+_l(t)),

(82)\rho(l(t) ~|~ {\mathcal E}(x_c(t), \underline{X}^-_l(t))) = \rho(l(t) ~|~ \underline{{\mathcal X}}_{CL}(t, t_0, {\mathcal E}(x_0,X_0))) = \rho(l(t) ~|~ {\mathcal E}(x_c(t), \underline{X}^+_l(t))) .

The shape matrix of the external ellipsoid for maxmin reach set is determined from

(83)\hat{X}^+_l(t) & = (1+\overline{\pi}_l(t))A(t)\overline{X}^+_l(t)A^T(t) + \left(1+\frac{1}{\overline{\pi}_l(t)}\right) B(t)P(t)B^T(t), \\

(84)\overline{X}^+_l(t+1) & = \left((\hat{X}^+_l(t))^{1/2} + \overline{S}_l(t)(G(t)Q(t)G^T(t))^{1/2}\right)^T \times \nonumber \\ & \left((\hat{X}^+_l(t))^{1/2} + \overline{S}_l(t)(G(t)Q(t)G^T(t))^{1/2}\right),\\

(85)\overline{X}^+_l(t_0) & = X_0,

wherein

\overline{\pi}_l(t) = \frac{\langle l(t+1), B(t)P(t)B^T(t)l(t+1)\rangle^{1/2}}{\langle l(t), \overline{X}^+_l(t)l(t)\rangle^{1/2}},

and the orthogonal matrix \overline{S}_l(t) is determined by the equation

\begin{aligned} & & \overline{S}_l(t)(G(t)Q(t)G^T(t))^{1/2}l(t+1) = \\ & & \frac{\langle l(t+1), G(t)Q(t)G^T(t)l(t+1)\rangle^{1/2}}{\langle l(t+1), \hat{X}^+_l(t)l(t+1)\rangle^{1/2}}(\hat{X}^+_l(t))^{1/2}l(t+1).\end{aligned}

Equation (84) is valid only if {\mathcal E}(0,G(t)Q(t)G^T(t))\subseteq{\mathcal E}(0,\hat{X}^+_l(t)), otherwise the maxmin CLRS \overline{{\mathcal X}}_{CL}(t,t_0,{\mathcal E}(x_0,X_0)) is empty.

The shape matrix of the external ellipsoid for minmax reach set is determined from

(86)\breve{X}^+_l(t) & = \left((A(t)\underline{X}^+_l(t)A^T(t))^{1/2} + \underline{S}_l(t)(G(t)Q(t)G^T(t))^{1/2}\right)^T \times \nonumber \\ &\left((A(t)\underline{X}^+_l(t)A^T(t))^{1/2} + \underline{S}_l(t)(G(t)Q(t)G^T(t))^{1/2}\right)\\

(87)\underline{X}^+_l(t+1) & = (1+\underline{\pi}_l(t))\breve{X}^+_l(t) + \left(1+\frac{1}{\underline{\pi}_l(t)}\right) B(t)P(t)B^T(t),\\

(88)\underline{X}^+_l(t_0) & = X_0,

where

\underline{\pi}_l(t) = \frac{\langle l(t+1), B(t)P(t)B^T(t)l(t+1)\rangle^{1/2}}{\langle l(t+1), \breve{X}^+_l(t)l(t+1)\rangle^{1/2}},

and \underline{S}_l(t) is orthogonal matrix determined from the equation

\begin{aligned} & \underline{S}_l(t)(G(t)Q(t)G^T(t))^{1/2}l(t+1) = \\ & \frac{\langle l(t+1), G(t)Q(t)G^T(t)l(t+1)\rangle^{1/2}}{\langle l(t), \underline{X}^+_l(t)l(t)\rangle^{1/2}}(A(t)\underline{X}^+_l(t)A^T(t))^{1/2}l(t+1).\end{aligned}

Equations (86), (87) are valid only if {\mathcal E}(0,G(t)Q(t)G^T(t)\subseteq{\mathcal E}(0,A(t)\underline{X}^+_l(t)A^T(t)), otherwise minmax CLRS \underline{{\mathcal X}}_{CL}(t,t_0,{\mathcal E}(x_0,X_0)) is empty.

The shape matrix of the internal ellipsoid for maxmin reach set is determined from

(89)\hat{X}^-_l(t) & = \left((A(t)\overline{X}^-_l(t)A^T(t))^{1/2} + \overline{T}_l(t)(B(t)P(t)B^T(t))^{1/2}\right)^T \times \nonumber \\ & \left((A(t)\overline{X}^-_l(t)A^T(t))^{1/2} + \overline{T}_l(t)(B(t)P(t)B^T(t))^{1/2}\right)\\

(90)\overline{X}^-_l(t+1) & = (1+\overline{\eta}_l(t))\hat{X}^-_l(t) + \left(1+\frac{1}{\underline{\eta}_l(t)}\right) G(t)Q(t)G^T(t), \\

(91)\overline{X}^-_l(t_0) & = X_0,

where

\overline{\eta}_l(t) = \frac{\langle l(t+1), G(t)Q(t)G^T(t)l(t+1)\rangle^{1/2}}{\langle l(t+1), \hat{X}^-_l(t)l(t+1)\rangle^{1/2}},

and \overline{T}_l(t) is orthogonal matrix determined from the equation

\begin{aligned} & \overline{T}_l(t)(B(t)P(t)B^T(t))^{1/2}l(t+1) = \\ & \frac{\langle l(t+1), B(t)P(t)B^T(t)l(t+1)\rangle^{1/2}}{\langle l(t), \overline{X}^-_l(t)l(t)\rangle^{1/2}}(A(t)\overline{X}^-_l(t)A^T(t))^{1/2}l(t+1).\end{aligned}

Equation (90) is valid only if {\mathcal E}(0,G(t)Q(t)G^T(t)\subseteq{\mathcal E}(0,\hat{X}^-_l(t)).

The shape matrix of the internal ellipsoid for the minmax reach set is determined by

(92)\breve{X}^-_l(t) & = (1+\underline{\eta}_l(t))A(t)\underline{X}^-_l(t)A^T(t) + \left(1+\frac{1}{\underline{\eta}_l(t)}\right) G(t)Q(t)G^T(t),\\

(93)\underline{X}^-_l(t+1) & = \left((\breve{X}^-_l(t))^{1/2} + \underline{T}_l(t)(B(t)P(t)B^T(t))^{1/2}\right)^T \times \nonumber \\ &\left((\breve{X}^-_l(t))^{1/2} + \underline{T}_l(t)(B(t)P(t)B^T(t))^{1/2}\right),\\

(94)\underline{X}^-_l(t_0) & = X_0,

wherein

\underline{\eta}_l(t) = \frac{\langle l(t+1), G(t)Q(t)G^T(t)l(t+1)\rangle^{1/2}}{\langle l(t), \underline{X}^-_l(t)l(t)\rangle^{1/2}},

and the orthogonal matrix \underline{T}_l(t) is determined by the equation

\begin{aligned} &\underline{T}_l(t)(B(t)P(t)B^T(t))^{1/2}l(t+1) = \\ & \frac{\langle l(t+1), B(t)P(t)B^T(t)l(t+1)\rangle^{1/2}}{\langle l(t+1), \breve{X}^-_l(t)l(t+1)\rangle^{1/2}}(\breve{X}^-_l(t))^{1/2}l(t+1).\end{aligned}

Equations (92), (93) are valid only if {\mathcal E}(0,G(t)Q(t)G^T(t)\subseteq{\mathcal E}(0,A(t)\underline{X}^-_l(t)A^T(t)).

The point where the external and the internal ellipsoids both touch the boundary of the maxmin CLRS is

x_l^+(t) = x_c(t) + \frac{\overline{X}^+_l(t)l(t)}{\langle l(t), \overline{X}^+_l(t)l(t)\rangle^{1/2}} ,

and the bounday point of minmax CLRS is

x_l^-(t) = x_c(t) + \frac{\overline{X}^-_l(t)l(t)}{\langle l(t), \overline{X}^-_l(t)l(t)\rangle^{1/2}} .

Points x^{\pm}_l(t), t\geqslant t_0, form extremal trajectories. In order for the system to follow the extremal trajectory specified by some vector l_0, the initial state must be

(95)x_l^0 = x_0 + \frac{X_0l_0}{\langle l_0, X_0l_0\rangle^{1/2}}.

When there is no disturbance (G(t)=0 or Q(t)=0), \overline{X}^+_l(t)=\underline{X}^+_l(t) and \overline{X}^-_l(t)=\underline{X}^-_l(t), and the open-loop control that steers the system along the extremal trajectory defined by l_0 is

(96)u_l(t) = p(t) + \frac{P(t)B^T(t)l(t+1)}{\langle l(t+1), B(t)P(t)B^T(t)l(t+1)\rangle^{1/2}}.

Each choice of l_0 defines an external and internal approximation. If \overline{{\mathcal X}}_{CL}(t,t_0,{\mathcal E}(x_0,X_0)) is nonempty,

\bigcup_{\langle l_0,l_0\rangle=1}{\mathcal E}(x_c(t),\overline{X}^-_l(t)) = \overline{{\mathcal X}}_{CL}(t,t_0,{\mathcal E}(x_0,X_0)) = \bigcap_{\langle l_0,l_0\rangle=1}{\mathcal E}(x_c(t),\overline{X}^+_l(t)).

Similarly for \underline{{\mathcal X}}_{CL}(t,t_0,{\mathcal E}(x_0,X_0)),

\bigcup_{\langle l_0,l_0\rangle=1}{\mathcal E}(x_c(t),\underline{X}^-_l(t)) = \underline{{\mathcal X}}_{CL}(t,t_0,{\mathcal E}(x_0,X_0)) = \bigcap_{\langle l_0,l_0\rangle=1}{\mathcal E}(x_c(t),\underline{X}^+_l(t)).

Similarly, tight ellipsoidal approximations of maxmin and minmax CLBRS with terminating conditions (t_1, {\mathcal E}(y_1,Y_1)) can be obtained for those directions l(t) satisfying

(97)l(t) = \Phi^T(t_1,t)l_1,

with some fixed l_1, for which they exist.

With boundary conditions

(98)y_c(t_1)=y_1, ~~~ \overline{Y}^+_l(t_1)=\overline{Y}^-_l(t_1)=\underline{Y}^+_l(t_1)=\underline{Y}^-_l(t_1)=Y_1,

external and internal ellipsoids for maxmin CLBRS \overline{{\mathcal Y}}_{CL}(t_1,t,{\mathcal E}(y_1,Y_1)) at time t, {\mathcal E}(y_c(t),\overline{Y}^+_l(t)) and {\mathcal E}(y_c(t),\overline{Y}^-_l(t)), are computed as external and internal ellipsoidal approximations of the geometric sum-difference

A^{-1}(t)\left( {\mathcal E}(y_c(t+1),\overline{Y}^+_l(t+1)) \oplus B(t){\mathcal E}(-p(t),P(t)) \dot{-}G(t){\mathcal E}(-q(t),Q(t)) \right)

and

A^{-1}(t)\left( {\mathcal E}(y_c(t+1),\overline{Y}^-_l(t+1)) \oplus B(t){\mathcal E}(-p(t),P(t)) \dot{-}G(t){\mathcal E}(-q(t),Q(t)) \right)

in direction l(t) from (97). Section Geometric Sum-Difference describes the operation of geometric sum-difference for ellipsoids.

External and internal ellipsoids for minmax CLBRS \underline{{\mathcal Y}}_{CL}(t_1,t,{\mathcal E}(y_1,Y_1)) at time t, {\mathcal E}(y_c(t),\underline{Y}^+_l(t)) and {\mathcal E}(y_c(t),\underline{Y}^-_l(t)), are computed as external and internal ellipsoidal approximations of the geometric difference-sum

A^{-1}(t)\left( {\mathcal E}(y_c(t+1),\underline{Y}^+_l(t+1)) \dot{-}G(t){\mathcal E}(-q(t),Q(t)) \oplus B(t){\mathcal E}(-p(t),P(t)) \right)

and

A^{-1}(t)\left( {\mathcal E}(y_c(t+1),\underline{Y}^-_l(t+1)) \dot{-}G(t){\mathcal E}(-q(t),Q(t)) \oplus B(t){\mathcal E}(-p(t),P(t)) \right)

in direction l(t) from (97). Section Geometric Difference-Sum describes the operation of geometric difference-sum for ellipsoids.

[1]In discrete-time case t assumes integer values.
[2]We are being general when giving the basic definitions. However, it is important to understand that for any specific continuous-time dynamical system it must be determined whether the solution exists and is unique, and in which class of solutions these conditions are met. Here we shall assume that function f is such that the solution of the differential equation (1) exists and is unique in Fillipov sense. This allows the right-hand side to be discontinuous. For discrete-time systems this problem does not exist.
[3]Minkowski sum of sets {\mathcal W}, {\mathcal Z}\subseteq {\bf R}^n is defined as {\mathcal W}\oplus {\mathcal Z}= \{w+z ~|~ w\in{\mathcal W}, ~ z\in{\mathcal Z}\}. Set {\mathcal W}\oplus{\mathcal Z} is nonempty if and only if both, {\mathcal W} and {\mathcal Z} are nonempty. If {\mathcal W} and {\mathcal Z} are convex, set {\mathcal W}\oplus{\mathcal Z} is convex.
[4]{\mathcal M} is weakly invariant with respect to the target set {\mathcal Y}_1 and times t_0 and t, if for every state x_0\in{\mathcal M} there exists a control u(\tau, x(\tau))\in{\mathcal U}(\tau), t_0\leqslant\tau< t, that steers the system from x_0 at time t_0 to some state in {\mathcal Y}_1 at time t. If all controls in {\mathcal U}(\tau), t_0\leqslant\tau<t steer the system from every x_0\in{\mathcal M} at time t_0 to {\mathcal Y}_1 at time t, set {\mathcal M} is said to be strongly invariant with respect to {\mathcal Y}_1, t_0 and t.
[5]There exists f^{-1}(t,x,u) such that x(t)=f^{-1}(t, x(t+1), u, v).
[6]

Hausdorff semidistance between compact sets {\mathcal W}, {\mathcal Z}\subseteq {\bf R}^n is defined as

{\bf dist}({\mathcal W}, {\mathcal Z}) = \min\{\langle w-z, w-z\rangle^{1/2} ~|~ w\in{\mathcal W}, \; z\in{\mathcal Z}\},

where \langle\cdot, \cdot\rangle denotes inner product.
[7]The Minkowski difference of sets {\mathcal W}, {\mathcal Z}\in{\bf R}^n is defined as {\mathcal W}\dot{-}{\mathcal Z}= \left\{\xi\in{\bf R}^n ~|~ \xi \oplus {\mathcal Z}\subseteq {\mathcal W}\right\}. If {\mathcal W} and {\mathcal Z} are convex, {\mathcal W}\dot{-}{\mathcal Z} is convex if it is nonempty.
[8]So-called verification problems often consist in ensuring that the system is unable to reach an ‘unsafe’ target set within a given time interval.
[9]

The case when A(t) is singular is described in [VAR2007]. The idea is to substitute A(t) with the nonsingular A_\delta(t) = A(t) + \delta U(t)W(t), in which U(t) and W(t) are obtained from the singular value decomposition

A(t) = U(t)\Sigma(t)V(t) .

The parameter \delta can be chosen based on the number of time steps for which the reach set must be computed and the required accuracy. The issue of inverting ill-conditioned matrices is also addressed in [VAR2007].
[10]Note that for (78) A(t) must be invertible.

References

[VAR2007](1, 2) P. Varaiya A. A. Kurzhanskiy. Ellipsoidal techniques for reachability analysis of discrete-time linear systems. IEEE Transactions on Automatic Control, 52(1):26–38, 2007. linear systems. IEEE Transactions on Automatic Control, 52(1):26–38, 2007.

Table Of Contents

Previous topic

Ellipsoidal Calculus

Next topic

Installation

This Page

Navigation

© Copyright 2011-2014 Moscow State University, Faculty of Computational Mathematics and Computer Science, System Analysis Department, 2004-2011 The Regents of the University of California. Created using Sphinx 1.2.2.