text:
# Approximate Computation of Reach Sets in Hybrid Systems

## 1 Introduction

Hybrid systems combine discrete and continuous dynamics. The dynamics of the continuous variables within a discrete state are specified by differential equations or differential inclusions. An important problem in the analysis and sysnthesis of such systems is the so called reachability problem, which asks, for two sets of configurations of a given hybrid system, say $`𝒳_1`$ and $`𝒳_2`$, where a configuration consists of discrete and continuous components, whether or not there is a hybrid trajectory with the initial configuration in $`𝒳_1`$ and the final configuration in $`𝒳_2`$. A hybrid trajectory may be described as a trajectory of configurations consisting of discrete state jumps and smooth arcs, where each arc evolves according to the continuous dynamics of a discrete state, with the starting and end points of each arc satisfying the jump conditions of discrete transitions. A more precise definition is presented in Sec. 2, which provides a concise introduction to hybrid systems and the reachability problem.

The reachability problem is undecidable for certain classes of linear hybrid systems (i.e., hybrid systems having linear trajectories within each discrete state and linear or constant reset maps, also called constant slope hybrid systems) , although in some cases decidability results have been obtained  (see also ). Therefore, for a general hybrid system, a reasonable alternative appears to be to find semi-decision procedures for the reachability problem . A computational approach to this problem also requires finding the set of states reached by the continuous variables evolving according to the dynamics of a discrete state. In this paper, we consider the problem of computing and representing the reach sets of the continuous variables within a discrete state when the dynamics of the continuous variables are specified by differential equations with initial conditions belonging to a specified initial set.

In this context, various methods have been proposed in the literature for finding reach sets of continuous variables . In Sec. 3, we describe a method for computing the reach sets based on the idea that a subset of the boundary of the initial set may be found, such that it is sufficient to compute the solutions with the initial points lying in that set. This method is similar to that in  (and also to some extent to that described in ). In particular, we present a schematic algorithm, which is somewhat more general in its scope than that in  and simpler compared to that in .

Besides these algorithms for reach set computation, an equally important issue the representation of the reach sets for manipulating the sets efficiently. This requires representation of the reach sets in terms of more convenient sets such as, for example polyhedral or subalgebraic sets, that are simple to represent and easy to handle for practical purposes. However, since the representing class of sets may not contain a member that exactly equals the reach set, we may have to find approximations to the reach sets by those that belong to the representing class. (See, for example,  for approximation with polyhedra, and  for approximation using ellipsoids). Typically, over-approximations may be used for verifying whether a safety requirement may potentially be violated by any of the behaviours starting from a given initial set, while under-approximations are needed for characterizing a set of states from which a desirable property is always achievable. We describe in Sec. 4 a method for over-approximating the reach sets by polyhedra when the dynamics are specified by linear differential equations and the initial set is a polyhedron. These algorithms are illustrated with some simple examples in Sec. 5, while Sec. 6 concludes the paper.

## 2 A Brief Survey of Hybrid Systems

In this section, we present a brief introduction to hybrid systems, and provide motivation for the remaining sections.

### 2.1 Preliminary Definitions

We begin with a somewhat detailed but general definition of hybrid systems.

###### Definition 1

A hybrid system is a tuple $`H=(Q,X,\mathrm{\Sigma },G,E,Init,f)`$, where

1. $`Q`$ is a finite set of discrete states (also called locations).
2. $`X=^n`$, $`n1`$, is the set of continuous states, where $``$ denotes the set of real numbers. We demote the continuous state variable by $`𝐱`$.
3. $`\mathrm{\Sigma }`$ is a finite set of discrete events or environment actions. Some events in $`\mathrm{\Sigma }`$ are controllable, while others not. Hence, it is convenient to assume that $`\mathrm{\Sigma }=\mathrm{\Sigma }_c\mathrm{\Sigma }_d`$, where $`\mathrm{\Sigma }_c`$ is the set of controllable events, and $`\mathrm{\Sigma }_d`$ is the set of uncontrollable or disturbance events.
4. $`GQ\times X`$ is a set of state invariance conditions. When the system is in state $`q`$, the continuous variables belong to $`G(q)=\{𝐱X:(q,𝐱)G\}`$.
5. $`EQ\times 𝒫(^n)\times \mathrm{\Sigma }\times \{^n^n\}\times Q`$ is the set of transition edges. An edge $`eE`$, where $`e=(q_e,X_e,\sigma _e,r_e,q_e^{})`$, is interpreted as follows:

* If the continuous state is in $`X_e`$ and the event $`\sigma _e`$ occurs, then transition edge $`e`$ is enabled in state $`q_e`$. Thus $`X_e`$ is the set of switching points of the continuous variables from state $`q_e`$ to $`q_e^{}`$. The set $`X_e`$ is specified by a predicate, and is also called a guard. In case there are many transition edges simultaneously enabled, then the system may select one of the edges nondeterministically.
    * If a transition edge $`e`$ is selected by the system, then the continuous state is reset using the function $`r_e`$ when the system enters state $`q_e^{}`$. The reset values obey the state invariance condition, $`r_e(X_e)G(q_e^{})`$.
6. $`InitG`$ is a set of initial conditions.
7. $`f`$ is an $`n`$-dimensional vector field with real-valued components governing the dynamics of the continuous state $`𝐱`$. The domain of definition of the function $`f`$ is a set $`𝒟`$, where $`𝒟=Q\times X`$ (if there are no continuous control variables) or $`𝒟=Q\times X\times U`$ (if $`U^m`$, $`m1`$, is the range of the continuous control variable, which is a function $`𝐮:[0,T_q)U`$, where $`T_q>0`$ is as large as may be necessary, depending on the discrete state $`qQ`$). When in state $`q`$, the continuous variables evolve according to the dynamical law

$$\frac{d𝐱}{dt}=f(q,𝐱)$$

or according to the dynamical law

$$\frac{d𝐱}{dt}=f(q,𝐱,𝐮)$$

depending on the presence or absence of the continuous control variables.111Another possibility is to specify the system dynamics of the continuous variables in terms of differential inclusions. In this case, $`f(q,𝐱)`$ is a set, and the continuous variables evolve according to
$$\frac{d𝐱}{dt}f(q,𝐱).$$
But we do not consider this case here.

The initial conditions of the continuous evolution are specified either by the reset maps when the system takes a discrete transition and reaches the state $`q`$, or by the (nondeterministic) initial conditions of the system given as the set $`Init`$.

If, in addition to the set of initial values, $`Init`$, the set of final or accepting values in $`Q\times X`$, denoted by $`Final`$, are also specified, then the hybrid system is called a hybrid automaton.

Referring to the system dynamics in the above definition, for our purposes, we assume that there are no continuous control variables. Therefore, in the sequel, we assume that the dynamics of the continuous variables are specified in the form

$$\frac{d𝐱}{dt}=f(q,𝐱)$$

without any continuous control variables.

###### Definition 2

Let $`H=(Q,X,\mathrm{\Sigma },G,E,Init,f)`$ be a hybrid system. Each point $`(q,𝐱)Q\times X`$ is called a configuration of the hybrid system $`H`$.

In Definitions 3 \- 6, we introduce certain terminology, starting with the notion of a step of a hybrid system, leading upto the notions of predecessors and successors of a configuration $`(q,𝐱)`$. But before that, we mention that predecessors and successors may be defined irrespective of initial conditions, whereas for defining the executions of a hybrid system, we require the specification of the initial set $`Init`$, as will be seen later (refer to Definitions 4, 5 and 6) .

###### Definition 3

Let $`H=(Q,X,\mathrm{\Sigma },G,E,Init,f)`$ be a hybrid system, and let $`(q,𝐱)`$ and $`(q^{},𝐱^{})`$ be two configurations of $`H`$. Then, the pair of configurations $`((q,𝐱),(q^{},𝐱^{}))`$ is called

1. a time-step, if $`q=q^{}`$ and, for some $`t0`$, there is a function $`𝐲:[0,t]X`$, satisfying $`\dot{𝐲}(s)=f(q,𝐲(s))`$, $`𝐲(s)G(q)`$, for $`s[0,t]`$, $`𝐲(0)=𝐱`$ and $`𝐲(t)=𝐱^{}`$.
2. an edge-step if there is an edge $`e=(q_e,X_e,\sigma _e,r_e,q_e^{})E`$, such that $`q=q_e`$, $`q^{}=q_e^{}`$, $`𝐱X_e`$ and $`𝐱^{}=r_e(𝐱)`$.
3. a $`\sigma `$-step, where $`\sigma \mathrm{\Sigma }`$, if there is an edge $`\mathrm{}=(q_{\mathrm{}},X_{\mathrm{}},\sigma _{\mathrm{}},r_{\mathrm{}},q_{\mathrm{}}^{})E`$ with $`\sigma =\sigma _{\mathrm{}}`$, such that $`q=q_{\mathrm{}}`$, $`q^{}=q_{\mathrm{}}^{}`$, $`𝐱X_{\mathrm{}}`$ and $`𝐱^{}=r_{\mathrm{}}(𝐱)`$.

A step of the hybrid system $`H`$ is a pair of configurations $`((q,𝐱),(q^{},𝐱^{}))`$, such that $`((q,𝐱),(q^{},𝐱^{}))`$ is either a time-step for some $`t0`$, or an edge-step for some $`eE`$, or a $`\sigma `$-step for some $`\sigma \mathrm{\Sigma }`$.

Of course, every edge-step is a $`\sigma `$-step for some $`\sigma \mathrm{\Sigma }`$, and conversely, every $`\sigma `$-step is an edge-step for some edge $`eE`$. Hence these definitions may seem redundant. However, the distinction between the two types of transitions would become obvious if there is a discrete state controller, defined as a function $`C:Q\times X\mathrm{\Sigma }_c`$, that triggers discrete controllable events to enable a particular transition depending on the present configuration (see for instance  for the case of timed automata). But in this work, we shall not have occasion to discuss about discrete state controllers.

###### Definition 4

(Trajectories and Executions of a Hybrid System) A hybrid trajectory or simply a trajectory of $`H`$ is a sequence of configurations $`(q_1,𝐱_1),(q_2,𝐱_2),(q_3,𝐱_3),\mathrm{}`$, where for each $`i1`$, $`((q_i,𝐱_i),(q_{i+1},𝐱_{i+1}))`$ is a step. A trajectory is

1. finite, if the number of steps is finite;
2. an execution, if $`(q_1,𝐱_1)Init`$, where $`(q_1,𝐱_1)`$ is the initial configuration; and
3. a finite execution, if it is both.

We now define, for each $`\sigma \mathrm{\Sigma }`$, the set-valued successor function, $`Post_\sigma :Q\times X2^Q\times 2^X`$. As will be briefly mentioned later (after Definition 6), our definitions lead in a natural way to extract a labeled transition system (see ) from a hybrid system.

###### Definition 5

($`\sigma `$-Successors of a Configuration)
Let $`H=(Q,X,\mathrm{\Sigma },G,E,Init,f)`$ be a hybrid system, $`\sigma \mathrm{\Sigma }`$ and $`(q,𝐱)Q\times X`$. Then, we define $`Post_\sigma (q,𝐱)`$ to be the union of the two sets $`S_1,S_2^\sigma Q\times X`$, where

1. $`S_1=\{(q^{},𝐱^{}):((q,𝐱),(q^{},𝐱^{}))\text{ is a time-step for some }t0\}`$; and
2. $`S_2^\sigma =\{(q^{\prime \prime },𝐱^{\prime \prime }):((q^{},𝐱^{}),(q^{\prime \prime },𝐱^{\prime \prime }))\text{ is a }\sigma \text{-step for some }(q^{},𝐱^{})S_1\text{ }\}`$.

Further, if $`(q_s,𝐱_s)Post_\sigma (q,𝐱)`$, then $`(q_s,𝐱_s)`$ is called a $`\sigma `$-successor of $`(q,𝐱)`$.

The set $`S_1`$, as in the first part of this definition, is of main interest to us in the later sections of this paper. Specifically, we will deal with the set $`Reach_q^G\left(X_0\right)`$, where $`qQ`$ and $`X_0G(q)`$, defined as

$$Reach_q^G\left(X_0\right)=\{𝐱:((q,𝐱_0),(q,𝐱))\text{ is a time-step for some}\text{ }t0\text{ and }𝐱_0X_0\text{ }\}.$$
(1)
In the subsequent sections, we will be concerned more with this and a related set. But for now, we shall proceed with our discussion with the following definition:

###### Definition 6

($`\sigma `$-Predecessors of a Configuration)
Let $`H=(Q,X,\mathrm{\Sigma },G,E,Init,f)`$ be a hybrid system, $`\sigma \mathrm{\Sigma }`$ and $`(q,𝐱)Q\times X`$. Then, we define $`Pre_\sigma (q,𝐱)`$ to be the union of the two sets $`P_1,P_2^\sigma Q\times X`$, where

1. $`P_1=\{(q^{},𝐱^{}):((q^{},𝐱^{}),(q,𝐱))\text{ is a time-step for some }t0\}`$; and
2. $`P_2^\sigma =\{(q^{\prime \prime },𝐱^{\prime \prime }):((q^{\prime \prime },𝐱^{\prime \prime }),(q^{},𝐱^{}))\text{ is a }\sigma \text{-step for some }(q^{},𝐱^{})P_1\text{ }\}`$.

Further, if $`(q_p,𝐱_p)Pre_\sigma (q,𝐱)`$, then $`(q_p,𝐱_p)`$ is called a $`\sigma `$-predecessor of $`(q,𝐱)`$.

Note that the notion of $`\sigma `$-successor generalizes the notion of $`\sigma `$-step by including time-steps. The set valued function $`Post_\sigma (q,𝐱)`$ defines, in a natural way, a transition relation, $`\stackrel{\sigma }{}(Q\times X)\times (Q\times X)`$, as follows: $`(q,𝐱)\stackrel{\sigma }{}(q^{},𝐱^{})`$, if $`(q^{},𝐱^{})Post_\sigma (q,𝐱)`$. The transition relation $`\stackrel{\sigma }{}`$ is also called $`\sigma `$-transition relation. Similarly, the set valued function $`Pre_\sigma (q,𝐱)`$ defines another transition relation, $`\stackrel{\sigma }{}(Q\times X)\times (Q\times X)`$, as follows: $`(q,𝐱)\stackrel{\sigma }{}(q^{},𝐱^{})`$, if $`(q,𝐱)Pre_\sigma (q^{},𝐱^{})`$.

For two configurations $`(q_1,𝐱_1)`$ and $`(q_2,𝐱_2)`$, if $`((q_1,𝐱_1),(q_2,𝐱_2))`$ is a $`\sigma `$-step, then $`(q_2,𝐱_2)Post(q_1,𝐱_1)`$ and $`(q_1,𝐱_1)Pre(q_2,𝐱_2)`$; hence, in this case, we have $`(q_1,𝐱_1)\stackrel{\sigma }{}(q_2,𝐱_2)`$ and $`(q_1,𝐱_1)\stackrel{\sigma }{}(q_2,𝐱_2)`$. This may mislead us to get the false impression that $`\stackrel{\sigma }{}`$ and $`\stackrel{\sigma }{}`$ are the same; and to avoid any such possible confusion, we emphasize that, in general, it is not true that $`(q_1,𝐱_1)\stackrel{\sigma }{}(q_2,𝐱_2)`$ implies $`(q_1,𝐱_1)\stackrel{\sigma }{}(q_2,𝐱_2)`$; and the same with the converse statement. Hence, the two relations $`\stackrel{\sigma }{}`$ and $`\stackrel{\sigma }{}`$ are not the same.

The definitions of $`Post_\sigma `$ and $`Pre_\sigma `$ can be extended straightforwardly to subsets of $`Q\times X`$, as follows: for $`SQ\times X`$, define

$$Post_\sigma (S)=\underset{(q,𝐱)S}{}Post_\sigma (q,𝐱),\text{ and }Pre_\sigma (S)=\underset{(q,𝐱)S}{}Pre_\sigma (q,𝐱).$$

Finally, define

$$Post(S)=\underset{\sigma \mathrm{\Sigma }}{}Post_\sigma (S),\text{ and }Pre(S)=\underset{\sigma \mathrm{\Sigma }}{}Pre_\sigma (S).$$

We sometimes refer to $`Post`$ and $`Pre`$ as the 1-step transition functions222This should not be confused with the notion of step as defined in Definition 3.. More generally, the $`k`$-step transition functions $`Post^k`$ and $`Pre^k`$ are defined inductively as follows: For $`SQ\times X`$,

$$Post^1(S)=Post(S),\text{ and for }k2\text{}Post^k(S)=Post(Post^{k1}(S)).$$

Similarly,

$$Pre^1(S)=Pre(S),\text{ and for }k2\text{}Pre^k(S)=Pre(Pre^{k1}(S)).$$

Finally, define $`Post^{}`$ and $`Pre^{}`$ as

$$Post^{}(S)=\underset{k1}{}Post^k(S),\text{ and }Pre^{}(S)=\underset{k1}{}Pre^k(S).$$

### 2.2  Reachability Problem for Hybrid Systems

In the notation just discussed, the reachability problem for a hybrid system $`H=(Q,X,\mathrm{\Sigma },G,E,Init,f)`$, where $`Init`$ is not necessarily specified in advance, may be posed as follows:

| ReachProblem1:For two subsets $`S_1`$ and $`S_2`$ of $`Q\times X`$, is there a finite |
| --- |
| trajectory, $`(q_1,𝐱_1),(q_2,𝐱_2),\mathrm{},(q_N,𝐱_N)`$, for some |
| $`N1`$, such that $`(q_1,𝐱_1)S_1`$ and $`(q_N,𝐱_N)S_2`$ ? |

This may also be rephrased as

| ReachProblem2: | For two subsets $`S_1`$ and $`S_2`$ of $`Q\times X`$, |
| --- | --- |
|  | whether $`Post^{}(S_1)S_2\mathrm{}`$ ? |
|  | equivalently:  whether $`S_1Pre^{}(S_2)\mathrm{}`$ ? |

If the initial set $`Init`$ is specified in advance, then, with $`S_1=Init`$, the reachability problem ReachProblem1 becomes

| ReachProblem3: | For a subset $`S`$ of $`Q\times X`$, |
| --- | --- |
|  | is there a finite execution of $`H`$ |
|  | with final configuration in $`S`$ ? |

There is a counterpart to the reachability problem, called the avoidance problem, which may be posed as follows:

| AvoidProblem: | For two subsets $`S_1`$ and $`S_2`$ of $`Q\times X`$, |
| --- | --- |
|  | whether $`Post^{}(S_1)S_2=\mathrm{}`$ ? |
|  | equivalently:  whether $`S_1Pre^{}(S_2)=\mathrm{}`$ ? |

In the sequel, we restrict our attention to ReachProblem1 or ReachProblem2. We observe that the answers to these questions depend, in general, not only on the hybrid system $`H`$, but also on the class $`𝒞`$ of subsets of $`Q\times X`$ that are under consideration. Informally, the class $`𝒞`$ is required to

1. include a specified class of sets $`𝒮`$, where $`𝒮`$ may consist of the initial set, sets defined by state invariance conditions, those defined by guard conditions, and domains and ranges of reset maps of edges (and also final sets, if specified),
2. be closed under the boolean set operations of union and complimentation, and under the functions $`Post_\sigma `$ and $`Pre_\sigma `$, and
3. have an effective decision procedure for answering questions such as, for two sets $`S_1,S_2𝒞`$, whether $`S_1=S_2`$ or not.

A formal presentation of these notions is beyond the scope of this work, although a brief discussion may be found in Appendix A. For more details, the interested reader may refer to the references on Model Theory, such as . But, for our purposes, we will be content with the following (somewhat informal) definition:

###### Definition 7

An algorithm for the reachability problem is said to be

1. a decision procedure, if the algorithm stops after a finite number of steps with the correct answer, where the answer can be either yes or no.
2. a semi-decision procedure, if the algorithm

1. never stops with an incorrect answer, and
    2. always stops after a finite number of steps with the correct answer, whenever the answer is yes.

A hybrid system $`H`$ is decidable, if there is a decision procedure for the reachability problem for $`H`$.

Obviously, for a hybrid system $`H`$, if there are semi-decision procedures for both reachability problem and avoidance problem, then $`H`$ is decidable.

We describe below a schematic semi-decision procedure for the reachability problem of a hybrid system:

In the later part of this section, we shall be mainly concerned with the computational aspects of the $`Post`$-operator appearing in the while-loop of the above schematic.

### 2.3 Computation of the $`Post`$-operator

We describe here a computational approach to finding $`Post(S)`$, $`SQ\times ^n`$, as in the semi-decision procedure for ReachProblem2 discussed in Sec. 2.2. Recall that, within a discrete state $`qQ`$, the continuous dynamics of $`𝐱`$ are specified by

$$\frac{d𝐱}{dt}=f(q,𝐱)$$
(2)
with the initial conditions $`𝐱(0)=𝐱_0X_0X=^n`$. We begin our discussion with the following definition.

###### Definition 8

Let $`^+`$ denote the set of nonnegative real numbers. Then, for each discrete state $`qQ`$, the flow associated with equation (5) is a function $`\varphi _q:^n\times ^+^n`$, defined as $`\varphi _q(𝐱,t)=\gamma _𝐱(t)`$, $`t0`$ and $`𝐱^n`$, where the function $`\gamma _𝐱:^+^n`$ satisfies

$$\frac{d\gamma _𝐱(t)}{dt}=f(q,\gamma _𝐱(t)),$$

with the initial condition $`\gamma _𝐱(0)=𝐱`$.

Recall the set of reachable phases $`Reach_q^G`$ as defined in (1). In the following definition, the same is defined in terms of the flow function $`\varphi _q`$:

###### Definition 9

Let $`qQ`$ and $`X_0^n`$. Then, the set of reachable continuous phases in the state $`q`$, is the set $`Reach_q^G\left(X_0\right)`$ defined as

$`Reach_q^G(X_0)`$ $`=`$ $`\{𝐱^n:𝐱_0X_0t0\text{ }𝐱=\varphi _q(𝐱_0,t)\text{ and}`$ (3)
$`s0st\varphi _q(𝐱_0,s)G(q)\}`$
We now define another operator, called the projection operator, as follows:

###### Definition 10

Let $`SQ\times X`$. Then, the projection of S on X, is the set $`\pi __X(S)=\{𝐱X:(q,𝐱)S\}`$.

In terms of the operators $`\pi __X`$ and $`Reach_q^G`$, the operator $`Post`$ may be rewritten as

$$Post(S)=\underset{e=(q_e,X_e,\sigma _e,r_e,q_e^{})E,q=q_e}{}_{X_e{\scriptscriptstyle [Reach_q^G\pi __X(S)]}\mathrm{}}\{q_e^{}\}\times [r_eReach_q^G\pi __X(S)],$$
(4)
where, for two operators $`g`$ and $`h`$, $`gh`$ is the composition of the two operators. With regard to computational issues, the set operators $`\pi __X`$ and $`r_e`$ may not pose difficulties (provided $`r_e`$ is suitably specified). Hence the problem reduces to that of computation of $`Reach_q^G`$, and effectively representing the resulting set. Referring to the definition of $`Reach_q^G`$ as in (3), the significant and challenging task is the elimination of quantifiers, whenever possible. But, as shown in , not all theories of the real number system may admit quantifier elimination (see also the discussion presented in  and in Appendix A). This fact provides a motivation for the study of alternative methods (without requiring quantifier elimination method) for computating the set $`Reach_q^G(X_0)`$.

Specifically, in the remaining part of the work, we shall be concerned with computation of the operator $`Reach_q^G`$ and with representation of the resulting sets. Below is a concise description of the questions that we shall be interested in and of the organization of the rest of the paper:

1. How to effectively compute the set operator $`Reach_q^G`$, at least when the global invariance condition is not imposed but a time bound is specified? In this case, we are interested in the set $`Reach_q(X_0,[0,\tau ])`$, $`\tau 0`$, $`X_0^n`$, defined as

$$Reach_q(X_0,[0,\tau ])=\{𝐱^n:𝐱_0X_0t\mathrm{\hspace{0.33em}0}t\tau \text{ and }𝐱=\varphi _q(𝐱_0,t)\}.$$

In Sec. 3, we present a schematic algorithm for computing $`Reach_q(X_0,[0,\tau ])`$. This algorithm is based on a generalization of the method described in .
2. How to compute – either exactly or approximately – the set operator $`Reach_q^G`$ when the global invariance condition is imposed (without time bound). In this case, we have to deal with the set operator $`Reach_q^G`$ as defined in (3), and obviously, it would be best to find an algorithm based on quantifier elimination. But, as mentioned earlier, we do not assume that quantifier elimination method is feasible, hence we may have to find an algorithm that computes an approximation to $`Reach_q^G`$. In Sec. 3, we indicate how to extend the algorithm for computing $`Reach_q(X_0,[0,\tau ])`$, as mentioned above, to compute an under-approximation or an over-approximation (depending on which one is preferred) to the set $`Reach_q^G(X_0)`$.
3. Finally, how to represent the reach sets obtained by the operators $`Reach_q`$ and $`Reach_q^G`$, i.e. the sets $`Reach_q(X_0,[0,\tau ])`$ or $`Reach_q^G(X_0)`$ as the case may be, such that boolean set operations such as union and complementation can be performed efficiently. This requires represention of the reach sets in terms of sets that are simple and easy to handle, such as polyhedral sets and subalgebraic sets. In , an algorithm algorithm for over-approximating the reach sets using polyhedral sets is presented. In Sec. 4, we describe another algorithm for over-approximating the reach sets with polyhedral sets. The algorithm presented in this paper differs from that of , but if these two methods – that described in Sec. 4 and that of  – are used together, then better results of over-approximation may be obtained.

The remaining part of the paper is orgainized as follows. In Sec. 3 and Sec. 4, we restrict our attention to these issues. We illustrate these algorithms with simple examples in Sec. 5, while Sec. 6 summarizes the paper.

## 3 Computational Approaches for Finding Reach Sets

### 3.1 Preliminary Discussion on Reach Sets

Let $`H=(Q,X,\mathrm{\Sigma },G,E,Init,f)`$ be a hybrid system. For $`qQ`$, consider the equation

$$\frac{d𝐱}{dt}=f(q,𝐱)$$
(5)
with the initial conditions $`𝐱(0)=𝐱_0X_0X=^n`$. It is customary to assume that $`f(q,)`$ is Lipschitz continuous in the second variable, in order to ensure existence and uniquness of solutions to the system of differential equations in (5). Further, we assume that the initial set $`X_0`$ is closed. Let, as before, the flow (refer to Defintion 8) associated with (5) be $`\varphi _q(𝐱,t)`$ and let $`Reach_q(X_0,t)=\{\varphi _q(𝐱_0,t):𝐱_0X_0\}`$ be the set of phases reached at time $`t0`$ in state $`q`$ with initial conditions in $`X_0`$. Further, let

$`Reach_q(X_0,[0,\tau ])`$ $`=`$ $`{\displaystyle \underset{0t\tau }{}}Reach_q(X_0,t),\tau 0,`$
$`Reach_q(X_0,[0,\mathrm{}))`$ $`=`$ $`{\displaystyle \underset{\tau >0}{}}Reach_q(X_0,[0,\tau ]).`$
We shall now define, for a set $`X_q`$, such that $`X_0X_q^n`$, another set $`Reach_q^{}(X_0,X_q,[0,\mathrm{}))`$. For this purpose, we first define, for each $`𝐱X_0`$, a number $`\tau _𝐱`$ as follows:

$$\tau _𝐱=inf\{t>0:\varphi _q(𝐱,t)X_q\},$$
(6)
where we assume that if $`\tau _𝐱=\mathrm{}`$, then $`\varphi _q(𝐱,\mathrm{})`$ denotes the set of all $`\omega `$-limit points (see, for instance, ) of $`\varphi _q(𝐱,t)`$, for $`t0`$, and

$$\{\varphi _q(𝐱,t):0t\mathrm{}\}=\{\varphi _q(𝐱,t):0t<\mathrm{}\}\varphi (𝐱,\mathrm{}).$$

Further, let $`\mathrm{\Theta }_q(𝐱,X_q)^n`$ be defined as follows:

$$\mathrm{\Theta }_q(𝐱,X_q)=\{\begin{array}{c}\{\varphi _q(𝐱,t):0t\tau _𝐱\},\text{if }\tau _𝐱<\mathrm{}\text{ and }\varphi (𝐱,\tau _𝐱)X_q\text{,}\hfill \\ \{\varphi _q(𝐱,t):0t\mathrm{}\},\text{if }\tau _𝐱=\mathrm{}\text{ and }\varphi (𝐱,\mathrm{})X_q\text{,}\hfill \\ \{\varphi _q(𝐱,t):0t<\tau _𝐱\},\text{otherwise.}\hfill \end{array}$$

Then the set $`Reach_q^{}(X_0,X_q,[0,\mathrm{}))`$ is defined as

$$Reach_q^{}(X_0,X_q,[0,\mathrm{}))=\underset{𝐱X_0}{}\mathrm{\Theta }_q(𝐱,X_q).$$

Now, by taking $`X_q=G(q)`$ in the above definition, we have

$$Reach_q^G\left(X_0\right)=Reach_q^{}(X_0,X_q,[0,\mathrm{})),$$

where $`Reach_q^G\left(X_0\right)`$ is defined as in (3), reproduced below for convenience:

$`Reach_q^G(X_0)`$ $`=`$ $`\{𝐱^n:𝐱_0X_0t0\text{ }𝐱=\varphi _q(𝐱_0,t)\text{ and}`$
$`s0st\varphi _q(𝐱_0,s)G(q)\}`$
Throughout the rest of the discussion, we fix a state $`q`$, and consider the problem of computing the sets $`Reach_q(X_0,[0,\tau ])`$ and $`Reach_q^{}(X_0,X_q,[0,\mathrm{}))`$.

### 3.2 A Computational Method: Generalized Face-Lifting Algorithm

When $`X_0`$ is suitably specified (such as, for example, a rectangle), this problem may be solved by finding the solutions to (5) with $`𝐱(0)`$ on the boundary of $`X_0`$. This results in evolving the boundary of $`X_0`$. More precisely, let $`S_0`$ be the boundary of $`X_0`$, and $`X(\tau )=X_0_{0t\tau }\{\varphi _q(𝐱_0,t):𝐱_0S_0\}`$. Then $`X(\tau )=Reach_q(X_0,[0,\tau ])`$ (see Appendix A for a proof). This is the idea underlying the computational approach, called face lifting method, described in . In this section, we shall study this in considerably more general setting, and describe an algorithm, called generalized face lifting method. It may be noted that the method described in  does not assume that a global invarience requirement is imposed; hence the algorithm of  computes only $`Reach_q(X_0,[0,\tau ])`$. However, we shall extend our method for computing $`Reach_q(X_0,[0,\tau ])`$ for computing an approximation to the set $`Reach_q^G(X_0)`$, where the approximation can chosen to be either under-approximation or over-approximation.

Further, let $`S_0^+S_0`$ be the set of those boundary points of $`X_0`$, such that the solution with initial point in $`S_0^+`$ extends into $`X_0^c\stackrel{\mathrm{\Delta }}{=}^n\backslash X_0`$ , i.e.,

$$S_0^+=\{𝐱_0S_0:ϵ=ϵ(𝐱_0)>0\text{ such that }\varphi _q(𝐱_0,t)X_0^c,t(0,ϵ)\},$$
(7)
and define $`X^+(\tau )=X_0_{0t\tau }\{\varphi _q(𝐱_0,t):𝐱_0S_0^+\}`$. It turns out that $`X^+(\tau )=X(\tau )=Reach_q(X_0,[0,\tau ])`$ (see Appendix A). Therefore, in order to find $`Reach_q(X_0,[0,\tau ])`$, we have to find only those solutions for which the initial conditions are in $`S_0^+`$.

If $`S_0^+`$ can be found explicitly (by inspection of $`𝐟`$ and $`S_0`$), then the problem reduces justifiably to finding the solutions with initial conditions in $`S_0^+`$. Otherwise, we may have to find a means of obtaining an outer approximation to $`S_0^+`$, i.e., a set $`S_1^+`$ satisfying $`S_0^+S_1^+S_0`$. We suggest a way to find such a set. For this purpose, we assume that $`X_0`$ is specified as follows: there is a continuously differentiable function, $`\mathrm{}:^n`$, such that if $`𝐱\stackrel{}{X_0}`$ then $`\mathrm{}(𝐱)<0`$ (where $`\stackrel{}{X_0}`$ denotes the interior of $`X_0`$, i.e., the largest open set contained in $`X_0`$), and if $`𝐱X_0^c`$ then $`\mathrm{}(𝐱)>0`$. Hence $`\mathrm{}(𝐱)=0`$ on $`S_0`$. Define

$$S_1^+=\{𝐱S_0:\mathrm{}(𝐱)f(q,𝐱)0\}.$$

It may be easily shown that $`S_0^+S_1^+`$ (see Appendix A). Futher, we observe that the set of points for which $`\mathrm{}(𝐱)f(q,𝐱)>0`$, constitutes an inner approximation to $`S_0^+`$. The method of finding an $`S_1^+`$ as described above can be extended easily to the situation where $`X_0`$ is specified as the set of intersection of finitely many sets of the form $`\mathrm{}_i(𝐱)0`$, with the strict inequality in $`\stackrel{}{X_0}`$.

With this notation, we proceed to describe a schematic algorithm for finding $`X^+(\tau )`$ for $`\tau 0`$. In order to exploit the semigroup property of the reach set, i.e., $`Reach_q(X_0,[0,t])=Reach_q(Reach_q(X_0,[0,s]),[0,ts])`$, for any $`t`$ and $`s`$, with $`0st`$, we consider a sequence of time intervals $`[0,\tau _1]`$, $`[\tau _1,\tau _2]`$, $`[\tau _2,\tau _3]`$, etc., where $`0<\tau _1<\tau _2<\tau _3<\mathrm{}`$, and $`\tau _i\mathrm{}`$, as $`i\mathrm{}`$. It may be convenient to choose for some $`\tau >0`$, $`\tau _i=i\tau `$, $`i=0,1,2,\mathrm{}`$, although we do not require such an assumption in the algorithm.

| Procedure for $`Reach_q(X_0,[0,\tau ])`$ |
| --- |
| initialize: $`X^+(0):=X_0`$, $`F_0:=S_1^+`$ and $`i:=0`$ |
| while $`\tau >\tau _i`$ and $`F_i\mathrm{}`$ do |
| if $`\tau <\tau _{i+1}`$ then $`\mathrm{\Delta }_i:=\tau \tau _i`$ else $`\mathrm{\Delta }_i:=\tau _{i+1}\tau _i`$ end if |
| $`T_i:=\{\varphi _q(𝐱,t):t[0,\mathrm{\Delta }_i],𝐱F_i\}`$ |
| /\* this computational step requires special attention! \*/ |
| $`X^+(\tau _i+\mathrm{\Delta }_i):=X^+(\tau _i)T_i`$ |
| $`F_{i+1}:=\{\varphi (𝐱,\mathrm{\Delta }_i):𝐱F_i\}\backslash X^+(\tau _i)`$ |
| $`i:=i+1`$ |
| end while |

In the above schematic, the initialization step “$`F_0:=S_1^+`$” could be replaced with “$`F_0:=S_0^+`$”, as it is not necessary to initialize $`F_0`$ to $`S_1^+`$. However, before proceeding further, it must be mentioned that with reference to this algorithm, we assume that the computational steps “$`T_i:=\{\varphi _q(𝐱,t):t[0,\mathrm{\Delta }_i],𝐱F_i\}`$” and “$`F_{i+1}:=\{\varphi (𝐱,\mathrm{\Delta }_i):𝐱F_i\}\backslash X^+(\tau _i)`$” can be performed effectively.

We now describe a method for extending this procedure to another procedure that computes an approximation to $`Reach_q^{}(X_0,X_q,[0,\mathrm{}))`$. But since $`Reach_q^G(X_0)=Reach^{}(X_0,G(q),[0,\mathrm{}))`$, the procedure to be described below finds an approximation to $`Reach_q^G(X_0)`$, when $`X_q=G(q)`$. The approximation can be chosen to be either under-approximation or over-approximation, depeding on a flag “$`under\mathrm{\_}approximate`$”, passed as input to the algorithm (refer to the schematic algorithm described below).

To facilitate the discussion, we consider the flow $`\psi _q(𝐱,t)`$, $`t0`$, $`𝐱X=^n`$, corresponding to the differential equation

$$\frac{d𝐱}{dt}=f(q,𝐱)$$
(8)
with the initial conditions $`𝐱(0)=𝐱_0X=^n`$. Now since the flow $`\psi _q`$ has the opposite direction to that of $`\varphi _q`$, for any $`t0`$ and two subsets $`X_1`$ and $`X_2`$ of $`^n`$, we have, $`X_2=\varphi _q(X_1,t)`$ if and only if $`X_1=\psi _q(X_2,t)`$. The function $`\psi _q`$ will be used in the procedure for computing an approximation to $`Reach^{}(X_0,X_q,[0,\mathrm{}))`$, for finding, in each iteration, the set of initial conditions in $`F_i`$ corresponding to which the solutions for the time interval $`[0,\mathrm{\Delta }_i]`$ may violate the global invariance requirement. With this, we present the schematic algorithm as follows:

| Procedure for $`Reach_q^{}(X_0,X_q,[0,\mathrm{}))`$ |
| --- |
| precondition: $`X_0X_q`$ |
| boolean input flag: $`under\mathrm{\_}approximate`$ |
| /\* $`under\mathrm{\_}approximate=1`$, if under-approximation is preferred \*/ |
| /\* by default, the procedure over-approximates the reach set \*/ |
| initialize: $`X^+(0):=X_0`$, $`F_0:=S_1^+`$ and $`i:=0`$ |
| while $`F_i\mathrm{}`$ do |
| $`\mathrm{\Delta }_i:=\tau _{i+1}\tau _i`$ |
| $`T_i:=\{\varphi _q(𝐱,t):t[0,\mathrm{\Delta }_i],𝐱F_i\}`$ |
| if $`T_iX_q`$ then /\* global invariance not violated \*/ |
| $`X^+(\tau _i+\mathrm{\Delta }_i):=X^+(\tau _i)T_i`$ |
| $`F_{i+1}:=\{\varphi _q(𝐱,\mathrm{\Delta }_i):𝐱F_i\}\backslash X^+(\tau _i)`$ |
| else /\* global invariance violated by at least one trajectory \*/ |
| $`U_i:=T_i\backslash X_q`$ |
| $`V_i:=\{\psi _q(𝐱,t):t[0,\mathrm{\Delta }_i],𝐱U_i\}`$ |
| $`F_i^{}=F_i\backslash V_i`$ |
| $`T_i^{}:=\{\varphi _q(𝐱,t):t[0,\mathrm{\Delta }_i],𝐱F_i^{}\}`$ |
| if $`under\mathrm{\_}approximate=1`$ then |
| $`X^+(\tau _i+\mathrm{\Delta }_i):=X^+(\tau _i)U_i`$ |
| else |
| $`X^+(\tau _i+\mathrm{\Delta }_i):=X^+(\tau _i)T_i^{}`$ |
| end if |
| $`F_{i+1}:=\{\varphi _q(𝐱,\mathrm{\Delta }_i):𝐱F_i^{}\}\backslash X^+(\tau _i)`$ |
| end if |
| $`i:=i+1`$ |
| end while |

As with the previous algorithm, the initialization step “$`F_0:=S_1^+`$” could be replaced with “$`F_0:=S_0^+`$”, depending on convenience. Termination of this procedure may be guaranteed, if certain assumptions are satisfied. The required assumptions are as follows:

1. the set $`Y_q=\overline{X_q}`$ is compact;
2. $`X_0`$ is closed in $`^n`$; and
3. for every $`𝐱_0X_0`$, $`Y_q`$ does not contain any $`\omega `$-limit points of the trajectory of $`\varphi _q(𝐱_0,t)`$, $`t0`$, more precisely,

$$\varphi _q(𝐱_0,\mathrm{})Y_q=\mathrm{},𝐱_0X_0.$$

Under these assumptions, it can be shown that (see Appendix C) there is a $`\tau _{\mathrm{max}}>0`$, such that for every $`𝐱X_0`$, there is a $`\tau (𝐱)>0`$ with $`0<\tau (𝐱)\tau _{max}`$ and $`\varphi _q(𝐱,\tau (𝐱))X_q`$. Obviously, with such a $`\tau _{\mathrm{max}}`$, we have $`\tau _𝐱\tau _{\mathrm{max}}`$, for each $`𝐱X_0`$, where $`\tau _𝐱`$ is as in (6). Termination of the algorithm may be deduced from this as follows: Consider a sequence of sets $`F_0^{(i)}`$, $`i0`$, defined as $`F_0^{(i)}=\{𝐱S_0:\varphi _q(𝐱,\tau _i)F_i\}`$. It is easy to check that $`F_0^{(i)}`$ is a nonincreasing sequence of sets, i.e., $`F_0^{(i+1)}F_0^{(i)}`$. Let $`𝐱X_0`$, and let $`k`$ be a nonnegative integer such that $`\tau _k\tau _𝐱<\tau _{k+1}`$. Now, during the $`k`$’th iteration of the while-loop, the global invariance condition (the condition that $`T_kX_q`$) is violated, since $`\varphi _q(𝐱,t)X_q`$, for some $`t`$ such that $`\tau _𝐱t<\tau _{k+1}`$. Let $`\tau _𝐱^{}[\tau _𝐱,\tau _{k+1})`$ be a time instant such that $`\varphi _k(𝐱,\tau _𝐱^{})X_0`$, and let $`t_𝐱=\tau _𝐱^{}\tau _k`$. With this choice of $`t_𝐱`$, we have both $`0t_𝐱<\mathrm{\Delta }_k`$ and $`\varphi _q(𝐱,\tau _k+t_𝐱)X_q`$. Thus, $`\varphi _q(𝐱,\tau _k+t_𝐱)U_k=T_k\backslash X_q`$ (refer to the schematic above). Let $`𝐲=\varphi _q(𝐱,\tau _k+t_𝐱)`$, so $`\varphi _q(𝐱,\tau _k)=\psi _q(𝐲,t_𝐱)`$. But, since $`𝐲U_k`$, $`\psi _q(𝐲,t_𝐱)V_k`$ (refer to the schematic). Hence $`\varphi _q(𝐱,\tau _k)V_k`$ and $`\varphi _q(𝐱,\tau _k)F_k^{}=F_k\backslash V_k`$. Therefore $`\varphi _q(𝐱,\tau _{k+1})F_{k+1}`$, implying that $`𝐱F_0^{(k+1)}`$. By the monotonicity of the sequence of sets $`F_0^{(i)}`$, we have $`𝐱F_0^{(i)}`$, for any $`ik+1`$. Finally, let $`m0`$ be an integer such that $`\tau _m\tau _{\mathrm{max}}<\tau _{m+1}`$. Now since $`\tau _𝐱\tau _{\mathrm{max}}`$, $`𝐱F_0^{(m+1)}`$, for any $`𝐱X_0`$. Hence $`F_0^{(m+1)}=\mathrm{}`$, so $`F_{m+1}=\mathrm{}`$, satisfying the terminating condition of the while-loop of the algorithm, after at most $`m`$ number of iterations.

In the rest of the paper, we discuss approximation of the set $`T_i`$ appearing in the second line of the while-loop of the above schematic (see also the schematic for $`Reach(X_0,[0,\tau ])`$), when $`f(q,𝐱)`$ is of the form $`f(q,𝐱)=A_q𝐱`$, where $`A_q`$ is a constant $`n\times n`$ matrix with real entries and the initial set $`X_0`$ is a polyhedron. For convenience, we drop the subscript $`q`$ in $`A_q`$, and deal with the linear system

$$\frac{d𝐱(t)}{dt}=A𝐱,$$

where $`A`$ is a constant $`n\times n`$ matrix with real valued entries.

## 4 Representation of the Reach Sets

### 4.1 Preliminary Discussion

In this section, we discuss representation of the reach set $`Reach_q(X_0,[0,\tau ])`$, such that boolean set operations such as union and complementation can be performed efficiently. The most convenient representation schemes appear to be representation in terms of the following classes of subsets of $`^n`$:

1. Polyhedral Sets : These are the sets which can be written as union of finitely many polyhedra. A polyhedron in $`^n`$ is a set which may be expressed as the intersection of finitely many closed half spaces , i.e., a finite intersection of sets of the form $`\{𝐱^n:𝐚^T𝐱b\}`$, where $`𝐚^n`$ and $`b`$. It may be observed this class of sets corresponds to (and includes) the class of definable sets in the theory of linear inequalities of $``$, i.e., the theory obtained when $`0`$ and $`1`$ are the constant elements, $`+`$ and $``$ are the (binary) function symbols, and $`<`$ is the (binary) relation symbol (see Appendix A).
2. Subalgebraic Sets : These are the sets which can be written as union of finite number of sets defined by polynomial inequalities, i.e., sets that can be expressed as union of finitely many sets, each of which is the intersection of finitely many sets the form $`\{𝐱^n:p(𝐱)c\}`$, where $`p:^n`$ is a polynomial function with real-valued coefficients and $`c`$. This class of sets corresponds to (and includes) the class of sets that are definable in the theory of $``$ when viewed as an ordered field. The theory of ordered field of $``$ is the theory obtained by extending the theory of linear inequalities of $``$ by including a binary function symbol for representing the product of two real numbers, denoted by $``$ (see Appendix A).

It may be observed that, for a general flow function $`\varphi _q`$, the reach set $`Reach_q(X_0,[0,\tau ])`$ may not be exactly representable in any one of these classes, even if $`X_0`$ is. Hence, as an alternative, we may have to settle for an approximation of the reach set by sets that belong to the respective classes. In this context, we shall restrict our attention to over-approximation of the reach set with polyhedral sets.

Specifically, we describe an algorithm for reach set over-approximation with polyhedral sets, when the system dynamics (within a discrete state $`q`$) are specified as

$$\frac{d𝐱}{dt}=A𝐱,𝐱(0)X_0,$$
(9)
where the initial set $`X_0`$ is a polyhedron and $`A`$ is a constant real valued $`n\times n`$ matrix. Based on the discussion of the last section, we consider the problem of approximating the reach set of the solutions of the equation

$$\frac{d𝐱}{dt}=A𝐱,𝐱(0)F_0,$$
(10)
where $`F_0`$ is a face of $`X_0`$.

The main references on approximate reach set computation using polyhedral sets appear to be  and . In , the initial set is assumed to be convex (not necessarily a polyhedron), and for each time instance, each slice of the “reach tube” (i.e., for each $`t[\tau _i,\tau _{i+1}]`$, the set $`\varphi (F_i,t)`$) is approximated by polyhedra. The method described in  assumes $`F_i`$ to be a polyhedron, and approximates the tube for the time interval $`[\tau _i,\tau _{i+1}]`$, i.e., the set $`T_i=_{0t\tau _{i+1}\tau _i}\varphi (F_i,t)`$, by first constructing a polyhedral over-approximation of the reach tube for the time interval and further over-approximating the resulting polyhedron by a “griddy” polyhedron, i.e., a set that may be expressed as a union of unit hypercubes with integer left-most vertices .

Being a convex set, any polyhedron that over-approximates the reach set contains the convex hull of $`T_i`$. An algorithm based on this observation (as given in ) is to find a “bloated” convex hull of $`F_i`$ and $`F_{i+1}`$. This method seems to be simple and straightforward, and gives reasonable approximations in a rather short time (refer to ). But we may note that those faces of the convex hull of $`F_i`$ and $`F_{i+1}`$, for which the solution set $`T_i`$ lies on one side, need not be bloated. By shrinking the convex hull, an under-approximation is obtained by the same method.

If the time steps are constant, i.e., $`\tau _k=k\tau `$, for some $`\tau >0`$, as in , then we have to find approximation $`P_0`$ only for $`T_0`$, since symbolically, $`T_k=e^{(A\tau )}T_{k1}`$, for $`k1`$, and therefore, $`T_k=e^{(kA\tau )}T_0=e^{(A\tau _k)}T_0`$. Thus, $`P_k=e^{(A\tau _k)}P_0`$ is the required approximation, and if $`\lambda _j`$, $`j=1,2,\mathrm{},N`$, are the normals of the faces $`P_0`$, then the normals of the faces of $`P_k`$ are given by $`e^{(A^T\tau _k)}\lambda _j`$. From this, we may infer that in the subsequent iterations of the algorithm of , the convex hull of $`F_i`$ and $`F_{i+1}`$ need not be computed.

### 4.2 Over-approximating the Reach Set with Polyhedra

In this section, we consider the problem of over-approximating the reach sets, when the dynamics are specified as

$$\frac{d𝐱(t)}{dt}=A𝐱(t),t0,$$

with the initial conditions $`𝐱(0)X_0`$. The solution is given explicitly by $`𝐱(t)=e^{(At)}𝐱(0)`$. We assume that the initial set $`X_0`$ is a polyhedron, defined by

$$\lambda _i^T𝐱h_i0,i=1,2,\mathrm{},\kappa ,$$
(11)
where $`\kappa `$ is a fixed positive integer. We assume that $`\lambda _i^n`$, $`|\lambda _i|=1`$ and $`h_i`$, for $`1i\kappa `$. Then the boundary set of $`X_0`$ consists of sets $`C_j`$, $`1j\kappa `$, of the form

$`\lambda _j^T𝐱h_j`$ $`=`$ $`0,\text{and}`$
$`\lambda _i^T𝐱h_i`$ $``$ $`0,\text{ for }1i\kappa \text{ and }ij\text{ }.`$
Further, we consider only those polyhedral initial sets $`X_0`$ and matrices $`A`$ that satisfy the following assumptions:

1. $`X_0`$ is compact.
2. The set $`S_0^+`$, as defined in (7) corresponding to the flow function $`\varphi _q(𝐱,t)=e^{(At)}𝐱`$, can be expressed as $`S_0^+=C_{j_1}C_{j_2}\mathrm{}C_{j_m}`$, where $`1j_1<j_2<\mathrm{}<j_m\kappa `$, such that for each $`i`$ with $`1im`$, $`\lambda _{j_i}^TA𝐱h_{j_i}\delta _i>0`$, for every $`𝐱C_{j_i}`$.

Referring to these assumptions, it may be mentioned that the main limitation of this method appears to be the restriction imposed by assumption $`𝒜_2`$ on the initial set $`X_0`$ and the matrix $`A`$. In contrast, the method described in  does not require such an assumption. But for our purposes, we need this assumption.

Now, let $`j`$ be an index such that $`C_jS_0^+`$, as in assumption ($`𝒜_2`$). For simplicity and definiteness, let us assume that $`j=\kappa `$ and put $`F_0=C_\kappa `$. Hence, $`F_0`$ is described by

$`\lambda _i^T𝐱h_i`$ $``$ $`0,\text{ }1i\kappa 1\text{,   and}`$
$`\lambda _\kappa ^T𝐱h_\kappa `$ $`=`$ $`0.`$
After some manipulations and rearrangements, if necessary, the above system of linear constraints can be transformed into an equivalent system of linear constraints, which also describes $`F_0`$ but is of the following form:

$$\begin{array}{ccc}𝐚_i^T𝐱b_i\hfill & \hfill & 0,\text{ }1ik1\text{,   and }\hfill \\ 𝐚_k^T𝐱b_k\hfill & =\hfill & 0,\hfill \end{array}\}$$
(12)
where $`k\kappa `$, $`𝐚_k=\lambda _\kappa `$, and for $`1i<k`$, $`|𝐚_i|=1`$ and $`𝐚_i^T𝐚_k=0`$. It may be observed that, since $`|\lambda _k|=1`$ initially, we have $`|𝐚_k|=1`$ as well. Also, by assumption ($`𝒜_2`$), there is a $`\delta >0`$, such that for all the vectors $`𝐱F_0`$, $`𝐚_k^TA𝐱\delta `$. Further, we assume that the system of linear constraints in (12) is consistent and that each of the constraints is linearly independent from the remaining constraints.

With this construction, our objective is to describe an algorithm to find a polyhedron $`P_0`$ as an over-approximation to the set

$$T_0=\{e^{(At)}𝐱_0:𝐱_0F_0,t[0,\mathrm{\Delta }]\},$$
(13)
where $`\mathrm{\Delta }`$ is a small positive number for which the the following condition holds:

1. for some $`\delta _0>0`$, $`𝐚_k^TAe^{(At)}𝐱_0\delta _0`$, for all $`𝐱_0F_0`$ and for all $`t[\mathrm{\Delta },\mathrm{\Delta }]`$.

The existence of such a $`\mathrm{\Delta }`$ can be assured by our assumptions ($`𝒜_1`$) and ($`𝒜_2`$). Later, we shall derive an estimate for such a $`\mathrm{\Delta }`$, for any fixed $`\delta _0>0`$ such that $`\delta _0<\delta `$. In comparison, the method of  does not impose any such restrictions on $`\mathrm{\Delta }`$. But, it may be mentioned that, in both cases, the accuracy of approximation of either method depends on how small a value is chosen for the parameter $`\mathrm{\Delta }`$. The two conditions ($`C_2`$) and ($`C_3`$) stated below follow from condition ($`C_1`$):

1. $`𝐚_k^Te^{(At)}𝐱_0b_k>0`$, for all $`𝐱_0F_0`$ and for all $`t`$ such that $`0<t\mathrm{\Delta }`$; and
2. $`𝐚_k^Te^{(At)}𝐱_0b_k<0`$, for all $`𝐱_0F_0`$ and for all $`t`$ such that $`\mathrm{\Delta }t<0`$.

Before proceeding to describe our algorithm, we note that the set $`F_\mathrm{\Delta }=\{e^{(A\mathrm{\Delta })}𝐱_0:𝐱_0F_0\}`$ may be described by the system of linear constraints

$`𝐚_i^T(\mathrm{\Delta })𝐱b_i`$ $``$ $`0,i=1,2,\mathrm{},k1,`$
$`𝐚_k^T(\mathrm{\Delta })𝐱b_k`$ $`=`$ $`0.`$
where $`𝐚_i(\mathrm{\Delta })=e^{(A^T\mathrm{\Delta })}𝐚_i`$, $`1ik`$. Dividing throughout the last equation in the above system by $`|𝐚_k(\mathrm{\Delta })|`$, we get

$`𝐚_i^T(\mathrm{\Delta })𝐱b_i`$ $``$ $`0,i=1,2,\mathrm{},k1,`$
$`𝐛_k^T𝐱b_k^{}`$ $`=`$ $`0,`$
where $`𝐛_k=\frac{𝐚_k(\mathrm{\Delta })}{|𝐚_k(\mathrm{\Delta })|}`$ and $`b_k^{}=\frac{b_k}{|𝐚_k(\mathrm{\Delta })|}`$. Now, we observe that, since $`e^{(A^T\mathrm{\Delta })}`$ is invertible and $`|𝐚_i|=1`$, each $`𝐚_i(\mathrm{\Delta })`$ is nonzero, and since $`𝐚_k`$ is orthogonal to $`𝐚_i`$, for $`1ik1`$, $`𝐚_k(\mathrm{\Delta })`$ and $`𝐚_i(\mathrm{\Delta })`$ are pair-wise independent. So, after subtracting from each of the inequalities of the above system an appropriate constant times the last equation, and after normalization, the above system of linear constraints can be transformed into the following system of linear constraints:

$$\begin{array}{ccc}𝐛_i^T𝐱b_i^{}\hfill & \hfill & 0,i=1,2,\mathrm{},k1,\hfill \\ 𝐛_k^T𝐱b_k^{}\hfill & =\hfill & 0,\hfill \end{array}\}$$
(14)
where, for $`1ik`$, $`|𝐛_i|=1`$, and for $`1ik1`$, $`𝐛_i^T𝐛_k=0`$. Now, by the condition ($`C_3`$), for $`0t<\mathrm{\Delta }`$ and $`𝐱_0F_0`$, $`𝐚_k^T(\mathrm{\Delta })e^{(At)}b_k=𝐚_k^Te^{\left(A(t\mathrm{\Delta })\right)}b_k<0`$. Hence, we also have, for $`0t<\mathrm{\Delta }`$ and $`𝐱_0F_0`$, $`𝐛_k^Te^{(At)}𝐱_0b_k^{}<0`$.

In this notation, we now describe a schematic algorithm, the output of which is a set of $`4k`$ many parameters, consisting of pairs of vectors and constants, $`(\mu _i,c_i)`$ and $`(\nu _i,d_i)`$, $`i=1,2,\mathrm{},2k`$, where for each $`i`$, $`\mu _i,\nu _i^n`$ and $`c_i,d_i`$, such that the polyhedron $`P_0`$ of intersection of the $`4k`$ half-spaces defined by $`L_i=\{𝐱^n:\mu _i^T𝐱c_i0\}`$ and $`L_i^{}=\{𝐱^n:\nu _i^T𝐱d_i0\}`$ is an over-approximating polyhedron for the set $`T_0`$ as defined in (13).

| Schematic Algorithm for Finding Over-approximating Polyhedron |
| --- |
| Reach Set : $`T_0=\{e^{(At)}𝐱_0:t[0,\mathrm{\Delta }],𝐱_0F_0\}`$ (refer to (13)) |
| Output : $`4k`$ vectors $`\mu _i`$, $`\nu _i`$ and $`4k`$ real numbers $`c_i`$, $`d_i`$, $`1i2k`$ |
| /\* the $`4k`$ half-spaces are $`L_i`$ and $`L_i^{}`$, $`1i2k`$, where \*/ |
| /\* $`L_i=\{𝐱^n:\eta _i(𝐱)=\mu _i^T𝐱c_i0\}`$, and  \*/ |
| /\* $`L_i^{}=\{𝐱^n:\eta _i^{}(𝐱)=\nu _i^T𝐱d_i0\}`$  \*/ |
| for $`i=1,2,\mathrm{},k1`$, do /\* to find ($`\mu _i`$, $`c_i`$) and ($`\nu _i`$, $`d_i`$) \*/ |
| $`l_i:=inf\{l:(𝐚_i^T𝐱b_i)l(𝐚_k^T𝐱b_k)0,`$ $`𝐱T_0\}`$; |
| $`\mu _i:=𝐚_il_i𝐚_k`$; $`c_i:=b_il_ib_k`$; |
| $`l_i^{}:=inf\{l^{}:(𝐛_i^T𝐱b_i^{})+l^{}(𝐛_k^T𝐱b_k^{})0,`$ $`𝐱T_0\}`$; |
| $`\nu _i:=𝐛_i+l_i^{}𝐛_k`$; $`d_i=b_i^{}+l_i^{}b_k^{}`$; |
| end for |
| $`\mu _k:=𝐚_k^T`$; $`c_k=b_k`$; |
| $`\nu _k:=𝐛_k^T`$; $`d_k=b_k^{}`$; |
| $`l_k:=inf\{l:(𝐚_k^T𝐱b_k)l0,`$ $`𝐱T_0\}`$; |
| $`\mu _k:=𝐚_k`$; $`c_{k+1}:=b_k+l_k`$; |
| $`l_k^{}:=inf\{l:(𝐛_k^T𝐱b_k^{})l0`$, $`𝐱T_0\}`$; |
| $`\nu _k:=𝐛_k`$; $`d_{k+1}:=b_k^{}+l_k^{}`$; |
| for $`i=1,2,\mathrm{},k1`$, do |
| /\* to find ($`\mu _{i+k+1},c_{i+k+1}`$) and ($`\nu _{i+k+1},d_{i+k+1}`$) \*/ |
| $`l_{i+k}=inf\{l:(𝐚_i^T𝐱b_i)l0,`$ $`𝐱T_0\}`$; |
| $`\mu _{i+k}=𝐚_i`$; $`c_{i+k+1}=b_i+l_{i+k}`$; |
| $`l_{i+k}^{}=inf\{l^{}:(𝐛_i^T𝐱b_i^{})l^{}0,`$ $`𝐱F_0\}`$; |
| $`\nu _{i+k}=𝐛_i`$; $`d_{i+k+1}=b_i^{}+l_{i+k}^{}`$; |
| end for |

We first observe that both $`l_i>\mathrm{}`$ and $`l_i^{}>\mathrm{}`$. If $`𝐚_k^TA𝐱\delta >0`$, as in assumption $`𝒜_2`$, then both $`l_i<\mathrm{}`$ and $`l_i^{}<\mathrm{}`$, as we shall see in Sec. 4.3 below. If $`\mathrm{\Delta }`$ is small, such that the set $`T_0`$ is nearly a polyhedron, then the above method gives reasonable results. Referring to Step 1 of the first for-loop of the above algorithm, the reason for choosing $`l_i`$ to be the infimum over all $`l`$ for which $`(𝐚_i^T𝐱b_i)l(𝐚_k^T𝐱b_k)0`$, $`𝐱T_0`$, is obvious: if $`P_1`$ and $`P_2`$ are two polyhedra such that $`P_1`$ is obtained by the above algorithm and $`P_2`$ is obtained by replacing a constraint $`\mu _i^T𝐱c_i0`$, where $`\mu _i`$ and $`c_i`$ are as in Step 2, with another constraint $`\lambda _i^T𝐱h_i0`$, where $`\lambda _i=𝐚_il𝐚_k`$ and $`h_i=b_ilb_k`$, for some $`l>l_i`$, then $`P_1P_2`$.333This follows from the observation that if both $`(𝐚_i^T𝐱b_i)l(𝐚_k^T𝐱b_k)0`$ and $`(𝐚_k^T𝐱b_k)=\mu _k^T𝐱c_k0`$, where $`l`$ and $`𝐱^n`$, $`𝐱`$ not necessarily in $`T_0`$, then, for any $`h>l`$, $`(𝐚_i^T𝐱b_i)h(𝐚_k^T𝐱b_k)0`$. (An analogous statement holds for each of the remaining parameters, $`l_i^{}`$, $`l_k`$, $`l_k^{}`$, $`l_{k+i}`$ and $`l_{k+i}^{}`$, as in the algorithm.)

It can be shown that the polyhedron included in the $`k+1`$ half-spaces, specified by $`L_1,L_2,\mathrm{},L_k,L_{k+1}`$, is a bounded polyhedron (see Appendix D). We may also note that the hyperplanes $`\mu _i^T𝐱c_i=0`$ and $`\nu _i^T𝐱d_i=0`$, $`1ik1`$, are obtained by rotating the hyperplanes $`𝐚_i^T𝐱b_i=0`$ and $`𝐛_i^T𝐱b_i^{}=0`$ about their corresponding intersection with the hyperplanes $`𝐚_k^T𝐱b_k=0`$ and $`𝐛_k^T𝐱b_k^{}=0`$, respectively; whereas the hyperplanes $`\mu _{k+1}^T𝐱c_{k+1}=0`$ and $`\nu _{k+1}^T𝐱d_{k+1}=0`$ are obtained by translating the hyperplanes $`𝐚_k^T𝐱b_k=0`$ and $`𝐛_k^T𝐱b_k=0`$, respectively. Similarly, for $`i=1,2,\mathrm{},k1`$, the hyperplanes $`\mu _{k+i}^T𝐱c_{k+i}=0`$ and $`\nu _{k+i}^T𝐱d_{k+i}=0`$ are translations of the hyperplanes $`𝐚_i^T𝐱b_i=0`$ and $`𝐛_i^T𝐱b_i^{}=0`$, respectively.

In remaining part of the section (in Sec. 4.3 below), we shall derive upper bounds for the numbers $`l_i`$ and $`l_i^{}`$, $`1i2k1`$, that are defined in the schematic. If all the $`l_i`$ and $`l_i^{}`$ are replaced with their corresponding upper bounds, $`\widehat{l}_i`$ and $`\widehat{l}_i^{}`$ respectively, then we obtain conservative estimates – $`\widehat{\mu }_i`$, $`\widehat{c}_i`$, $`\widehat{\nu }_i`$ and $`\widehat{d}_i`$ respectively – for $`\mu _i`$, $`c_i`$, $`\nu _i`$ and $`d_i`$. After this replacement, we obtain a modified algorithm for over-approximation with polyhedra. It may be observed that these estimates overcome the difficulties in finding the infima required for obtaining $`l_i`$’s and $`l_i^{}`$’s. But unfortunately, the conservative upper bounds that we derive for $`l_i`$’s and $`l_i^{}`$’s (as in Sec. 4.3) may turn out to be very large. However, better accuracy may be obtained if the modified alogrithm is used in conjunction with that of  (see Fig. 5 in Sec. 5). More precisely, the intersection of the polyhedron obtained by the method described here with that obtained by the method of  gives a smaller over-approximating polyhedron, as will be discussed in Sec. 5, while illustrating these algorithms with simple examples.

### 4.3  Upper Bounds for $`l_i`$ and $`l_i^{}`$

In this section, we derive upper bounds for the numbers $`l_i`$ and $`l_i^{}`$ appearing in the schematic algorithm. We begin with the following result:

1. Claim 1. For $`i=1,2,\mathrm{},k`$, with $`𝐚_i`$ and $`b_i`$ as in (12) and $`𝐱_0F_0`$ and $`t(0,\mathrm{\Delta }]`$, we have

$$\left|\frac{(𝐚_i^Te^{(At)}𝐱_0𝐚_i^T𝐱_0)}{t}\right|M_0𝐚_i^TAe^{(A\mathrm{\Delta })},$$

where $`M_0=\mathrm{max}_{𝐱_0F_0}\{𝐱_0\}`$.

To prove the claim, we first note that $`(𝐚_i^Te^{(At)}𝐱_0𝐚_i^T𝐱_0)`$ $`=t\times 𝐚_i^TAe^{(A\theta )}𝐱_0`$, for some $`\theta `$ with $`0<\theta <t`$, where $`\theta `$ may depend on $`t`$. Hence, we have

$`\left|{\displaystyle \frac{(𝐚_i^Tee^{(At)}𝐱_0𝐚_i^T𝐱_0)}{t}}\right|`$ $`=`$ $`|𝐚_i^TAe^{(A\theta )}𝐱_0|,\text{where }\theta \text{ is such that }0<\theta <t`$
$``$ $`|𝐚_i^TAe^{(A\theta )}|\times 𝐱_0`$
$``$ $`M_0𝐚_i^TAe^{(A\mathrm{\Delta })}.`$
With $`𝐛_i`$ and $`b_i^{}`$, we have an analogous result, but with a slight modification, as in the following:

1. Claim 2. For $`i=1,2,\mathrm{},k`$, with $`𝐛_i`$ and $`b_i^{}`$ as in (12) and $`𝐱_0F_0`$ and $`t[0,\mathrm{\Delta })`$, we have

$$\left|\frac{(𝐛_i^Te^{(At)}𝐱_0𝐛_i^Te^{(A\mathrm{\Delta })}𝐱_0)}{(\mathrm{\Delta }t)}\right|M_0𝐛_i^TAe^{(A\mathrm{\Delta })},$$

where $`M_0=\mathrm{max}_{𝐱_0F_0}\{𝐱_0\}`$.

We have $`(𝐛_i^Te^{(At)}𝐱_0𝐛_i^Te^{(A\mathrm{\Delta })}𝐱_0)`$ $`=(t\mathrm{\Delta })\times 𝐛_i^TAe^{(A\theta )}𝐱_0`$, for some $`\theta `$ with $`t<\theta <\mathrm{\Delta }`$. Hence

$`\left|{\displaystyle \frac{(𝐛_i^Te^{(At)}𝐱_0𝐛_i^Te^{(A\mathrm{\Delta })}𝐱_0)}{(\mathrm{\Delta }t)}}\right|`$ $`=`$ $`|𝐛_iAe^{(A\theta )}𝐱_0|,\text{ for some }\theta \text{ with }t<\theta <\mathrm{\Delta }`$
$``$ $`𝐛_i^TAe^{(A\theta )}\times 𝐱_0`$
$``$ $`M_0𝐛_i^TAe^{(A\mathrm{\Delta })}.`$
Finally, we shall make the following claim, before deriving upper bounds for $`l_i`$’s and $`l_i^{}`$’s:

1. Claim 3. Assume, for some $`\delta >0`$, $`𝐚_k^TA𝐱_0\delta `$, for every $`𝐱_0F_0`$, and let $`\delta _0`$ be such that $`0<\delta _0<\delta `$. Then there is a $`\mathrm{\Delta }>0`$ such that $`𝐚_k^TAe^{(At)}𝐱_0\delta _0`$, for all $`𝐱_0F_0`$ and for all $`t[\mathrm{\Delta },\mathrm{\Delta }]`$.

In order to prove the claim, we shall derive a conservative estimate for a $`\mathrm{\Delta }>0`$ for which the claim holds. Fix a $`\delta _0`$ such that $`0<\delta _0<\delta `$. We have

$`𝐚_k^TAe^{(At)}𝐱_0`$ $`=`$ $`𝐚_k^TA𝐱_0+{\displaystyle _0^t}𝐚_k^TA^2e^{(As)}𝐱_0ds,\text{ where }𝐱_0F_0\text{,}`$
$``$ $`𝐚_k^TA𝐱_0{\displaystyle _0^{|t|}}𝐚_k^TA^2e^{(As)}𝐱_0ds`$
$``$ $`𝐚_k^TA𝐱_0M_0A(e^{A|t|}1),`$
where $`M_0=\mathrm{max}_{𝐱_0F_0}\{𝐱_0\}`$. Therefore, if $`\mathrm{\Delta }`$ is such that

$$M_0A(e^{A\mathrm{\Delta }}1)\delta \delta _0,$$

then the claim holds.

#### Upper Bounds for $`l_i`$ and $`l_i^{}`$, $`1ik1`$.

Recall that, for $`1ik1`$,

$$l_i=inf\{l:(𝐚_i^T𝐱b_i)l(𝐚_k^T𝐱b_k)0,𝐱T_0\}.$$
(15)
Referring to the right hand side of (15), for every $`𝐱_0F_0`$, since $`𝐚_k^T𝐱_0b_k=0`$, for any real number $`l`$, we have $`(𝐚_i^T𝐱_0b_i)l(𝐚_k^T𝐱_0b_k)=(𝐚_i^T𝐱_0b_i)0`$. Therefore, we may assume $`𝐱T_0\backslash F_0`$, and $`(𝐚_k^T𝐱b_k)>0`$. Thus $`l_i`$ has to be chosen such that

$$l_i=\underset{𝐱T_0\backslash F_0}{sup}\frac{(𝐚_i^T𝐱b_i)}{(𝐚_k^T𝐱b_k)}.$$

Rewriting the above, we have

$$l_i=\underset{0<t\mathrm{\Delta }}{sup}_{𝐱_0F_0}\frac{(𝐚_i^Te^{(At)}𝐱_0b_i)}{(𝐚_k^Te^{(At)}𝐱_0b_k)}.$$

Dividing the numerator and denominator by $`t`$, this may be written as

$$l_i=\underset{0<t\mathrm{\Delta }}{sup}_{𝐱_0F_0}\frac{\alpha _i(𝐱_0,t)}{\beta _i(𝐱_0,t)},$$

where $`\alpha _i(𝐱_0,t)`$ and $`\beta _i(𝐱_0,t)`$ are given by

$$\alpha _i(𝐱_0,t)=\frac{(𝐚_i^Te^{(At)}𝐱_0b_i)}{t}\text{ and }\beta _i(𝐱_0,t)=\beta (𝐱_0,t)=\frac{(𝐚_k^Te^{(At)}𝐱_0b_k)}{t}.$$

But, since $`𝐚_i^T𝐱_0b_i0`$, for every $`𝐱_0F_0`$, we have

$$\frac{(𝐚_i^Te^{(At)}𝐱_0b_i)}{t}\frac{(𝐚_i^Te^{(At)}𝐱_0𝐚_i^T𝐱_0)}{t}\left|\frac{(𝐚_i^Te^{(At)}𝐱_0𝐚_i^T𝐱_0)}{t}\right|.$$

By Claim 1, we have

$$\underset{0<t\mathrm{\Delta }}{sup}_{𝐱_0F_0}\alpha _i(𝐱_0,t)M_0𝐚_i^TAe^{(A\mathrm{\Delta })}.$$

As to the denominator, since $`𝐚_k^T𝐱_0b_k=0`$, for every $`𝐱_0F_0`$, we have

$$\frac{(𝐚_k^Te^{(At)}𝐱_0b_k)}{t}=\frac{(𝐚_k^Te^{(At)}𝐱_0𝐚_k^T𝐱_0)}{t}=𝐚_k^TAe^{(A\theta )}𝐱_0\delta _0,$$

where, $`\theta `$ is some number in the interval $`(0,t)`$, and the last inequality is due to Claim 3. Hence

$$\underset{0<t\mathrm{\Delta }}{inf}_{𝐱_0F_0}\beta (𝐱_0,t)\delta _0>0.$$

Combining both inequalities, we have

$$l_i\frac{sup_{0<t\mathrm{\Delta },𝐱_0F_0}\alpha _i(𝐱_0,t)}{inf_{0<t\mathrm{\Delta },𝐱_0F_0}\beta _i(𝐱_0,t)}\widehat{l}_i=\frac{M_0𝐚_i^TAe^{(A\mathrm{\Delta })}}{\delta _0}.$$

Likewise, $`l_i^{}`$ must be chosen such that

$$l_i^{}=\underset{0t<\mathrm{\Delta }}{sup}_{𝐱_0F_0}\frac{\alpha _i^{}(𝐱_0,t)}{\beta _i^{}(𝐱_0,t)},$$

where $`\alpha _i^{}(𝐱_0,t)`$ and $`\beta _i^{}(𝐱_0,t)`$ are given by

$$\alpha _i^{}(𝐱_0,t)=\frac{(𝐛_i^Te^{(At)}𝐱_0b_i^{})}{(\mathrm{\Delta }t)}\text{ and }\beta _i^{}(𝐱_0,t)=\beta ^{}(𝐱_0,t)=\frac{(𝐛_k^Te^{(At)}𝐱_0b_k^{})}{(\mathrm{\Delta }t)}.$$

Now, since $`𝐛_i^Te^{(A\mathrm{\Delta })}𝐱_0b_i^{}0`$, for $`𝐱_0F_0`$,

$`{\displaystyle \frac{(𝐛_i^Te^{(At)}𝐱_0b_i^{})}{(\mathrm{\Delta }t)}}`$ $``$ $`{\displaystyle \frac{(𝐛_i^Te^{(At)}𝐱_0𝐛_i^Te^{(A\mathrm{\Delta })}𝐱_0)}{(\mathrm{\Delta }t)}}`$
$``$ $`\left|{\displaystyle \frac{(𝐛_i^Te^{(At)}𝐱_0𝐛_i^Te^{(A\mathrm{\Delta })})}{(\mathrm{\Delta }t)}}\right|`$
$``$ $`M_0𝐛_i^TAe^{(A\mathrm{\Delta })},\text{by }\text{Claim 2}.`$
For the denominator function, since $`𝐛_k^Te^{(A\mathrm{\Delta })}𝐱_0b_k^{}=0`$, for every $`𝐱_0F_0`$, we have

$$\frac{(𝐛_k^Te^{(At)}𝐱_0b_k^{})}{(\mathrm{\Delta }t)}=\frac{(𝐛_k^Te^{(At)}𝐱_0𝐛_k^Te^{(A\mathrm{\Delta })}𝐱_0)}{(\mathrm{\Delta }t)}=𝐛_k^TAe^{(A\theta )}𝐱_0,$$

where $`\theta `$ is some number in the interval $`(t,\mathrm{\Delta })`$. Now,

$$𝐛_k^TAe^{(A\theta )}𝐱_0=\frac{𝐚_k^Te^{(A\mathrm{\Delta })}Ae^{(A\theta )}𝐱_0}{𝐚_k^Te^{(A\mathrm{\Delta })}}\delta _1=\frac{\delta _0}{𝐚_k^Te^{(A\mathrm{\Delta })}}.$$

Hence, we have the following upper bound for $`l_i^{}`$:

$$l_i^{}\frac{sup_{0<t\mathrm{\Delta },𝐱_0F_0}\alpha _i^{}(𝐱_0,t)}{inf_{0<t\mathrm{\Delta },𝐱_0F_0}\beta _i^{}(𝐱_0,t)}\widehat{l}_i^{}=\frac{M_0𝐚_i^TAe^{(A\mathrm{\Delta })}}{\delta _1}.$$

Thus, we obtain the following upper bounds:

$`l_i`$ $``$ $`\widehat{l}_i={\displaystyle \frac{M_0𝐚_i^TAe^{(A\mathrm{\Delta })}}{\delta _0}},1ik1`$
$`l_i^{}`$ $``$ $`\widehat{l}_i^{}={\displaystyle \frac{M_0𝐚_i^TAe^{(A\mathrm{\Delta })}}{\delta _1}},1ik1,\text{ where}`$
$`\delta _1={\displaystyle \frac{\delta _0}{𝐚_k^Te^{(A\mathrm{\Delta })}}}.`$
#### Upper Bounds for $`l_k`$, $`l_k^{}`$, $`l_{k+i}`$ and $`l_{k+i}^{}`$, $`1ik1`$.

We have $`l_{i+k}=inf\{l:(𝐚_i^T𝐱b_i)l0,𝐱T_0\}`$, and $`l_k=inf\{l:(𝐚_k^T𝐱b_k)l0,𝐱T_0\}`$. Let $`l_{2k}=l_k`$, so we may consider $`l_{k+i}`$, for $`1ik`$. We have to choose $`l_{k+i}`$ such that

$$l_{k+i}=\underset{𝐱T_0}{sup}(𝐚_i^T𝐱b_i)=\underset{0<t\mathrm{\Delta }}{sup}_{𝐱_0F_0}(𝐚_i^Te^{(At)}𝐱_0b_i).$$

Now, since $`(𝐚_i^T𝐱_0b_i)0`$, for every $`𝐱_0F_0`$, we have

$`\underset{0<t\mathrm{\Delta }}{sup}_{𝐱_0F_0}(𝐚_i^Te^{(At)}𝐱_0b_i)`$ $``$ $`\underset{0<t\mathrm{\Delta }}{sup}_{𝐱_0F_0}(𝐚_i^Te^{(At)}𝐱_0𝐚_i^T𝐱_0)`$
$``$ $`\underset{0<t\mathrm{\Delta }}{sup}_{𝐱_0F_0}t\times {\displaystyle \frac{(𝐚_i^Te^{(At)}𝐱_0𝐚_i^T𝐱_0)}{t}}`$
$``$ $`\mathrm{\Delta }\times \underset{0<t\mathrm{\Delta }}{sup}_{𝐱_0F_0}|{\displaystyle \frac{(𝐚_i^Te^{(At)}𝐱_0𝐚_i^T𝐱_0)}{t}}|`$
$``$ $`\mathrm{\Delta }\times M_0𝐚_i^TAe^{(A\mathrm{\Delta })}.`$
Similarly, for $`1ik1`$, we have to choose $`l_{k+i}^{}`$ as follows

$$l_{k+i}^{}=\underset{𝐱T_0}{sup}(𝐛_i^T𝐱b_i^{})=\underset{0t<\mathrm{\Delta }}{sup}_{𝐱_0F_0}(𝐛_i^Te^{(At)}𝐱_0b_i^{}).$$

A calculation similar to the above shows

$$\underset{0t<\mathrm{\Delta }}{sup}_{𝐱_0F_0}(𝐛_i^Te^{(At)}𝐱_0b_i^{})\mathrm{\Delta }\times M_0𝐛_i^TAe^{(A\mathrm{\Delta })}.$$

Finally, for $`l_k^{}`$, we have

$$l_k^{}=\underset{𝐱T_0}{sup}[(𝐛_k^T𝐱b_k^{})]=\underset{0t<\mathrm{\Delta }}{sup}_{𝐱_0F_0}(𝐛_k^Te^{(A\mathrm{\Delta })}𝐱_0𝐛_k^Te^{(At)}𝐱_0)$$

Another sequence of similar calculations shows

$$\underset{0t<\mathrm{\Delta }}{sup}_{𝐱_0F_0}(𝐛_k^Te^{(A\mathrm{\Delta })}𝐱_0𝐛_k^Te^{(At)}𝐱_0)\mathrm{\Delta }\times M_0𝐛_k^TAe^{(A\mathrm{\Delta })}.$$

So, to collect all the estimates, we have

$`l_k`$ $``$ $`\widehat{l}_k=M_0\mathrm{\Delta }𝐚_k^TAe^{(A\mathrm{\Delta })}`$
$`l_k^{}`$ $``$ $`\widehat{l}_k^{}=M_0\mathrm{\Delta }𝐛_k^TAe^{(A\mathrm{\Delta })}`$
$`l_{k+i}`$ $``$ $`\widehat{l}_{k+i}=M_0\mathrm{\Delta }𝐚_i^TAe^{(A\mathrm{\Delta })},1ik1`$
$`l_{k+i}^{}`$ $``$ $`\widehat{l}_{k+i}^{}=M_0\mathrm{\Delta }𝐛_i^TAe^{(A\mathrm{\Delta })},1ik1`$
## 5  Illustration

### 5.1 Example 1

We first illustrate the schematic algorithm presented in Sec. 3 with the help of an example taken from . Let

$$\begin{array}{ccc}\dot{x}\hfill & =& a\hfill \\ \dot{y}\hfill & =& b\hfill \end{array}\}a,b>0,(x(0),y(0))X_0=[0,1]\times [0,1].$$

So, for $`t0`$, $`x(t)=x(0)+at`$ and $`y(t)=y(0)+bt`$. It is easy to check that $`S_0^+=\{1\}\times [0,1][0,1]\times \{1\}`$. Therefore

$$X^+(\tau )=X_0\underset{0t\tau }{}\{(1+at,c+bt):0c1\}\underset{0t\tau }{}\{(c+at,1+bt):0c1\}.$$

This is illustrated in Fig. 1.

### 5.2 Example 2

In this example, we illustrate the over-approximation algorithm. Let $`𝐱(t)=[x_1(t),x_2(t)]^T^2`$ satisfy

$$\begin{array}{ccc}\dot{x}_1\hfill & =& x_2\hfill \\ \dot{x}_2\hfill & =& x_1\hfill \end{array}\}(x_1(0),x_2(0))F_0=[1,\sqrt{2}]\times \{0\}.$$

The solution is given explicitly by

$$\left[\begin{array}{c}x_1(t)\\ x_2(t)\end{array}\right]=\left[\begin{array}{cc}\mathrm{cos}t& \mathrm{sin}t\\ \mathrm{sin}t& \mathrm{cos}t\end{array}\right]\left[\begin{array}{c}x_1(0)\\ x_2(0)\end{array}\right]$$

With $`\mathrm{\Delta }=\pi /6`$, for the time interval $`[0,\pi /6]`$, the solution set in a parametric form is $`T_0=\{(a\mathrm{cos}t,a\mathrm{sin}t):a[1,\sqrt{2}],\mathrm{\hspace{0.33em}0}t\pi /6\}`$. For this example, in the notation of Sec. 4, we have

$$𝐚_1=\left[\begin{array}{c}1\\ 0\end{array}\right],𝐚_2=\left[\begin{array}{c}1\\ 0\end{array}\right],𝐚_3=\left[\begin{array}{c}0\\ 1\end{array}\right],b_1=\sqrt{2},b_2=1,\text{ and }b_3=0,$$

and

$$𝐛_1=\left[\begin{array}{c}\frac{\sqrt{3}}{2}\\ 0.5\end{array}\right],𝐛_2=\left[\begin{array}{c}\frac{\sqrt{3}}{2}\\ 0.5\end{array}\right],𝐛_3=\left[\begin{array}{c}0.5\\ \frac{\sqrt{3}}{2}\end{array}\right],b_1^{}=\sqrt{2},b_2^{}=1,\text{ and }b_3^{}=0.$$

The set $`T_0`$ and the result of the algorithm with $`l_i`$ and $`l_i^{}`$ exactly found as in the schematic of Sec. 4 are shown in Fig. 2, where the horizontal axis corresponds to $`x_1`$ and the vertical axis corresponds to $`x_2`$.

Now, recalling the notation of Sec. 4, we may take $`\delta =1`$, $`\delta _0=\frac{\sqrt{3}}{2}`$. It is easy to check that $`M_0=\sqrt{2}`$ and $`𝐚_i^TA=𝐛_i^TA=1`$, for $`1i3`$. Hence

$$l_i\widehat{l}=\frac{\sqrt{2}\times 1\times e^{\pi /6}}{\sqrt{3}/2}=\frac{2\sqrt{2}\times e^{\pi /6}}{\sqrt{3}}2.7566424$$

and the same upper bound holds for $`l_i^{}`$, since $`\delta _1=\delta _0`$ in this example. For $`l_k`$ and $`l_k^{}`$ we have the following bounds:

$$l_k=l_k^{}\widehat{l}_k=\sqrt{2}\times \frac{\pi }{6}\times e^{\pi /6}1.249999$$

In this example, the same upper bound, given by $`\widehat{l}_k`$, holds for $`l_{k+i}`$ and $`l_{k+i}^{}`$. Combining all this, the polyhedron is the intersection of the half-spaces given by $`\eta _i(𝐱)0`$ and $`\eta _i^{}(𝐱)0`$, $`1i6`$, where

$`\eta _1(𝐱)=x_12.7566424x_2\sqrt{2},`$ $`\eta _1^{}(𝐱)=0.5122958x_1+2.8873223x_2\sqrt{2},`$
$`\eta _2(𝐱)=x_12.7566424x_2+1.0,`$ $`\eta _2^{}(𝐱)=2.2443466x_1+1.8873223x_2+1.0,`$
$`\eta _3(𝐱)=x_2,`$ $`\eta _3^{}(𝐱)=x_1+\sqrt{3}x_2,`$
$`\eta _4(𝐱)=x_21.249999,`$ $`\eta _4^{}(𝐱)=0.5x_1{\displaystyle \frac{\sqrt{3}}{2}}x_21.249999.`$
$`\eta _5(𝐱)=x_1\sqrt{2}1.249999,`$ $`\eta _5^{}(𝐱)={\displaystyle \frac{\sqrt{3}}{2}}x_1+0.5x_2\sqrt{2}1.249999,`$
$`\eta _6(𝐱)=x_1+0.249999,`$ $`\eta _6^{}(𝐱)={\displaystyle \frac{\sqrt{3}}{2}}x_10.5x_2+0.249999.`$
In this example, $`\eta _4^{}(𝐱)0`$, $`\eta _6(𝐱)0`$ and $`\eta _6^{}(𝐱)0`$, are redundant.

Fig. 3 illustrates the result of the algorithm with the first $`8`$ half-spaces defined by $`\eta _i(𝐱)0`$ and $`\eta _i^{}(𝐱)0`$, $`1i4`$, where the dashed lines correspond to the lines $`\eta _i(𝐱)=0`$ and $`\eta _i^{}(𝐱)=0`$, and the polyhedron included in their intersection is shown in thick lines. The point of intersection of $`\eta _1(𝐱)=0`$ and $`\eta _4(𝐱)=0`$ is $`(4.8600138,\mathrm{\hspace{0.17em}1.249999})`$, of $`\eta _4(𝐱)=0`$ and $`\eta _1^{}(𝐱)=0`$ is $`(4.2845099,\mathrm{\hspace{0.17em}1.249999})`$, and of $`\eta _2(𝐱)=0`$ and $`\eta _2^{}(𝐱)=0`$ is $`(0.575162,\mathrm{\hspace{0.17em}0.154114})`$. So the vertices in the counter-clockwise order are given by $`(\sqrt{2},\mathrm{\hspace{0.17em}0}),`$ $`(4.8600138,\mathrm{\hspace{0.17em}1.249999}),`$ $`(4.2845099,\mathrm{\hspace{0.17em}1.249999}),`$ $`(\sqrt{3}/\sqrt{2},\mathrm{\hspace{0.17em}1}/\sqrt{2}),`$ $`(\sqrt{3}/2,\mathrm{\hspace{0.17em}0.5}),`$ $`(0.575162,\mathrm{\hspace{0.17em}0.154114})`$ and $`(1,0)`$.

Fig. 4 shows the result when the remaining two half-spaces defined by $`\eta _5(𝐱)0`$ and $`\eta _5^{}(𝐱)0`$ are also used. Finally, Fig. 5 shows the intersection of the polyhedron obtained by the method reported in this paper with that which may possibly be obtained by the method of , the latter being shown in dotted lines. The polyhedron corresponding to the method of  is obtained as follows: the convex hull of the sets $`[1,\sqrt{2}]\times \{0\}`$ and $`\{(a\mathrm{cos}\frac{\pi }{6},a\mathrm{sin}\frac{\pi }{6}):a[1,\sqrt{2}]\}`$ is the polygon with vertices in the counter-clockwise order $`(\sqrt{2},0)`$, $`(\sqrt{3}/\sqrt{2},1/\sqrt{2})`$, $`(\sqrt{3}/2,1/2)`$ and $`(1,0)`$. Therefore the polygon is the interesction of the half-spaces $`\zeta _i(𝐱)0`$ and $`\zeta _i^{}(𝐱)0`$, $`i=1,2`$, where

$`\zeta _1(𝐱)=x_2,`$ $`\zeta _1^{}(𝐱)=x_1+\sqrt{3}x_2,\text{ and}`$
$`\zeta _2(𝐱)={\displaystyle \frac{1}{\sqrt{2}}}(x_1\sqrt{2})\left({\displaystyle \frac{\sqrt{3}}{\sqrt{2}}}\sqrt{2}\right)x_2,`$ $`\zeta _2^{}(𝐱)={\displaystyle \frac{1}{2}}(x_11)+\left({\displaystyle \frac{\sqrt{3}}{2}}1\right)x_2.`$
Now if $`\widehat{ϵ}`$ is an upper bound for the bloating parameter $`ϵ`$, then the half-spaces of the over-approximating polyhedron corresponding to the method of  are given by $`\widehat{\zeta }_i(𝐱)0`$ and $`\widehat{\zeta }_i^{}(𝐱)0`$, $`i=1,2`$, where

$`\widehat{\zeta }_1(𝐱)`$ $`=`$ $`x_2+\widehat{ϵ},`$
$`\widehat{\zeta }_1^{}(𝐱)`$ $`=`$ $`x_1+\sqrt{3}x_2+\widehat{ϵ},`$
$`\widehat{\zeta }_2(𝐱)`$ $`=`$ $`0.70710678x_1+0.18946869x_21+\widehat{ϵ},\text{ and}`$
$`\widehat{\zeta }_2^{}(𝐱)`$ $`=`$ $`0.70710678x_10.1339746x_2+0.70710678+\widehat{ϵ}`$
The upper bound for the bloating parameter $`ϵ`$ as given in  works out to be

$`ϵ`$ $``$ $`M_0\left(e^{(A\mathrm{\Delta })}1A\mathrm{\Delta }{\displaystyle \frac{3}{8}}A^2\mathrm{\Delta }^2\right)`$
$`=`$ $`\sqrt{2}\times \left(e^{(\pi /6)}1{\displaystyle \frac{\pi }{6}}{\displaystyle \frac{3\pi ^2}{8\times 36}}\right)`$
$`=`$ $`\sqrt{2}\times 0.061684640.087235255<0.09.`$
Fig. 5 shows the results when an upper bound for the bloating parameter is chosen to be $`\widehat{ϵ}0.2`$, where the dotted lines correspond to the line $`\widehat{\zeta }_1(𝐱)=0`$, $`\widehat{\zeta }_1^{}(𝐱)=0`$, $`\widehat{\zeta }_2(𝐱)=0`$ and $`\widehat{\zeta }_2^{}(𝐱)=0`$. As may be expected, the polyhedron of intersection of the two polyhedra – one polyhedron bounded by the dashed lines corresponding to the method described here and another bounded by the dotted lines corresponding to the method of  – gives better results of over-approximation.

## 6 Discussion and Conclusion

An important issue of hybrid systems appears to be computability of the reach sets of the continuous variables. From a computational point, both the problems of computation and efficient representation of the reach sets of the continuous varibles are difficult, in general, owing to the limitations of quantifier elimination method. In this context, approximation of the reach sets by more convenient sets, such as polyhedra and subalgebraic sets, is discussed in the literature. In this paper, along with a method for finding the reach sets, an algorithm for over-approximation of the reach sets with polyhedra when dynamics of the continuous variables are specified by linear differential equations and the inital set is a polyhedron. A practical version of the over-approximation algorithm is also discussed in this paper. However, it seems that better results of over-approximation may be obtained by taking the intersection of the polyhedron obtained by the method reported here with that obtained by the method given in . It is hoped that the over-approximation method presented here may be extended to systems with more general dynamics and initial sets.

## Appendix A. Decidable Classes of Subsets of $`Q\times ^n`$

We introduce here model theoretic concepts for formalizing the notion of decidable class $`𝒞`$ of subsets of $`Q\times ^n`$ that we were interested in Sec. LABEL:ReachProblems. Since $`Q`$ is a finite set, for a decidable class $`𝒞`$ of subsets of $`Q\times ^n`$, we may take, for each $`qQ`$, a decidable class of subsets of $`^n`$, $`𝒮_q`$, $`qQ`$, and choose $`𝒞`$ to be the product class $`_{qQ}\{\{q\}\times S_q:S_q𝒮_q\}`$. Therefore, we restrict our attention to the discussion of classes of subsets of $`^n`$. Most parts of this section are taken from  (see also \[Laff00\], which in turn cites the references ).

###### Definition 11

A language is a tuple of three sets, $`=(L_r,L_f,L_c)`$, where

1. $`L_r`$ is a set of relations,
2. $`L_f`$ is a set of functions, and
3. $`L_c`$ are a set of constants.

###### Definition 12

A model of a language $`=(L_r,L_f,L_c)`$ consists of a nonempty set $`𝒮`$, together with an interpretation of the relations, functions and constants.

We denote a model by $`(𝒮,L_r,L_f,L_c)`$, where the interpretation is not made explicit. In the following, let $`𝒱=\{x,y,z,x_0,x_2,\mathrm{}\}`$ denote a countable set of variables.

###### Definition 13

A term of a language, $`=(L_r,L_f,L_c)`$, is inductively defined as follows:

1. each variable $`\theta 𝒱`$ is a term,
2. each constant $`cL_c`$ is a term, and
3. for an $`m`$-ary function $`gL_f`$, where $`m1`$, and $`m`$ terms, $`\theta _1,\theta _2,\mathrm{},\theta _m`$, $`g(\theta _1,\theta _2,\mathrm{},\theta _m)`$ is a term.

###### Definition 14

An atomic formula of a language $``$ is either

1. $`\theta _1=\theta _2`$, where $`\theta _1`$ are $`\theta _2`$ two terms of $``$, or
2. $`\rho (\theta _1,\mathrm{},\theta _n)`$, where $`n1`$ and $`\rho L_r`$ is an $`n`$-ary relation.

###### Definition 15

A first order formula, or simply a formula, of a language $``$ is recursively defined as one of the following:

1. an atomic formula, or
2. $`\neg \varphi `$, where $`\varphi `$ is formula and $`\neg `$ is logical negation, or
3. $`\varphi _1\varphi _2`$, where $`\varphi _1`$ and $`\varphi _2`$ are formulas and $``$ is logical and, or
4. $`x:\varphi `$ or $`x:\varphi `$, where $`\varphi `$ is a formula, $`x`$ is a variable, and $``$ (there exists) and $``$ (for all) are quantifiers; in this case, each occurrence of the variable $`x`$ in the formula $`\varphi `$ is called a bound occurrence.

###### Definition 16

The occurrence of a variable in a formula is free, if it is not bound. A sentence in a language $``$ is a formula with no free variable. A theory of $``$ is a subset of sentences.

For a model $`𝒮`$ of the language $``$, we shall be particularly interested in the theory defined by the set of all sentences that are true in $`𝒮`$. To emphasize this, we refer to this theory as the theory of $`(𝒮,L_r,L_f,L_c)`$.

###### Definition 17

Let $``$ be a language and $`𝒮`$ be a model of $``$. A set $`X𝒮^n`$ is said to be definable in the language $``$, if there is an $`n`$-ary formula $`\varphi (x_1,\mathrm{},x_n)`$ such that $`X`$ can be written as $`X=\{(x_1,\mathrm{},x_n)S^n:\varphi (x_1,\mathrm{},x_n)\}`$.

###### Definition 18

Let $``$ be a language, $`𝒮`$ be a model of $``$, and $`𝒞`$ be the class of definable sets. Then $`𝒞`$ is said to be decidable, if the theory of $`(𝒮,L_r,L_f,L_c)`$ is decidable, i.e., there is a decision procedure that, given an $``$-sentence $`\varphi `$, decides whether $`\varphi `$ belongs to the theory of $`(𝒮,L_r,L_f,L_c)`$ or not.

Examples

1. The theory $`(,\{<\},\{+,\},\{0,1\})`$ is the theory of linear constraints with integer coefficients, denoted by $`Lin()`$. The sets defined by these formulas are called polyhedral sets.
2. The theory $`(,\{<\},\{+,,\},\{0,1\})`$ is the theory of polynomial constraints with integer coefficients, denoted by $`OF()`$. The sets defined by these formulas are called subalgebraic sets.

The result stated below is due to A. Tarski :

###### Theorem 1

The first order theory $`OF()`$ is decidable.

###### Definition 19

Let $``$ be a language and let $`𝒮`$ be a model of $``$. We say the theory of $`(𝒮,L_r,L_f,L_c)`$ admits quantifier elimination if every first order formula of $`(𝒮,L_r,L_f,L_c)`$ is equivalent to a formula of $`(𝒮,L_r,L_f,L_c)`$ without quantifiers.

Examples: Decidability and Quantifier Elimination

1. The theory $`OF()`$ consisting of $`(,\{<\},\{+,,\},\{0,1\})`$ admits quantifier elimination and is also decidable.
2. Let $`OF_{exp}()`$ be the theory consisting of $`(,\{<\},\{+,,,exp\},\{0,1\})`$, where $`exp`$, representing the exponential function, is a new function symbol. This theory does not admit quantifier elimination, and it is not known whether this theory is decidable. (See ).

## Appendix B: $`Reach_q(X_0,[0,\tau ])=X(\tau )=X^+(\tau )`$

We assume that $`X_0`$ is closed and its boundary, denoted by $`S_0`$, consists of a finite union of smooth surfaces. Let

$$X(\tau )=X_0\underset{0t\tau }{}\{\varphi _q(𝐱_0,t):𝐱_0S_0\}.$$

Also let

$$S_0^+=\{𝐱S_0:ϵ=ϵ(𝐱)>0\text{ such that }\varphi (𝐱,t)X_0^c,t(0,ϵ)\},$$
(16)
and $`X^+(\tau )=X_0_{0t\tau }\{\varphi _q(𝐱_0,t):𝐱_0S_0^+\}`$.

###### Proposition 1

$`X^+(\tau )=X(\tau )=Reach(X_0,[0,\tau ])`$.

Proof. Note that $`X^+(\tau )X(\tau )Reach(X_0,[0,\tau ])`$. So we have to show that $`Reach(X_0,[0,\tau ])X^+(\tau )`$. Let $`𝐳Reach(X_0,t)`$, $`t>0`$. If $`𝐳X_0`$, then $`𝐳X^+(\tau )`$. Now suppose $`𝐳X_0`$. So there is a $`𝐳_0X_0`$ such that $`𝐳=\varphi (𝐳_0,s)`$, for some $`s`$ with $`0<s\tau `$. If $`𝐳_0S_0^+`$, then $`𝐳X^+(\tau )`$. Otherwise let $`\tau ^{}=sup\{t:ts\text{ and }\varphi (𝐳_0,t)X_0\}`$, so $`0<\tau ^{}<s\tau `$. Now $`𝐳^{}=\varphi (𝐳_0,\tau ^{})S_0^+`$. Therefore, $`𝐳=\varphi (𝐳^{},s\tau ^{})X^+(\tau )`$.

Now suppose that $`X_0`$ is specified as $`X_0=\{𝐱^n:\mathrm{}(𝐱)0\}`$, where $`\mathrm{}:^n`$ is continuously differentiable. Further, we assume that if $`𝐱\stackrel{}{X_0}`$ then $`\mathrm{}(𝐱)<0`$ (where $`\stackrel{}{X_0}`$ denotes the interior of $`X_0`$, i.e., the largest open set contained in $`X_0`$) So, obviously, if $`𝐱S_0`$ then $`\mathrm{}(𝐱)=0`$, and if $`𝐱X_0^c`$ then $`\mathrm{}(𝐱)>0`$. Let $`S_1^+S_0`$ be defined as

$$S_1^+=\{𝐱S_0:\mathrm{}(𝐱)f(q,𝐱)0\}.$$

We now show that $`S_0^+S_1^+`$. To this end, we show that $`S_0\backslash S_1^+S_0\backslash S_0^+`$. Let $`𝐱S_0\backslash S_1^+`$, so $`\mathrm{}(𝐱)f(q,𝐱)<0`$. Now, since $`\mathrm{}(\varphi _q(𝐱,0))=\mathrm{}(𝐱)=0`$, and since at $`t=0`$, $`\frac{d\mathrm{}(\varphi _q(𝐱,t))}{dt}=\mathrm{}(𝐱)f(q,𝐱)<0`$, we have, for a sufficiently small $`ϵ>0`$, $`\mathrm{}(\varphi _q(𝐱,t))<0`$, whenever $`0<t<ϵ`$, which happens only if $`\varphi _q(𝐱,t)\stackrel{}{X_0}`$, $`t(0,ϵ)`$, implying that $`𝐱S_0^+`$; hence $`𝐱S_0\backslash S_0^+`$. Thus $`S_0\backslash S_1^+S_0\backslash S_0^+`$, as required.

## Appendix C: $`Reach_q^{}(X_0,X_q,[0,\mathrm{}))`$

Let $`X_q^n`$ with compact closure, and let $`Y_q=\overline{X_q}`$. Let $`𝐟:W^n`$ be a continuous function defined on an open set $`W`$ containing $`Y_q`$, satisfying a Lipschitz condition on $`W`$. Let $`X_0`$ be a closed subset of $`Y_q`$. We assume that for each $`𝐱X_0`$, a function $`\gamma _𝐱:^+^n`$ exists and satisfies the differential equation

$$\frac{d\gamma _𝐱(t)}{dt}=𝐟(\gamma _𝐱),t0,$$
(17)
with the initial condition $`\gamma _𝐱(0)=𝐱`$. Hence the flow $`\varphi (𝐱,t)`$ associated with equation (17) is defined for all $`t0`$ and $`𝐱X_0`$. Further, assume that the $`\omega `$-limit set of the flow $`\varphi (𝐱,t)`$, $`t0`$, does not intersect $`Y_q`$ for any point $`𝐱X_0`$. More precisely, we assume, for $`𝐱X_0`$,

$$L_\omega (𝐱)Y_q=\mathrm{},$$
(18)
where

$$L_\omega (𝐱)=\underset{t0}{}\overline{\varphi (𝐱,[t,\mathrm{}))}.$$

(See .) In what follows, we show that if (18) holds for every $`𝐱X_0`$, then $`\tau _{\mathrm{max}}>0`$ (depending on $`X_0`$), such that for every $`𝐱X_0`$, there is a $`\tau =\tau (𝐱)`$ with $`0<\tau (𝐱)\tau _{max}`$ and $`\varphi (𝐱,t(𝐱))Y_q`$. For this, let $`0<\tau _1<\tau _2<\mathrm{}`$ be an increasing sequence such that $`\tau _k\mathrm{}`$, as $`k\mathrm{}`$ (as in Sec. 3), and define the sets $`A(\tau _k,𝐱)=\overline{\varphi (𝐱,[\tau _k,\mathrm{}))}`$. We first observe that $`L_\omega (𝐱)=_kA(\tau _k,𝐱)`$.

###### Proposition 2

Let $`𝐱X_0`$. If (18) holds, then $`\tau =\tau (𝐱)>0`$ such that $`\varphi (𝐱,t)Y_q`$, $`t\tau `$.

Proof. We have $`L_\omega (𝐱)Y_q=`$ $`\left[_kA(\tau _k,𝐱)\right]Y_q`$ $`=_k\left[A(\tau _k,𝐱)Y_q\right]`$. Now let $`E(\tau _k,𝐱)=A(\tau _k,𝐱)Y_q`$. So, for a fixed $`𝐱X_0`$, the sets $`E(\tau _k,𝐱)`$, $`k=1,2,3,\mathrm{}`$, is a decreasing sequence of closed subsets of $`Y_q`$. Since $`Y_q`$ is compact, $`_kE(\tau _k,𝐱)=\mathrm{}`$ implies $`E(\tau _k,𝐱)=\mathrm{}`$, for all but finitely many $`k`$. Therefore, for some $`K1`$, $`kK`$, $`E(\tau _k,𝐱)=\mathrm{}`$. Hence, $`\varphi (𝐱,[t,\mathrm{}))Y_q=\mathrm{}`$, $`t\tau _K`$. Therefore, with $`\tau =\tau _K`$, $`\varphi (𝐱,t)Y_q`$, $`t\tau `$.

We also need the following proposition:

###### Proposition 3

Let $`𝐱X_0`$ be a point for which there is a $`\tau >0`$, such that $`\varphi (𝐱,\tau )Y_q`$. Then there is a $`\delta =\delta _𝐱>0`$, such that if $`|𝐱𝐲|<\delta _𝐱`$ and $`𝐲X_0`$, then $`\varphi (𝐲,\tau )Y_q`$.

Proof. Let $`|𝐟(𝐱)𝐟(𝐲)|C|𝐱𝐲|`$. Then $`|\varphi (𝐱,\tau )\varphi (𝐲,\tau )||𝐱𝐲|e^{C\tau }`$ (see, for example, ), so for a fixed $`\tau `$, $`\varphi (,\tau )`$ is continuous in the first variable. Let $`𝐳=\varphi (𝐱,\tau )`$. Now since $`W\backslash Y_q`$ is open, there is an $`ϵ>0`$ such that $`B(𝐳,ϵ)W\backslash Y_q`$. By the continuity of $`\varphi (,\tau )`$ at $`𝐱`$, there is a $`\delta >0`$, such that $`|\varphi (𝐱,\tau )\varphi (𝐲,\tau )|=|𝐳\varphi (𝐲,\tau )|<ϵ`$, whenever $`|𝐲𝐱|<\delta `$ and $`𝐲X_0`$. Therefore, $`\varphi (𝐲,\tau )B(𝐳,ϵ)W\backslash Y_q`$.

From the previous two propositions, we get the following theorem:

###### Theorem 2

Let $`Y_q`$ be compact, and $`𝐟`$ be a continuous function defined on an open set $`W`$ containing $`Y_q`$ satisfying a Lipschitz condition. Further assume condition (18) holds for every $`𝐱X_0`$. Then $`\tau _{\mathrm{max}}>0`$ independent of $`𝐱`$, such that for each $`𝐱X_0`$, there is a $`\tau (𝐱)`$ with $`0<\tau (𝐱)\tau _{max}`$ and $`\varphi (𝐱,\tau (𝐱))Y_q`$.

Proof. Let $`𝐱X_0`$. By Proposition 2, there is a $`\tau (𝐱)>0`$, such that $`\varphi (𝐱,\tau (𝐱))Y_q`$. By Proposition 3, there is a $`\delta _𝐱>0`$, such that for any $`𝐲X_0`$ with $`|𝐲𝐱|<\delta _𝐱`$, $`\varphi (𝐲,\tau (𝐱))Y_q`$. $`\{B(𝐱,\delta _𝐱):𝐱X_0\}`$ is an open cover of $`X_0`$, containing a finite subcover, say, $`\{B(𝐱_1,\delta _{𝐱_1}),\mathrm{},B(𝐱_m,\delta _{𝐱_m})\}`$. Let $`\tau _{\mathrm{max}}=\mathrm{max}\{\tau (𝐱_1),\mathrm{},\tau (𝐱_m)\}`$. To check wether this choice of $`\tau _{max}`$ is as in the theorem, let $`𝐲X_0`$ be an arbitrary point. Now $`𝐲B(𝐱_i,\delta _{𝐱_i})`$, for some $`i`$ with $`1im`$, and by the choice of $`\delta _{𝐱_i}`$, $`\varphi (𝐲,\tau (𝐱_i))Y_q`$, concluding the proof.

## Appendix D: Boundedness Results

In this section, we show that if the initial set $`F`$ is bounded then the polyhedron, $`P`$, enclosing $`F_{[0,\mathrm{\Delta }]}`$, as obtained by the algorithm described in Sec. 3 is bounded. We assume that $`F`$ is nonempty. Recall that $`F`$ is given as the set of points $`𝐱`$ in $`^n`$ which satisfy the following constraints:

$$\begin{array}{ccc}\hfill 𝐚_i^T𝐱b_i& & 0,i=1,2,\mathrm{},k1,\hfill \\ \hfill 𝐚_k^T𝐱b_k& =& 0,\hfill \end{array}\}$$
(19)
and $`P`$ is included in the set of points $`𝐱^n`$ satisfying

$$\begin{array}{ccc}\hfill 𝐚_i^T𝐱b_il_i(𝐚_k^T𝐱b_k)& & 0,i=1,2,\mathrm{},k1,\hfill \\ \hfill 𝐚_k^T𝐱b_k& & 0,\text{ and}\hfill \\ \hfill 𝐚_k^T𝐱b_kl_k& & 0,\hfill \end{array}\}$$
(20)
where $`l_i`$ and $`l_k>0`$. Looking at the constraints, one may visualize the set $`P_1`$ satisfying (20) as a prism or a truncated pyramid, with its bottom given by (19) and its top given by

$$\begin{array}{ccc}\hfill 𝐚_i^T𝐱b_il_i(𝐚_k^T𝐱b_k)& & 0,i=1,2,\mathrm{},k1,\hfill \\ \hfill 𝐚_k^T𝐱b_kl_k& =& 0.\hfill \end{array}\}$$
(21)
We assume that $`F`$ is nonempty, and wish to show that if $`F`$ is bounded then so is the set $`P_1`$, defined as the set of points satisfying (20). We start with the following proposition:

###### Proposition 4

Assume that $`F`$ is nonempty and $`P_1`$ is not bounded. Then there is a point $`𝐱_0F`$ and a vector $`\lambda ^n`$, with $`|\lambda |=1`$, such that $`𝐱_0+t\lambda P_1`$, for all $`t0`$.

Proof. Fix a point $`𝐱_0F`$, and suppose $`P_1`$ is not bounded. So for each positive integer $`k`$, there is a point $`𝐱_kP_1`$, such that $`|𝐱_k𝐱_0|k`$. Let $`\lambda _k=\frac{(𝐱_k𝐱_0)}{|𝐱_k𝐱_0|}`$. Now, since $`FP_1`$, $`𝐱_0P_1`$ and, by convexity, the entire line segment $`𝐱_0+t\lambda _kP_1`$, for $`0tk`$. In particular, the points on these line segments satisfy (20). Now, since for each $`k`$, $`|\lambda _k|=1`$, and the closed unit ball in $`^n`$ is compact, there is a convergent subsequence of $`\lambda _k`$’s – say, $`\lambda _{k_i}`$, $`i=1,2,3,\mathrm{}`$ – such that $`\lambda _{k_i}\lambda ^n`$, as $`i\mathrm{}`$; since $`|\lambda _{k_i}|=1`$, $`|\lambda |=1`$. We show that $`𝐱_0+t\lambda P_1`$ for $`t[0,\mathrm{})`$. First note that for any positive integer $`m`$, and for all $`im`$, $`𝐱_0+t\lambda _{k_i}`$ satisfies the constraints (20), for $`t[0,k_m]`$. Therefore, $`𝐱_0+t\lambda `$ satisfies (20), for $`t[0,k_m]`$, and so $`𝐱_0+t\lambda P_1`$, for $`t[0,k_m]`$. The proposition is concluded by letting $`m\mathrm{}`$.

We now show that such an $`\lambda `$, as in the previous proposition, must be parallel to the hyperplane passing through $`F`$, the normal of which is given by $`𝐚_k`$.

###### Proposition 5

If $`𝐱_0F`$ and a unit vector $`\lambda ^n`$ are such that $`𝐱_0+t\lambda P_1`$, for all $`t[0,\mathrm{})`$, then $`𝐚_k^T\lambda =0`$.

Proof. For a contradiction assume that $`h=𝐚_k^T\lambda >0`$. Now since $`𝐱+t\lambda `$ satisfies (20), we must have

$$𝐚_k^T(𝐱_0+t\lambda )b_kl_k0,$$

which holds only if $`t\frac{(b_k+l_k𝐚_k^T𝐱_0)}{h}`$, contrary to the hypothesis that, for all $`t[0,\mathrm{})`$, $`𝐱_0+t\lambda P_1`$. Similarly, if $`h=𝐚_k^T\lambda <0`$, then the constraint

$$𝐚_k^T(𝐱_0+t\lambda )b_k0,$$

does not hold for $`t>0`$. Therefore we must have $`𝐚_k^T\lambda =0`$.

We now show that if $`F`$ is bounded, then for any unit vector $`\lambda `$ orthogonal to $`𝐚_k`$, there is an $`i`$ with $`1ik1`$, such that $`h_i=𝐚_i^T\lambda >0`$.

###### Proposition 6

If $`F`$ is nonempty and bounded, then for any unit vector $`\lambda ^n`$, for which $`𝐚_k^T\lambda =0`$, there is an $`i`$, with $`1ik1`$, such that $`𝐚_i^T\lambda >0`$.

Proof. Let $`𝐱_0F`$, and suppose that there is an $`\lambda ^n`$, with $`|\lambda |=1`$, $`𝐚_k^T\lambda =0`$ and $`𝐚_i^T\lambda 0`$, for $`1ik1`$. Then for any $`t>0`$, the vector $`𝐱_0+t\lambda `$ satisfies the constraints (19), and so $`𝐱_0+t\lambda F`$, contrary to the assumption that $`F`$ is bounded.

We now combine all the previous propositions to get the following result:

###### Theorem 3

If $`F`$ is nonempty and bounded, then $`P_1`$ is bounded.

Proof. If $`P_1`$ is not bounded, then by Proposition 4, there is an $`𝐱_0F`$ and a unit vector $`\lambda ^n`$ such that $`𝐱_0+t\lambda P_1`$, for any $`t0`$. By Proposition 5, $`𝐚_k^T\lambda =0`$. Now, by the previous proposition, since $`F`$ is bounded, there is an $`i`$, with $`1ik1`$, such that $`h_i=𝐚_i^T\lambda >0`$. But then the constraint

$$𝐚_i^T(𝐱_0+t\lambda )b_il_i(𝐚_k^T(𝐱_0+t\lambda )b_k)0,$$

cannot hold for any $`t`$ with

$$t>\frac{𝐚_i^T𝐱_0b_i}{h_i},$$

contrary to the assumption that $`𝐱_0+t\lambda P_1`$, for any $`t0`$. Hence $`P_1`$ is bounded.