# Counting regions of the boxed threshold arrangement

Priyavrat Deshpande  Chennai Mathematical Institute<EMAIL_ADDRESS>,
Krishna Menon  Chennai Mathematical Institute<EMAIL_ADDRESS>and
Anurag Singh  Chennai Mathematical Institute<EMAIL_ADDRESS>
###### Abstract.

In this paper we consider the hyperplane arrangement in $\mathbb{R}^{n}$ whose
hyperplanes are $\\{x_{i}+x_{j}=1\mid 1\leq i<j\leq n\\}\cup\\{x_{i}=0,1\mid
1\leq i\leq n\\}$. We call it the _boxed threshold arrangement_ since we show
that the bounded regions of this arrangement are contained in an $n$-cube and
are in one-to-one correspondence with the labeled threshold graphs on $n$
vertices. The problem of counting regions of this arrangement was studied
earlier by Joungmin Song. He determined the characteristic polynomial of this
arrangement by relating its coefficients to the count of certain graphs. Here,
we provide bijective arguments to determine the number of regions. In
particular, we construct certain signed partitions of the set
$\\{-n,\dots,n\\}\setminus\\{0\\}$ and also construct colored threshold graphs
on $n$ vertices and show that both these objects are in bijection with the
regions of the boxed threshold arrangement. We independently count these
objects and provide a closed form formula for the number of regions.

###### Key words and phrases:

threshold graph, hyperplane arrangement, finite field method.

###### 2010 Mathematics Subject Classification:

52C35, 32S22, 05C30, 11B68

PD and AS are partially supported by a grant from the Infosys Foundation

##  1\. Introduction

A _hyperplane arrangement_ $\mathcal{A}$ is a finite collection of affine
hyperplanes (i.e., codimension $1$ subspaces and their translates) in
$\mathbb{R}^{n}$. Without loss of generality we assume that arrangements we
consider are _essential_ , i.e., the subspace spanned by the normal vectors is
the ambient vector space. Removing all the hyperplanes of $\mathcal{A}$ leaves
$\mathbb{R}^{n}$ disconnected; counting the number of connected components
using diverse combinatorial methods is an active area of research. A _flat_ of
$\mathcal{A}$ is a nonempty intersection of some of the hyperplanes in
$\mathcal{A}$; the ambient vector space is a flat since it is an intersection
of no hyperplanes. Flats are naturally ordered by reverse set inclusion; the
resulting poset is called the _intersection poset_ and is denoted by
$\mathrm{L}(\mathcal{A})$. A _region_ of $\mathcal{A}$ is a connected
component of $\mathbb{R}^{n}\setminus\bigcup\mathcal{A}$. The _characteristic
polynomial_ of $\mathcal{A}$ is defined as

$\chi_{\mathcal{A}}(t):=\sum_{x\in\mathrm{L}(\mathcal{A})}\mu(\hat{0},x)\,t^{\dim(x)}$

where $\mu$ is the Möbius function of the intersection poset and $\hat{0}$
corresponds to the flat $\mathbb{R}^{n}$. The characteristic polynomial is a
fundamental combinatorial and topological invariant of the arrangement and
plays a significant role throughout the theory of hyperplane arrangements.
Among other things, the polynomial encodes enumerative information about the
stratification of the space $\mathbb{R}^{n}$, induced by the arrangement. We
refer the reader to Stanley’s notes [9] for more on the enumerative aspects of
hyperplane arrangements. In particular we have the following seminal result by
Zaslavsky.

######  Theorem 1.1 ([10]).

Let $\mathcal{A}$ be an arrangement in $\mathbb{R}^{n}$. Then the number of
regions of $\mathcal{A}$ is given by

$r(\mathcal{A})=(-1)^{n}\chi_{\mathcal{A}}(-1)$

and the number of bounded regions is given by

$b(\mathcal{A})=(-1)^{n}\chi_{\mathcal{A}}(1).$

The finite field method, developed by Athanasiadis [1], converts the
computation of the characteristic polynomial to a point counting problem. A
combination of these two results allowed for the computation of the number of
regions of several arrangements of interest.

Another way to count the number of regions is to give a bijective proof. This
approach involves finding a combinatorially defined set whose elements are in
bijection with the regions of the given arrangement and are easier to count.
For example, the regions of the braid arrangement (whose hyperplanes are given
by the equations $x_{i}-x_{j}=0$ for $1\leq i<j\leq n$) correspond to the $n!$
permutations of order $n$. The regions of the Shi arrangement (whose
hyperplanes are given by the equations $x_{i}-x_{j}=0,1$ for $1\leq i<j\leq
n$) are in bijection with the parking functions on $[n]$, hence the number of
regions is $(n+1)^{n-1}$. Refer to Stanley’s notes [9, Lecture 5] for details.

In the present article we consider the following hyperplane arrangement

$\mathcal{BT}_{n}:=\\{X_{i}+X_{j}=1\mid 1\leq i<j\leq n\\}\cup\\{X_{i}=0,1\mid
1\leq i\leq n\\}.$

The problem of counting regions of this arrangement was recently solved by
Song in a series of papers. His first approach [6, 5] involved relating the
coefficients of the characteristic polynomial to the generating functions for
the number of certain graphs. These generating functions [5, Theorem 1] have
quite a complicated expression and consequently make it difficult to determine
a compact form for the characteristic polynomial. Later Song [7] used the
finite field method to get a slightly simpler expression for the
characteristic polynomial.

The main aim of this paper is to provide bijective proofs for the number of
regions of this arrangement. In particular, in Section 3 we construct certain
(signed) ordered partitions of the set
$\\{-n,-(n-1),\dots,n-1,n\\}\setminus\\{0\\}$ and show that they are in
bijection with the regions of $\mathcal{BT}_{n}$. We also show how to count
these partitions. In Section 4 we spell out a recipe to color the vertices of
a labeled threshold graph on $n$ vertices such that the number of such colored
threshold graphs is equal to the regions of $\mathcal{BT}_{n}$. However, we
begin the article by establishing a simpler form for the characteristic
polynomial in Section 2. There we also justify the term “ _boxed threshold_ ”.

##  2\. The characteristic polynomial

We first translate the hyperplanes in $\mathcal{BT}_{n}$ in order to obtain a
combinatorially isomorphic arrangement with the same characteristic
polynomial. Putting $X_{i}=x_{i}+\frac{1}{2}$ for every $i$ we get

(1) $\\{x_{i}+x_{j}=0\mid 1\leq i<j\leq
n\\}\cup\\{x_{i}=-\frac{1}{2},\frac{1}{2}\mid 1\leq i\leq n\\}.$

We will stick to the notation $\mathcal{BT}_{n}$ to denote the above
arrangement. The notation $[k]$ is used to denote the set $\\{1,\dots,k\\}$.
We begin by stating a generalization of the finite field method, given by
Athanasiadis [2], in our context.

######  Theorem 2.1 ([2, Theorem 2.1]).

If $\mathcal{A}$ is a sub-arrangement of the type $C$ arrangement in
$\mathbb{R}^{n}$, that is, a sub-arrangement of $\\{x_{i}\pm x_{j}=0\mid 1\leq
i<j\leq n\\}\cup\\{x_{i}=0\mid i\in[n]\\}$, there exists an integer $k$ such
that for all odd integers $q$ greater than $k$,

$\chi_{\mathcal{A}}(q)=\\#(\mathbb{Z}_{q}^{n}\setminus V_{\mathcal{A}})$

where $V_{\mathcal{A}}$ is the union of hyperplanes in $\mathbb{Z}_{q}^{n}$
obtained by reducing $\mathcal{A}$ mod q.

We use this result to prove the following relationship between the
characteristic polynomial of sub-arrangements of the type $C$ arrangement and
that of their “boxed” versions.

######  Proposition 2.2.

Let $\mathcal{A}$ be an arrangement in $\mathbb{R}^{n}$ that is a sub-
arrangement of the type $C$ arrangement and let
$\mathcal{BA}=\mathcal{A}\cup\\{x_{i}=-\frac{1}{2},\frac{1}{2}\mid
i\in[n]\\}$. Then

$\chi_{\mathcal{BA}}(t)=\chi_{\mathcal{A}}(t-2).$

###### Proof.

Let $q$ be any large odd number. Set
$D_{q}^{n}:=\\{(a_{1},\dots,a_{n})\in\mathbb{Z}_{q}^{n}\mid
a_{i}\neq\pm\frac{q-1}{2}\\}$. Define a bijection
$f:\mathbb{Z}_{q-2}\rightarrow\mathbb{Z}_{q}\setminus\\{\frac{q-1}{2},-\frac{q-1}{2}\\}$
as

$f(i)=i\quad\text{for }i\in[-\frac{q-3}{2},\frac{q-3}{2}].$

It is clear that for any $a,b\in\mathbb{Z}_{q-2}$

  1. (1)

$a+b=0\Leftrightarrow f(a)+f(b)=0$,

  2. (2)

$a-b=0\Leftrightarrow f(a)-f(b)=0$, and

  3. (3)

$a=0\Leftrightarrow f(a)=0$.

Using $f$, we can define a bijection $F:\mathbb{Z}_{q-2}^{n}\rightarrow
D_{q}^{n}$ as

$F(a_{1},\dots,a_{n})=(f(a_{1}),\dots,f(a_{n}))\quad\text{for
}(a_{1},\dots,a_{n})\in\mathbb{Z}_{q-2}^{n}.$

By the properties of $f$, we can see that $F$ induces a bijection between
those tuples in $\mathbb{Z}_{q-2}^{n}$ that do not satisfy the defining
equation of any hyperplane in $\mathcal{A}$ and those tuples in
$\mathbb{Z}_{q}^{n}$ that do not satisfy the defining equation of any
hyperplane in $\mathcal{BA}$. So, we get that for large odd numbers $q$,

$\chi_{\mathcal{BA}}(q)=\chi_{\mathcal{A}}(q-2).$

Since $\chi_{\mathcal{BA}}$ and $\chi_{\mathcal{A}}$ are polynomials, we get
the required result. ∎

Let $\mathcal{T}_{n}$ denote the threshold arrangement in $\mathbb{R}^{n}$,
i.e., $\mathcal{T}_{n}:=\\{x_{i}+x_{j}=0\mid 1\leq i<j\leq n\\}$. The reason
this is called the threshold arrangement is that its regions are in bijection
with labeled threshold graphs on $n$ vertices (see Section 4 for details).
This is clearly a sub-arrangement of the type $C$ arrangement.

######  Corollary 2.3.

The characteristic polynomials of $\mathcal{BT}_{n}$ and $\mathcal{T}_{n}$ are
related as follows:

$\chi_{\mathcal{BT}_{n}}(t)=\chi_{\mathcal{T}_{n}}(t-2).$

Consequently, the number of bounded regions of $\mathcal{BT}_{n}$ is equal to
the number of regions of $\mathcal{T}_{n}$. Moreover, these bounded regions
are contained in the cube (or a _box_)
$\left[-\frac{1}{2},\frac{1}{2}\right]^{n}$. Next, we derive a closed form
expression for $\chi_{\mathcal{T}_{n}}(t)$ using the finite field method.

######  Proposition 2.4.

The characteristic polynomial of the threshold arrangement $\mathcal{T}_{n}$
is

$\chi_{\mathcal{T}_{n}}(t)=\sum_{k=1}^{n}(S(n,k)+n\cdot
S(n-1,k))\prod_{i=1}^{k}(t-(2i-1)).$

Here $S(n,k)$ are the Stirling numbers of the second kind.

###### Proof.

Using Theorem 2.1, we see that the characteristic polynomial of
$\mathcal{T}_{n}$ satisfies, for large odd values of $q$,

$\chi_{\mathcal{T}_{n}}(q)=|\\{(a_{1},\dots,a_{n})\in\mathbb{Z}_{q}^{n}\mid
a_{i}+a_{j}\neq 0\text{ for all }1\leq i<j\leq n\\}|.$

This means that we need to count the functions
$f:[n]\rightarrow\mathbb{Z}_{q}$ such that

  1. (1)

$|f^{-1}(0)|\leq 1$.

  2. (2)

$f$ can take at most one value from each of the sets

$\\{1,-1\\},\\{2,-2\\},\dots,\\{\frac{q-1}{2},-\frac{q-1}{2}\\}.$

We split the count into the two cases. If 0 is not attained by $f$, then all
values must be from

$\\{1,-1\\}\cup\\{2,-2\\}\cup\cdots\cup\\{\frac{q-1}{2},-\frac{q-1}{2}\\}.$

with at most one value attained in each set. So, there are

$\binom{\frac{q-1}{2}}{k}\cdot 2^{k}\cdot k!\cdot S(n,k)$

ways for $f$ to attain values from exactly $k$ of these sets. This is because
we have $\binom{\frac{q-1}{2}}{k}\cdot 2^{k}$ ways to choose the $k$ sets and
which element of each set $f$ should attain and $k!S(n,k)$ ways to choose the
images of the elements of $[n]$ after making this choice. So the total number
of $f$ such that 0 is not attained is

$\sum\limits_{k=1}^{n}\binom{\frac{q-1}{2}}{k}\cdot 2^{k}\cdot k!\cdot
S(n,k).$

When 0 is attained, there are $n$ ways to choose which element of $[n]$ gets
mapped to 0 and using a similar logic for choosing the images of the other
elements, we get that the total number of $f$ where 0 is attained is

$n\cdot\sum\limits_{k=1}^{n-1}\binom{\frac{q-1}{2}}{k}\cdot 2^{k}\cdot k!\cdot
S(n-1,k).$

So we get that for large $q$,

$\displaystyle\chi_{\mathcal{T}_{n}}(q)$
$\displaystyle=\sum\limits_{k=1}^{n}\binom{\frac{q-1}{2}}{k}\cdot 2^{k}\cdot
k!\cdot S(n,k)+n\cdot\sum\limits_{k=1}^{n-1}\binom{\frac{q-1}{2}}{k}\cdot
2^{k}\cdot k!\cdot S(n-1,k)$ $\displaystyle=\sum_{k=1}^{n}(S(n,k)+n\cdot
S(n-1,k))\prod_{i=1}^{k}(q-(2i-1)).$

Since $\chi_{\mathcal{T}_{n}}$ is a polynomial, we get the required result. ∎

######  Remark 2.5.

Note that the absolute value of the coefficient of $t^{j}$ in
$(t-1)(t-3)\cdots(t-(2k-1))$ counts the number of signed permutations on $[k]$
with $j$ odd cycles (see A039757 in the OEIS [4]). Let us denote that number
by $a(k,j)$. Using this we get a compact expression for the coefficient of
$t^{j}$ in $\chi_{\mathcal{T}_{n}}(t)$ as

$\sum_{k=j}^{n}(S(n,k)+n\cdot S(n-1,k))a(k,j).$

######  Corollary 2.6.

The coefficient of $t^{j}$ in $\chi_{\mathcal{BT}_{n}}(t)$ is

$\sum_{k=j}^{n}(S(n,k)+n\cdot S(n-1,k))b(k,j).$

where $b(k,j)$ is the coefficient of $t^{j}$ in $(t-3)(t-5)\cdots(t-(2k+1))$.

###### Proof.

Using Corollary 2.3 we get

$\chi_{\mathcal{BT}_{n}}(t)=\sum_{k=1}^{n}(S(n,k)+n\cdot
S(n-1,k))\prod_{i=1}^{k}(t-(2i+1)).$

Expanding the product gives us the required expression. Note that
$b(k,j)=-\sum\limits_{i=0}^{j}a(k+1,i)$ where $a(k,j)$ is defined in Remark
2.5. ∎

######  Remark 2.7.

We can also derive an expression for the exponential generating function for
the characteristic polynomial. Using Problem $25(c)$ in Stanley’s notes [9,
Lecture 5] we get

$\sum_{n\geq
0}\chi_{\mathcal{BT}_{n}}(t)\frac{x^{n}}{n!}=(1+x)(2e^{x}-1)^{\frac{(t-3)}{2}}.$

The generating function for the number of regions is

$\sum_{n\geq
0}r(\mathcal{BT}_{n})\frac{x^{n}}{n!}=\frac{e^{2x}(1-x)}{(2-e^{x})^{2}}.$

For the sake of completeness we enumerate the coefficients of the
characteristic polynomial for smaller values of $n$ (see Table 1). Song [5]
also computed the characteristic polynomial for $n\leq 10$, however there are
typos in all the expressions for $n\geq 4$, consequently the region numbers
are wrong. The sequence of number of regions of $\mathcal{BT}_{n}$ can be
found at the entry A341769 in the OEIS [4].

$n$ |  $\chi_{\mathcal{BT}_{n}}(t)$ |  $r(\mathcal{BT}_{n})$
---|---|---
2 |  $t^{2}-5t+6$ |  12
3 |  $t^{3}-9t^{2}+27t-27$ |  64
4 |  $t^{4}-14t^{3}+75t^{2}-181t+165$ |  436
5 |  $t^{5}-20t^{4}+165t^{3}-695t^{2}+1480t-1263$ |  3624
6 |  $t^{6}-27t^{5}+315t^{4}-2010t^{3}+7320t^{2}-14284t+11559$ |  35516
7 |  $t^{7}-35t^{6}+546t^{5}-4865t^{4}+26460t^{3}-87010t^{2}+158753t-122874$ |  400544
8 |  $t^{8}-44t^{7}+882t^{6}-10402t^{5}+78155t^{4}-379666t^{3}+1154965t^{2}-1995487t+1486578$ |  5106180
9 |  $t^{9}-54t^{8}+1350t^{7}-20286t^{6}+200025t^{5}-1331022t^{4}+5932143t^{3}-16952157t^{2}+27979203t-20158695$ |  72574936
10 |  $t^{10}-65t^{9}+1980t^{8}-36840t^{7}+459585t^{6}-3986031t^{5}+24172575t^{4}-100548090t^{3}+272771475t^{2}-432836011t+302751327$ |  1137563980
Table 1. Characteristic polynomial and the number of regions of
$\mathcal{BT}_{n}$ for $n\leq 10$.

##  3\. The signed ordered partitions

In this section we prove a bijection between regions of $\mathcal{BT}_{n}$,
and certain ordered partitions of
$[-n,n]\setminus\\{0\\}\cup\\{-\frac{1}{2},\frac{1}{2}\\}$ (the notation
$[-n,n]$ is used for the set $\\{-n,-n+1,\dots,n-1,n\\}$). We will denote
$-x_{i}$ by $x_{-i}$ for all $i\in[n]$. Let $R$ be a region of
$\mathcal{BT}_{n}$. We write

  1. (1)

$i\prec_{R}-j$ if $x_{i}+x_{j}<0$ in $R$ where $i\neq j$ in $[n]$.

  2. (2)

$i\succ_{R}-j$ if $x_{i}+x_{j}>0$ in $R$ where $i\neq j$ in $[n]$.

  3. (3)

$i\succ_{R}\frac{1}{2}$ if $x_{i}>\frac{1}{2}$ in $R$ where
$i\in[-n,n]\setminus\\{0\\}$.

  4. (4)

Similarly define $i\prec_{R}\frac{1}{2}$, $i\succ_{R}-\frac{1}{2}$ and
$i\prec_{R}-\frac{1}{2}$ for any $i\in[-n,n]\setminus\\{0\\}$.

So $\prec_{R}$ is a relation on
$[-n,n]\setminus\\{0\\}\cup\\{-\frac{1}{2},\frac{1}{2}\\}$ where comparable
elements are

  1. (1)

Elements of $[-n,n]\setminus\\{0\\}$ of different signs and different absolute
values.

  2. (2)

Elements of $[-n,n]\setminus\\{0\\}$ and $\pm\frac{1}{2}$.

######  Remark 3.1.

The reader can consider the relation $i\prec_{R}-j$ as $\overset{+}{i}$
appearing before $\overset{-}{j}$ in a signed permutation on $[n]$ and
similarly for other such relations. The equivalence defined in the following
lemma corresponds to choosing a signed permutation representative for each
region. This is similar to the way Spiro [8] associated ‘threshold pairs in
standard form’ to labeled threshold graphs.

######  Lemma 3.2.

Let $\sim$ be the relation defined on the set $N=[-n,n]\setminus\\{0\\}$ as

$a\sim b$ if the following hold:

  1. (1)

$a$ and $b$ are of the same sign, and

  2. (2)

there does not exist any $c\in N\cup\\{-\frac{1}{2},\frac{1}{2}\\}$ such that
$a\prec_{R}c\prec_{R}b$ or $b\prec_{R}c\prec_{R}a$.

Then, $\sim$ is an equivalence relation on $N$.

###### Proof.

It is clear that $a\sim a$ for any $a\in N$ and that $a\sim b$ implies that
$b\sim a$ for any $a,b\in N$. Let $a\sim b$ and $b\sim c$ for some distinct
$a,b,c\in N$. So $a,b,c$ are of the same sign. By definition of $\prec_{R}$
and $\sim$, the only possible $d\in N\cup\\{-\frac{1}{2},\frac{1}{2}\\}$ such
that $a\prec_{R}d\prec_{R}c$ is $-b$. We must show that this is not possible
(a similar argument works if $c\prec_{R}-b\prec_{R}a$). Suppose
$a\prec_{R}-b\prec_{R}c$. We then have $-c\prec_{R}b\prec_{R}-a$. If
$a\prec_{R}-c$, we will have $a\prec_{R}-c\prec_{R}b$, which contradicts
$a\sim b$. So we must have $-c\prec_{R}a$. But this gives $-a\prec_{R}c$ and
hence $b\prec_{R}-a\prec_{R}c$, which is a contradiction to $b\sim c$. Hence
$\sim$ is an equivalence on $N$. ∎

######  Definition 3.3.

The equivalence classes of $\sim$ defined in Lemma 3.2 will be called boxed
threshold blocks (corresponding to the region $R$). Positive blocks refer to
those blocks that contain positive numbers and similarly we define negative
blocks.

######  Remark 3.4.

Since all the elements of a boxed threshold block are of the same sign, we
sometimes consider a block as a subset of $[n]$ with a sign ($+$ or $-$)
assigned to it.

The proof of the following lemma follows from the definition of the
equivalence.

######  Lemma 3.5.

If $B$ is a boxed threshold block then so is $-B=\\{-b:b\in B\\}$.

######  Proposition 3.6.

Define an order $<$ on the set of boxed threshold blocks along with
$\pm\frac{1}{2}$ as follows:

  1. (1)

$B<B^{\prime}$ where $B,B^{\prime}$ are boxed threshold blocks if there exists
a sequence $c_{0},c_{1},\dots,c_{k}$ of elements in
$N\cup\\{-\frac{1}{2},\frac{1}{2}\\}$ such that
$c_{0}\prec_{R}c_{1}\prec_{R}\cdots\prec_{R}c_{k}$, $c_{0}\in B$ and $c_{k}\in
B^{\prime}$.

  2. (2)

$-\frac{1}{2}<\frac{1}{2}$.

  3. (3)

$B<\frac{1}{2}$ if $b\prec_{R}\frac{1}{2}$ for some $b\in B$. Similarly define
other relations between blocks and $\pm\frac{1}{2}$.

This is a total order in all cases except when there is a unique $i\in[n]$
such that $-\frac{1}{2}\prec_{R}i\prec_{R}\frac{1}{2}$.

###### Proof.

The transitivity of this order is straightforward. If we show that $B<B$ is
not possible for any block, the order is well-defined. Suppose
$c_{0}\prec_{R}c_{1}\prec_{R}\cdots\prec_{R}c_{k}$ where $c_{0},c_{k}\in B$.
If $c_{1}\neq-c_{k}$ we get a contradiction to $c_{0}\sim c_{k}$ and similarly
if $c_{k-1}\neq-c_{1}$. So we must have $c_{1}=-c_{k}$ and $c_{k-1}=-c_{1}$.
But then we get $c_{1}\prec_{R}-c_{k}$ and $-c_{1}\prec_{R}c_{k}$, which is a
contradiction.

We will now show that this order is a total order in all cases except when
there is unique $i\in[n]$ such that
$-\frac{1}{2}\prec_{R}i\prec_{R}\frac{1}{2}$. The relationship of any block
with $\pm\frac{1}{2}$ is always defined. So we only have to check that any two
blocks are comparable. Let $B,B^{\prime}$ be boxed threshold blocks. If
$B,B^{\prime}$ are distinct blocks of the same sign, they are comparable by
definition of the equivalence relation. If they are of opposite sign and not
of the form $\\{i\\},\\{-i\\}$ for some $i\in[n]$, they are comparable since
there is some $a\in B$ and $b\in B^{\prime}$ of opposite signs such that
$|a|\neq|b|$. So we have to deal with the case $B=\\{i\\}$ and
$B^{\prime}=\\{-i\\}$ for some $i\in[n]$. If $i\succ_{R}\frac{1}{2}$ or
$i\prec_{R}-\frac{1}{2}$, $\\{-i\\}$ and $\\{i\\}$ are comparable.

Suppose $-\frac{1}{2}\prec_{R}i\prec_{R}\frac{1}{2}$ and $\\{i\\}$ and
$\\{-i\\}$ are blocks and they are not comparable. We have already seen that
all the positive blocks are comparable and so are all the negative blocks. Let
the order on the positive blocks be $B_{1}<B_{2}<\cdots<B_{k}$ and hence the
order on the negative blocks is $-B_{k}<\cdots<-B_{2}<-B_{1}$.

Suppose $B_{1}=\\{i\\}$. Since $i$ is not the only number satisfying
$-\frac{1}{2}\prec_{R}i\prec_{R}\frac{1}{2}$, we have $b\prec_{R}\frac{1}{2}$
for all $b\in B_{2}$. Taking some $b\in B_{2}$, since $i\notin B_{2}$, there
is some $k\in[n]$ such that $i\prec_{R}-k\prec_{R}b$. Since $\\{-i\\}$ is a
block, and $k\neq i$, there is some $l\in[n]$ such that
$-i\prec_{R}l\prec_{R}-k$ (we cannot have $-k\prec_{R}l\prec_{R}-i$ since this
would mean $\\{i\\}$ and $\\{-i\\}$ are comparable). But this gives
$k\prec_{R}-l\prec_{R}i$, which contradicts the fact that there is no positive
block less than $\\{i\\}$. We get a similar contradiction if $B_{k}=\\{i\\}$.

Suppose $B_{m}=\\{i\\}$ for some $m\in[2,k-1]$. Take some $b_{p}\in B_{m-1}$
and $b_{s}\in B_{m+1}$. There are three possibilities:

  1. (1)

There are $k_{p},k_{s}\in[n]$ such that

$b_{p}\prec_{R}-k_{p}\prec_{R}i\prec_{R}-k_{s}\prec_{R}b_{s}.$

Since $k_{p},k_{s}\neq i$, and $\\{-i\\}$ is a block, there are
$c_{p},c_{s}\in[n]$ such that

$\displaystyle-k_{p}$ $\displaystyle\prec_{R}c_{p}\prec_{R}-i$
$\displaystyle-i$ $\displaystyle\prec_{R}c_{s}\prec_{R}-k_{s}$

where the other possibilities are not possible since $\\{i\\}$ and $\\{-i\\}$
are not comparable. So we get

$b_{p}\prec_{R}-k_{p}\prec_{R}c_{p}\prec_{R}-i\prec_{R}c_{s}\prec_{R}-k_{s}\prec_{R}b_{s}.$

But this means that $c_{p}$ and $c_{s}$ are in blocks between $B_{m-1}$ and
$B_{m+1}$ and hence in $\\{i\\}$, which is a contradiction.

  2. (2)

There is no $k_{s}\in[n]$ such that $i\prec_{R}-k_{s}<b_{s}$. So we must have
$i\prec_{R}\frac{1}{2}\prec_{R}b_{s}$,
$-\frac{1}{2}\prec_{R}b_{p}\prec_{R}\frac{1}{2}$ and $k_{p}\in[n]$ such that

$b_{p}\prec_{R}-k_{p}\prec_{R}i\prec_{R}\frac{1}{2}\prec_{R}b_{s}.$

This is because $i$ satisfies $-\frac{1}{2}\prec_{R}i\prec_{R}\frac{1}{2}$ and
is not the only such number. Since $k_{p}\neq i$, and $\\{-i\\}$ is a block,
there is some $c_{p}\in[n]$ such that

$-k_{p}\prec_{R}c_{p}\prec_{R}-i$

where the other possibility is not possible since $\\{i\\}$ and $\\{-i\\}$ are
not comparable. So we get

$-\frac{1}{2}\prec_{R}b_{p}\prec_{R}-k_{p}\prec_{R}c_{p}\prec_{R}-i$

which, along with previous observations, gives

$i\prec_{R}-c_{p}\prec_{R}k_{p}\prec_{R}\frac{1}{2}\prec_{R}b_{s}.$

But this means $k_{p}$ is in a positive block after $\\{i\\}$ but before
$B_{m+1}$, which is a contradiction.

  3. (3)

The case when there is no $k_{p}\in[n]$ such that
$b_{p}\prec_{R}-k_{p}\prec_{R}i$ is handled similarly.

In the case when there is unique $i\in[n]$ such that
$-\frac{1}{2}\prec_{R}i\prec_{R}\frac{1}{2}$, the only order that is not
defined is between the blocks $\\{i\\}$ and $\\{-i\\}$ and the order is of the
form

(2)
$-B_{k}<\cdots<-B_{2}<-\frac{1}{2}<\\{i\\},\\{-i\\}<\frac{1}{2}<B_{2}<\cdots<B_{k}$

where $B_{2},\dots,B_{k}$ are blocks (not necessarily positive). This can be
proved using similar arguments as before. ∎

If there is a unique $i\in[n]$ such that
$-\frac{1}{2}\prec_{R}i\prec_{R}\frac{1}{2}$ and the blocks are ordered as in
(2), we make the convention that $\\{-i\\}<\\{i\\}$ if $B_{2}$ is a positive
block and $\\{i\\}<\\{-i\\}$ if $B_{2}$ is a negative block. Once this
convention is made, we get a total order on the boxed threshold blocks,
$\frac{1}{2}$ and $-\frac{1}{2}$ in all cases. The proof of the following
lemma follows from the definition of the order.

######  Lemma 3.7.

For any blocks $B,B^{\prime}$, $B<B^{\prime}$ implies that $-B^{\prime}<-B$.
Similarly, other usual properties of taking the negative hold; such as
$B<\frac{1}{2}$ implies that $-\frac{1}{2}<-B$.

######  Proposition 3.8.

The total order associated to a region (defined in Proposition 3.6) is always
of one of the following forms:

  1. (1)

$-B_{k}<\cdots<-B_{2}<-B_{1}<-\frac{1}{2}<\frac{1}{2}<B_{1}<B_{2}<\cdots<B_{k}$
where $B_{i}$ and $B_{i+1}$ are of opposite signs for all $i\in[k-1]$.

  2. (2)

$-B_{k}<\cdots<-B_{l+1}<-\frac{1}{2}<-B_{l}<\cdots<-B_{1}<B_{1}<\cdots<B_{l}<\frac{1}{2}<B_{l+1}<\cdots<B_{k}$
where the size of $B_{1}$ is greater than 1 and $B_{i}$ and $B_{i+1}$ are of
opposite signs for all $i\in[l-1]$ and $i\in[l+1,k-1]$.

  3. (3)

$-B_{k}<\cdots<-B_{2}<-\frac{1}{2}<-B_{1}<B_{1}<\frac{1}{2}<B_{2}<\cdots<B_{k}$
where the size of $B_{1}$ is 1, $B_{1}$ and $B_{2}$ are of the same sign, and
$B_{i}$ and $B_{i+1}$ are of opposite signs for all $i\in[2,k-1]$.

The association of such an order is a bijection between the regions of
$\mathcal{BT}_{n}$ and total orders of the types mentioned above.

###### Proof.

The fact that the first half of the order is the negative of the second
follows from Lemma 3.7. By the definition of blocks, for any two blocks of the
same sign which are in the same position with respect to $\pm\frac{1}{2}$,
there is some block of opposite sign between them.

The first form is when there is no $i\in[n]$ such that
$-\frac{1}{2}\prec_{R}i\prec_{R}\frac{1}{2}$. The third form is by the
convention we made at the end of the proof of Proposition 3.6, when there is a
unique $i\in[n]$ such that $-\frac{1}{2}\prec_{R}i\prec_{R}\frac{1}{2}$.

So we have to show that when there is more than one $i\in[n]$ such that
$-\frac{1}{2}\prec_{R}i\prec_{R}\frac{1}{2}$, the block $B_{1}$ has size
greater than 1. Suppose $B_{1}=\\{a\\}$ for some $a\in N$. Let $b\in B_{2}$.
$b$ has the same sign as $-a$ and is in a different block. Hence there is some
$k\in N$ of opposite sign to $-a$ and $b$ such that $-a\prec_{R}k\prec_{R}b$.
But this means $k$ is in a block between $\\{-a\\}$ and $B_{2}$ and hence in
$\\{a\\}$, which is a contradiction since $k\neq a$.

We will now show that the association of such an order is a bijection between
the regions of $\mathcal{BT}_{n}$ and total orders of the types mentioned
above. First note that, from the total order associated to a region, we can
get back the inequalities describing the region as follows: For any $i\neq j$
in $[n]$, $x_{i}+x_{j}>0$ if and only if the block containing $-j$ is before
that containing $i$ and the relationship between $x_{i}$ and $\pm\frac{1}{2}$
is obtained in the same way. On the other hand, given an order of one of the
forms given above, choosing some real numbers $0<c_{1}<c_{2}<\cdots<c_{k}$
such that $c_{i}<\frac{1}{2}$ if $B_{i}<\frac{1}{2}$ and $c_{i}>\frac{1}{2}$
if $B_{i}>\frac{1}{2}$ and putting $x_{a}=c_{i}$ for all $a\in B_{i}$ for all
$i\in[k]$, gives a point satisfying the required inequalities. It can also be
shown that different such orders correspond to different regions. ∎

Hence, to count the number of regions of $\mathcal{BT}_{n}$, we just have to
count the number of orders of the forms mentioned in Proposition 3.6. Note
that the first half of the order is the negative of the second, so we just
consider the second half. That is, we count orders of the following types,
where in all cases, $B_{1},\dots,B_{k}$ is a partition of $[n]$ with a sign
assigned to each block:

  1. (1)

$\frac{1}{2}<B_{1}<B_{2}<\cdots<B_{k}$ where $B_{i}$ and $B_{i+1}$ are of
opposite signs for all $i\in[k-1]$.

  2. (2)

$B_{1}<\cdots<B_{l}<\frac{1}{2}<B_{l+1}<\cdots<B_{k}$ where the size of
$B_{1}$ is greater than 1 and $B_{i}$ and $B_{i+1}$ are of opposite signs for
all $i\in[l-1]$ and $i\in[l+1,k-1]$.

  3. (3)

$B_{1}<\frac{1}{2}<B_{2}<\cdots<B_{k}$ where the size of $B_{1}$ is 1, $B_{1}$
and $B_{2}$ are of the same sign, and $B_{i}$ and $B_{i+1}$ are of opposite
signs for all $i\in[2,k-1]$.

######  Proposition 3.9.

The total number of orders of the forms mentioned above is

$4\cdot a(n)+\sum_{k=1}^{n}4(k!-(k-1)!)(k\cdot S(n,k)-n\cdot S(n-1,k-1)).$

Here $a(n)$ is the $n^{th}$ ordered Bell number.

###### Proof.

We will count the number of orders of each of the above forms.

  1. (1)

In the first case, we just have to define an ordered partition of $[n]$ and
assign alternating signs to them. The number of ways this can be done is

$\sum_{k=1}^{n}\sum\limits_{\begin{subarray}{c}(a_{1},\dots,a_{k})\\\
a_{1}+\cdots+a_{k}=n\end{subarray}}2\cdot\dfrac{n!}{a_{1}!\cdots
a_{k}!}=2\cdot a(n).$

  2. (2)

In the second case, we consider two sub-cases:

    1. (a)

There is no block after $\frac{1}{2}$. In this case, we have to define an
ordered partition of $[n]$ whose first part has size greater that $1$ and
assign alternating signs to them. The number of ways this can be done is

$\sum_{k=1}^{n-1}\sum\limits_{\begin{subarray}{c}(a_{1},\dots,a_{k})\\\
a_{1}+\cdots+a_{k}=n,\ a_{1}\neq 1\end{subarray}}2\cdot\dfrac{n!}{a_{1}!\cdots
a_{k}!}=2(a(n)-n\cdot a(n-1))$

where the equality is because the number of ordered partitions of $[n]$ with
first block having size $1$ is $n\cdot a(n-1)$.

    2. (b)

There is some block after $\frac{1}{2}$. In this case, we have to again define
an ordered partition of $[n]$ whose first part has size greater that $1$. But
we then have to choose a spot between two blocks to place $\frac{1}{2}$ and
then choose a sign for the first block and the block after $\frac{1}{2}$. The
number of ways this can be done is

$\sum_{k=1}^{n-1}\sum\limits_{\begin{subarray}{c}(a_{1},\dots,a_{k})\\\
a_{1}+\cdots+a_{k}=n,\ a_{1}\neq 1\end{subarray}}4(k-1)\dfrac{n!}{a_{1}!\cdots
a_{k}!}.$

Making the following substitution for all $k\in[n-1]$

$\displaystyle\sum\limits_{\begin{subarray}{c}(a_{1},\dots,a_{k})\\\
a_{1}+\cdots+a_{k}=n,\ a_{1}\neq 1\end{subarray}}\dfrac{n!}{a_{1}!\cdots
a_{k}!}$ $\displaystyle=\sum\limits_{\begin{subarray}{c}(a_{1},\dots,a_{k})\\\
a_{1}+\cdots+a_{k}=n\end{subarray}}\dfrac{n!}{a_{1}!\cdots a_{k}!}\hskip
5.69046pt-\sum\limits_{\begin{subarray}{c}(1,a_{2},\dots,a_{k})\\\
1+a_{2}+\cdots+a_{k}=n\end{subarray}}\dfrac{n!}{1!a_{2}!\cdots a_{k}!}$
$\displaystyle=k!\cdot S(n,k)-n\cdot(k-1)!\cdot S(n-1,k-1)$

we get that the initial expression is the same as

$\sum_{k=1}^{n}4(k!-(k-1)!)(k\cdot S(n,k)-n\cdot S(n-1,k-1)).$

  3. (3)

In the third case, we have to choose the element of $[n]$ in $B_{1}$ and then
define an ordered partition of the remaining $(n-1)$ elements and assign
alternating signs to them. Since we want $B_{1}$ and $B_{2}$ to have the same
sign, we just need to assign a sign to $B_{2}$. So, the number of orders of
the third type is

$n\cdot\sum_{k=1}^{n-1}\sum\limits_{\begin{subarray}{c}(a_{1},\dots,a_{k})\\\
a_{1}+\cdots+a_{k}=n-1\end{subarray}}2\cdot\dfrac{(n-1)!}{a_{1}!\cdots
a_{k}!}=n\cdot 2\cdot a(n-1).$

Adding up the counts made for each form gives us the required result. ∎

So, from the observations made above, we have proved the following theorem:

######  Theorem 3.10.

The number of regions of $\mathcal{BT}_{n}$ is

$\displaystyle 4\cdot a(n)+\sum_{k=1}^{n}4(k!-(k-1)!)(k\cdot S(n,k)-n\cdot
S(n-1,k-1)).$

Similar arguments can be applied to the threshold arrangement
$\mathcal{T}_{n}$.

######  Proposition 3.11.

The regions of $\mathcal{T}_{n}$ are in bijection with ordered partitions of
$[-n,n]\setminus\\{0\\}$ of the form

$-B_{k}<\cdots<-B_{2}<-B_{1}<B_{1}<B_{2}<\cdots<B_{k}$

where all elements of each block have the same sign, the size of $B_{1}$ is
greater than $1$ and $B_{i}$ and $B_{i+1}$ are of opposite signs for all
$i\in[k-1]$.

The bijection in Proposition 3.11 is defined just as was done for
$\mathcal{BT}_{n}$. That is, the region associated to such an order satisfies
$x_{i}+x_{j}>0$ if and only if $-j$ appears before $i$ in the order. Again,
such orders are completely specified by their second half, which are ordered
partition of $[n]$ with first block size greater than $1$ and a sign assigned
to the first block (and alternate signs for consecutive blocks). So, we get
the following theorem:

######  Theorem 3.12.

The number of regions of $\mathcal{T}_{n}$ is

$2(a(n)-n\cdot a(n-1)).$

######  Remark 3.13.

It can be shown that the orders considered in Proposition 3.11 are in
bijection with the set of threshold pairs (of size $n$) in standard form111A
pair $(\pi,\omega)$ is a _threshold pair_ (of size $n$) if $\pi$ is a
permutation of size $n$ and $\omega$ is a word in $\\{+1,-1\\}^{n}$. A
threshold pair $(\pi,\omega)$ of size $n\geq 2$ is in _standard form_ if
$\omega_{1}=\omega_{2}$ and if $\omega_{i}=\omega_{i+1}$ implies
$\pi_{i}<\pi_{i+1}$ for all $1\leq i<n$. considered by Spiro [8]. In fact, he
showed that the threshold pairs are in bijection with the labeled threshold
graphs.

The known formula for the number of labeled threshold graphs is in terms of
Eulerian numbers. Hence for the sake of completeness, we now show that the
formula in Theorem 3.12 is the same as the one containing Eulerian numbers. We
first recall a few identities related to Eulerian numbers $A(n,k)$ and the
ordered Bell number $a(n)$. A quick reference for these identities is Bóna’s
book [3, Section 1.1].

  * •

$a(n)=\sum\limits_{k=0}^{n-1}2^{k}\cdot A(n,k)$

  * •

$A(n,0)=1$, $A(n,n-1)=1$ and $A(n,k)=0$ for all $k\geq n$.

  * •

$A(n,k)=(n-k)A(n-1,k-1)+(k+1)A(n-1,k)$ for $k\geq 1$.

######  Proposition 3.14.

We have the following equality

$r(\mathcal{T}_{n})=\sum\limits_{k=1}^{n-1}2^{k}(n-k)A(n-1,k-1).$

###### Proof.

From Theorem 3.12, we have

$\begin{split}r(\mathcal{T}_{n})&=2(a(n)-n\cdot a(n-1))\\\
&=2\cdot\Bigg{(}\sum\limits_{k=0}^{n-1}2^{k}\cdot
A(n,k)-n\cdot\Big{(}\sum\limits_{k=0}^{n-2}2^{k}\cdot
A(n-1,k)\Big{)}\Bigg{)}\\\
\frac{r(\mathcal{T}_{n})}{2}&=A(n,0)+\sum\limits_{k=1}^{n-1}2^{k}\big{(}(n-k)A(n-1,k-1)+(k+1)A(n-1,k)\big{)}\\\
&\hskip 14.22636pt-n\Big{(}\sum\limits_{k=0}^{n-2}2^{k}\cdot
A(n-1,k)\Big{)}\\\
&=\sum\limits_{k=1}^{n-1}2^{k}(n-k)A(n-1,k-1)+\Big{(}A(n,0)+\sum\limits_{k=1}^{n-1}2^{k}(k+1)A(n-1,k)\Big{)}\\\
&\hskip 14.22636pt-\sum\limits_{k=0}^{n-2}n\cdot 2^{k}\cdot
A(n-1,k)\end{split}$

Since $A(n-1,n-1)=0$,

$\begin{split}\frac{r(\mathcal{T}_{n})}{2}&=\sum\limits_{k=1}^{n-1}2^{k}(n-k)A(n-1,k-1)+\sum\limits_{k=0}^{n-2}2^{k}(k+1)A(n-1,k)-\sum\limits_{k=0}^{n-2}n\cdot
2^{k}\cdot A(n-1,k)\\\
&=\sum\limits_{k=1}^{n-1}2^{k}(n-k)A(n-1,k-1)-\sum\limits_{k=0}^{n-2}(n-k-1)2^{k}\cdot
A(n-1,k)\end{split}$

Replacing $k+1$ by $t$ in the second sum, we get

$\begin{split}r(\mathcal{T}_{n})&=2\cdot\Bigg{(}\sum\limits_{k=1}^{n-1}2^{k}(n-k)A(n-1,k-1)-\sum\limits_{t=1}^{n-1}2^{t-1}(n-t)A(n-1,t-1)\Bigg{)}\\\
&=\sum\limits_{k=1}^{n-1}2^{k}(n-k)A(n-1,k-1).\end{split}$

∎

##  4\. The colored threshold graphs

Before defining the colored threshold graphs that are in bijection with the
regions of the boxed threshold arrangement, we recall the bijection between
regions of the threshold arrangement and labeled threshold graphs.

######  Definition 4.1.

A threshold graph is defined recursively as follows:

  1. (1)

The empty graph is a threshold graph.

  2. (2)

A graph obtained by adding an isolated vertex to a threshold graph is a
threshold graph.

  3. (3)

A graph obtained by adding a vertex adjacent to all vertices of a threshold
graph is a threshold graph.

######  Definition 4.2.

A labeled threshold graph is a threshold graph having $n$ vertices with
vertices labeled distinctly using $[n]$.

Such labeled threshold graphs can be specified by a signed permutation of
$[n]$, that is, a permutation of $[n]$ with a sign associated to each number.
The signed permutation $i_{1}\ i_{2}\ \cdots\ i_{n}$ would correspond to the
labeled threshold graph obtained by adding vertices labeled
$|i_{1}|,|i_{2}|,\dots,|i_{n}|$ in order where a positive $i_{k}$ means that
$|i_{k}|$ is added adjacent to all previous vertices and a negative $i_{k}$
means that it is added isolated to the previous vertices. A maximal string of
positive numbers or negative numbers in a signed permutation will be called a
block.

######  Example 4.3.

The labeled threshold graph associated to the signed permutation on $[5]$
given by
$\overset{-}{2}\overset{-}{3}\overset{+}{1}\overset{+}{4}\overset{-}{5}$ is
shown in Figure 1.

$2$4$\xrightarrow{\overset{-}{3}}$

$2$$3$4$\xrightarrow{\overset{+}{1}}$

$2$$3$$1$4$\xrightarrow{\overset{+}{4}}$

$2$$3$$1$$\xrightarrow{\overset{-}{5}}$$4$

$2$$3$$1$$4$$5$

Figure 1. Construction of threshold graph corresponding to
$\overset{-}{2}\overset{-}{3}\overset{+}{1}\overset{+}{4}\overset{-}{5}$.

The following facts can be verified:

  1. (1)

The sign of the first number in the permutation does not matter and hence we
can make the first block have size greater than 1.

  2. (2)

Elements in the same block can be reordered.

Hence, labeled threshold graphs can be specified by an ordered partition of
$[n]$ with first block size greater than 1 and alternating signs assigned to
the blocks. In fact, this association is a bijection.

Given a threshold graph $G_{1}$, we can obtain this alternating signed ordered
partition of $[n]$ as follows: Since $G_{1}$ is a threshold graph, it has at
least one isolated vertex or at least one vertex that is adjacent to all other
vertices. If it has an isolated vertex, set $D_{1}$ to be the set of all
isolated vertices, assign it a negative sign and set $G_{2}$ to be the graph
obtained by deleting all the vertices of $D_{1}$ from $G_{1}$. If $G_{1}$ has
at least one vertex adjacent to all other vertices, set $D_{1}$ to be the set
of all such vertices, assign it a positive sign and set $G_{2}$ to be the
graph obtained by deleting all the vertices of $D_{1}$ from $G_{1}$. We repeat
this process until we obtain a graph $G_{k}$ which is complete, in which case
we set $D_{k}$ to be all vertices of $G_{k}$ and assign it a positive sign, or
$G_{k}$ has no edges, in which case we set $D_{k}$ to be all vertices of
$G_{k}$ and assign it a negative sign. Then set $B_{i}=D_{k-i+1}$ and assign
it the same sign as $D_{k-i+1}$ for all $i\in[k]$. The signed ordered
partition $B_{1},\dots,B_{k}$ is the one associated to $G_{1}$.

######  Example 4.4.

Figure 2 shows an example of obtaining the signed blocks from a threshold
graph. The corresponding signed ordered partition for this example is
$\overset{-}{\\{2,3\\}}\overset{+}{\\{1,4\\}}\overset{-}{\\{5\\}}$.

$2$$3$$1$$4$$5$$\longrightarrow$$D_{2}=\overset{+}{\\{1,4\\}}$

$2$$3$$1$$\longrightarrow$$4$$D_{1}=\overset{-}{\\{5\\}}$

$2$$3$4$\longrightarrow$$D_{2}=\overset{+}{\\{1,4\\}}$

$D_{3}=\overset{-}{\\{2,3\\}}$

Figure 2. Obtaining blocks from a threshold graph.

Hence, regions of $\mathcal{T}_{n}$ and labeled threshold graphs on $n$
vertices are both in bijection with ordered partitions of $[n]$ with first
block size greater than 1 and alternating signs assigned to the blocks. Hence
we obtain a bijection between regions of $\mathcal{T}_{n}$ and labeled
threshold graphs on $n$ vertices. By combining the definitions of the two
bijections we see that to a labeled threshold graph on $n$ vertices we assign
the region where $x_{i}+x_{j}>0$ if and only if there is an edge between $i$
and $j$.

This can be proved as follows: If $-B_{k}<\cdots<-B_{1}<B_{1}<\cdots<B_{k}$ is
the threshold block order corresponding to some region $R$ of
$\mathcal{T}_{n}$, $-j\prec_{R}i$ for some $i\neq j$ in $[n]$ if and only if
one of the following holds:

  1. (1)

$-j$ and $i$ both appear in $B_{1},\dots,B_{k}$ with $-j$ appearing first.

  2. (2)

$-j$ appears in $-B_{k},\dots,-B_{1}$ and $i$ appears in $B_{1},\dots,B_{k}$.

  3. (3)

$-j$ and $i$ both appear in $-B_{k},\dots,-B_{1}$ with $-j$ appearing first.

Equivalently, one of the following holds:

  1. (1)

$-j$ and $i$ both appear in $B_{1},\dots,B_{k}$ with $-j$ appearing first.

  2. (2)

$i$ and $j$ both appear in $B_{1},\dots,B_{k}$.

  3. (3)

$-i$ and $j$ both appear in $B_{1},\dots,B_{k}$ with $-i$ appearing first.

But this is precisely the condition for there to be an edge between $i$ and
$j$ in the threshold graph corresponding to $B_{1}<\cdots<B_{k}$.

We now move on to the boxed threshold arrangement.

######  Definition 4.5.

A colored threshold graph is defined recursively as follows:

  1. (1)

The empty graph is a colored threshold graph.

  2. (2)

A graph obtained by adding an isolated vertex to a colored threshold graph is
a colored threshold graph. If there are colored vertices in the initial
colored threshold graph, the new vertex should be colored red. If not, the new
vertex can be left uncolored or colored red.

  3. (3)

A graph obtained by adding a vertex adjacent to all vertices of a colored
threshold graph is a colored threshold graph. If there are colored vertices in
the initial colored threshold graph, the new vertex should be colored blue. If
not, the new vertex can be left uncolored or colored blue.

######  Definition 4.6.

A labeled colored threshold graph is a colored threshold graph with $n$
vertices with the vertices labeled distinctly with elements of $[n]$.

Just as for threshold graphs, labeled colored threshold graphs can be
represented as a signed permutation. However, we also have to specify if and
when the coloring of the vertices starts. This is done by using the symbol
$\frac{1}{2}$. Having $\frac{1}{2}$ at the end of the signed permutation means
that none of the vertices should be colored.

######  Example 4.7.

The sequence
$\overset{+}{2}\frac{1}{2}\overset{+}{1}\overset{+}{3}\overset{-}{4}\overset{-}{5}$
corresponds to the graph shown in Figure 3.

$2$$1$$3$$4$$5$ Figure 3. Labeled colored threshold graph corresponding to
$\overset{+}{2}\frac{1}{2}\overset{+}{1}\overset{+}{3}\overset{-}{4}\overset{-}{5}$.

Using similar observations about these sequences associated to labeled colored
threshold graphs as done for labeled threshold graphs, we get that labeled
colored threshold graphs are in bijection with orders of the forms counted in
Proposition 3.9. Since these orders also correspond to region of
$\mathcal{BT}_{n}$, we get a bijection between labeled colored threshold
graphs with $n$ vertices and regions of $\mathcal{BT}_{n}$. Just as before,
the inequalities describing the region associated to a colored threshold graph
are as follows: $x_{i}+x_{j}>0$ if and only if there is an edge between $i$
and $j$, $-\frac{1}{2}<x_{i}<\frac{1}{2}$ if $i$ is not colored,
$x_{i}>\frac{1}{2}$ if $i$ is colored blue and $x_{i}<-\frac{1}{2}$ if $i$ is
colored red. Notice that the underlying labeled threshold graph corresponds to
the $\mathcal{T}_{n}$ region that the $\mathcal{BT}_{n}$ region lies in.

Also, we can see that the bounded regions of $\mathcal{BT}_{n}$ are in
bijection with the regions of $\mathcal{T}_{n}$. Both are represented by
labeled threshold graphs with $n$ vertices. The bounded region of
$\mathcal{BT}_{n}$ corresponding to a region of $\mathcal{T}_{n}$ is the one
satisfying the same inequalities between $x_{i}+x_{j}$ and $0$ for all $i\neq
j$ in $[n]$ and having $-\frac{1}{2}<x_{i}<\frac{1}{2}$ for all $i\in[n]$.

$x_{1}=-\frac{1}{2}$$x_{1}=\frac{1}{2}$$x_{2}=-\frac{1}{2}$$x_{2}=\frac{1}{2}$$x_{1}+x_{2}=0$122121212121212121212121
Figure 4. Regions of $\mathcal{BT}_{2}$ represented by labeled colored
threshold graphs

## References

  * [1] C. A. Athanasiadis, Characteristic polynomials of subspace arrangements and finite fields, Adv. Math. 122 (1996), 193–233.
  * [2] C. A. Athanasiadis, Extended Linial hyperplane arrangements for root systems and a conjecture of Postnikov and Stanley, J. Algebraic Combin. 10 (1999), 207–225.
  * [3] M. Bóna, Combinatorics of Permutations, Second Edition, CRC Press, 2012.
  * [4] N. J. A. Sloane et al., The On-Line Encyclopedia of Integer Sequences, 2021\. Available at https://oeis.org.
  * [5] J. Song, Enumeration of graphs and the characteristic polynomial of the hyperplane arrangements $\mathcal{J}_{n}$, J. Korean Math. Soc. 54 (2017), 1595–1604.
  * [6] J. Song, On certain hyperplane arrangements and colored graphs, Bull. Korean Math. Soc. 54 (2017), 375–382.
  * [7] J. Song, Characteristic polynomial of the hyperplane arrangements $\mathcal{J}_{n}$ via finite field method, Commun. Korean Math. Soc. 33 (2018), 759–765.
  * [8] S. Spiro, Counting labeled threshold graphs with Eulerian numbers, Australas. J. Combin. 77 (2020), 249–255.
  * [9] R. P. Stanley, An introduction to hyperplane arrangements, in E. Miller, V. Reiner, and B. Sturmfels, eds., Geometric Combinatorics, Amer. Math. Soc., 2007, pp. 389–496.
  * [10] T. Zaslavsky, Facing up to arrangements: face-count formulas for partitions of space by hyperplanes, Mem. Amer. Math. Soc. 1 (1975), No. 154.


11institutetext: Department of Computer Science, University of North Carolina
at Chapel Hill

# NeuralExplorer: State Space Exploration of Closed Loop Control Systems Using
Neural Networks

Manish Goyal  Parasara Sridhar Duggirala
<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a framework for performing state space exploration
of closed loop control systems. Our approach involves approximating
sensitivity and a newly introduced notion of inverse sensitivity by a neural
network. We show how the approximation of sensitivity and inverse sensitivity
can be used for computing estimates of the reachable set. We then outline
algorithms for performing state space exploration by generating trajectories
that reach a neighborhood. We demonstrate the effectiveness of our approach by
applying it not only to standard linear and nonlinear dynamical systems, but
also to nonlinear hybrid systems and also neural network based feedback
control systems.

###### Keywords:

S

tate space exploration, Sensitivity, Inverse Sensitivity, Neural Networks,
Testing, Approximation, Falsification

##  1 Introduction

The decrease in the cost of embedded microcontrollers and integration of
software into many control processes has lead to the deployment of
increasingly sophisticated control algorithms in many safety critical
scenarios. Additionally, the recent improvements in software and hardware
platforms for training and evaluation of neural networks have made it easier
to integrate them in embedded devices. This increased complexity of software
in control tasks makes testing and validation of the closed loop systems very
challenging.

Testing is the most commonly used technique for checking whether a system
satisfies its specification. In a typical work flow, after designing the
feedback function, the control designer generates a few test cases for the
closed loop system and checks if they satisfy the specification. Given that
the state space is continuous, finding the trajectory that violates the
specification is similar to searching for a needle in a haystack. For example,
consider a regulation application where the output of the control system is
stabilized to a set point $s$ with the error threshold of $\delta$. Therefore,
the output should remain in the interval $[s-\delta,s+\delta]$. If all the
test cases satisfy the specification, the control designer does not know the
next test case input that results in a higher value of error than observed in
the test cases.

While the designer can encode the property as a temporal logic formula and use
off-the-shelf falsification tools, such an approach has some drawbacks. First,
falsification tools are geared towards finding a trajectory that violates the
given specification, not necessarily to help the control designer in state
space exploration. Second, if the threshold of $\delta$ is increased, the
results from the falsification analysis are no longer useful. Finally, using
falsification tools would require the specification to be given in a temporal
logic such as signal temporal logic or metric temporal logic. While such
specification might be useful in the verification and deployment phase, they
are sometimes a hinderance during the design phase. What the control designer
needs, is a tool that helps her to systematically generate trajectories that
result in various values of error thresholds. Currently, there are no tools
for performing state space exploration based on the previous tests generated
or based on the input given by the control designer.

In this paper, we present _NeuralExplorer_ , a technique for performing state
space exploration of closed loop control systems using neural networks. After
the control designer generates a few test scenarios, she can use
NeuralExplorer to generate test case scenarios that reach a neighborhood. The
artifact that helps us in this endeavor is _sensitivity_. Informally,
sensitivity of a closed loop system is the change in the trajectory of the
system as a result of perturbing the initial condition. In this paper, we also
present a backward time notion of sensitivity called _inverse sensitivity_.
Given a sample set of trajectories in the domain of interest, we train a
neural network to learn the sensitivity and the inverse sensitivity functions.
This neural network is then used to generate a trajectory (or trajectories)
that reaches a destination (or a neighborhood around it).

Our framework has two primary advantages. First, since NeuralExplorer relies
only on the sample test cases, it does not require a model of the system and
can be applied to a black-box systems. Second, since sensitivity is a
fundamental property of the closed loop system, approximating it using a
neural network is generalizable to trajectories that are beyond the test cases
generated by the control designer. In evaluating our framework, we were not
only able to perform state space exploration for standard linear and nonlinear
dynamical systems, but also for nonlinear hybrid systems and neural network
based feedback control systems. We believe that NeuralExplorer is useful for
generating corner cases and supplements some of the existing testing and
reachable set computation procedures.

##  2 Related Work

Reachability analysis is often employed for proving the safety specification
of safety critical control system [11, 4]. Some of the recent works in this
domain are SpaceEx [22], Flow* [9], CORA [3] and HyLAA [7]. These techniques
use a symbolic representation for the reachable set of states. While these are
useful for proving that the safety specification is satisfied, generating
counterexamples using reachability analysis is still an area of research [24,
26].

For generating trajectories that violate a given safety specification,
falsification techniques are widely applied [19, 15]. In these techniques, the
required specification is specified in a temporal logic such as Metric
Temporal Logic (MTL) [31] or Signal Temporal Logic (STL) [35, 32]. Given a
specification, falsification techniques generate several sample trajectories
and use various heuristics [38, 2, 43, 48, 12, 23] for generating trajectories
that violate the specification. Prominent tools in this domain include
S-Taliro [5] and Breach [13].

Bridging falsification and reachability are simulation driven verification
methods [14, 16, 28, 20]. These methods compute an upper bound on the
sensitivity of the trajectories and compute an over-approximation of the
reachable set using sample trajectories. While these techniques bridge the gap
between falsification and verification, they suffer from curse of
dimensionality. That is, the number of trajectories generated might increase
exponentially with the dimensions. C2E2 [17], and DryVR [21] are some of the
well known tools in this domain.

Given the rich history of application of neural networks in control [36, 33,
37] and the recent advances in software and hardware platforms, neural
networks are now being deployed in various control tasks. As a result, many
verification techniques are now being developed for neural network based
control systems [30, 47, 44, 18]. Additionally, techniques for verification of
neural networks deployed in other domains have been proposed [46, 45, 27].

In this paper, we use neural networks to approximate an underlying property of
sensitivity. We believe that this is a valid approach because recently, many
neural network based frameworks for learning the dynamics or their properties
have been proposed [39, 40, 34, 42, 8, 41].

##  3 Preliminaries

We denote the elements of the state space as $x$ to be elements in
${\mathbb{R}}^{n}$. Vectors are denoted as $v$. We denote the dynamics of the
plant as

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

Where $x$ is the state space of the system that evolves in ${\mathbb{R}}^{n}$
and $u$ is the input space in ${\mathbb{R}}^{m}$.

######  Definition 1 (Unique Trajectory Feedback Functions)

A feedback function $u=g(x)$ is said to be unique trajectory feedback function
if the closed loop system $\dot{x}=f(x,g(x))$ is guaranteed existence and
uniqueness of the solution for the initial value problem for all initial
points $x_{0}\in{\mathbb{R}}^{n}$.

Notice that for a feedback function to give a unique trajectory feedback, it
need not be differentiable. From the sufficuent conditions of ODE solutions,
it is sufficient if $g(x)$ is continuous and is lipschitz.

######  Definition 2 (Trajectories of Closed Loop System)

Given a unique trajectory feedback function $u=g(x)$, a trajectory of closed
loop system $\dot{x}=f(x,g(x))$, denoted as $\xi_{g}(x_{0},t)$ ($t\geq 0$), is
the solution of the initial value problem of the differential equation
$\dot{x}=f(x,g(x))$ with initial condition $x_{0}$. We often drop the feedback
function $g$ when it is clear from the context.

We extend the notion of trajectory to include backward time trajectories as
well. Given $t>0$, the backward time trajectory $\xi_{g}(x_{0},-t)=x$ such
that $\xi_{g}(x,t)=x_{0}$. We denote backward time trajectory as
$\xi^{-1}(x_{0},t)$.

Given $x_{0},x_{1}\in{\mathbb{R}}^{n}$ and $t>0$ such that
$\xi(x_{0},t)=x_{1}$, then $\xi^{-1}(x_{1},t)=x_{0}$. It is trivial to observe
that $\xi(\xi^{-1}(x_{0},t),t)=x_{0}$.

######  Definition 3 (Sensitivity of Trajectories)

Given an initial state $x_{0}$, vector $v$, and time $t$, the sensitivity of
the trajectories, denoted as $\Phi(x_{0},v,t)$ is defined as.

$\Phi(x_{0},v,t)=\xi(x_{0}+v,t)-\xi(x_{0},t).$

Informally, sensitivity is the vector difference between the trajectories
starting from $x_{0}$ and $x_{0}+v$ after time $t$. We extend the definition
of sensitivity to backward time trajectories as

$\Phi^{-1}(x_{0},v,t)=\xi^{-1}(x_{0}+v,t)-\xi^{-1}(x_{0},t).$

We call $\Phi^{-1}(x_{0},v,t)$ as inverse sensitivity function. Informally,
inverse sensitivity function gives us the perturbation of the initial
condition that is required to displace the trajectory passing through $x_{0}$
by $v$. Observe that $\xi(\xi^{-1}(x_{0},t)+\Phi^{-1}(x_{0},v,t),t)=x_{0}+v$.

For general nonlinear differential equations, analytic representation of the
trajectories of the ODEs need not exist. If the closed loop system is a smooth
function, then the infinite series for the trajectories is given as

$\xi(x_{0},t)=x_{0}+{\cal L}_{f}(x_{0})t+{\cal
L}_{f}^{2}(x_{0})\frac{t^{2}}{2!}+{\cal
L}_{f}^{3}(x_{0})\frac{t^{3}}{3!}+\ldots.$ (2)

Where ${\cal L}_{f}^{i}$ is the $i^{th}$ order Lie-derivative over the field
$f(x,g(x))$ at the state $x_{0}$. Hence, one can write the sensitivity
function as

$\Phi(x_{0},v,t)=v+({\cal L}_{f}(x_{0}+v)-{\cal L}_{f}(x_{0}))t+({\cal
L}_{f}^{2}(x_{0}+v)-{\cal L}_{f}^{2}(x_{0}))\frac{t^{2}}{2!}+\ldots.$ (3)

$\Phi^{-1}(x_{0},v,t)$ is obtained by substituting $-f$ for $f$ in Equation 3.
When the closed loop dynamics is linear, i.e., $\dot{x}=Ax$, it is easy to
observe that $\Phi(x_{0},v,t)=e^{At}v$, $\Phi^{-1}(x_{0},v,t)=e^{-At}v$ where
$e^{At}$ ($e^{-At}$) is the matrix exponential of the matrix $At$ ($-At$).
Obseve that for linear systems, the inverse sensitivity function is
independent of the state $x_{0}$.

For nonlinear dynamical systems, one can truncate the infinite series upto a
specific order and obtain an approximation. However, for hybrid systems that
have state based mode switches, or for feedback functions where the closed
loop dynamics is not smooth or is discontinuous, such an infinite series
expansion is hard to compute. The central idea in this paper is to approximate
$\Phi$ and $\Phi^{-1}$ using a neural network and perform state space
exploration using such neural networks.

##  4 Neural Network Approximations of Sensitivity and Inverse Sensitivity

Given a domain of operation $D\subseteq{\mathbb{R}}^{n}$, one can generate a
finite set of trajectories for testing the system operation in $D$. Often,
these trajectories are generated using numerical ODE solvers which return
trajectories sampled at a regular time step. For approximating sensitivity and
inverse sensitivity, we generate a finite number of time bounded trajectories
where the step size, time bound, and the number of trajectories are specified
by the user. The trajectories can be generated either according to a
probability distribution specified by the user or from specific initial
configurations provided by her.

Given a sampling of a trajectory at regular time interval with step size $h$,
i.e., $\xi(x_{0},0)$, $\xi(x_{0},h)$, $\xi(x_{0},2h)$, $\ldots$,
$\xi(x_{0},kh)$, we make two observations. First, any prefix or suffix of this
sequence is also a trajectory, albeit, of a shorter duration. Hence, from a
given set of trajectories, one can generate more _virtual trajectories_ that
have shorter duration. Second, given two trajectories (real or virtual)
starting from initial states $x_{1}$ and $x_{2}$, ($x_{1}\neq x_{2}$), We have
the two following observations.

$\displaystyle\Phi(x_{1},x_{2}-x_{1},t)$ $\displaystyle=$
$\displaystyle\xi(x_{2},t)-\xi(x_{1},t)$ (4)
$\displaystyle\Phi^{-1}(\xi(x_{1},t),\xi(x_{2},t)-\xi(x_{1},t),t)$
$\displaystyle=$ $\displaystyle x_{2}-x_{1}.$ (5)

Given an initial set of trajectories, we generate virtual trajectories and use
Equations 4 and 5 for generating all tuples $\langle x_{0},v,t,v_{sen}\rangle$
and $\langle x_{0},v,t,v_{isen}\rangle$ such that $v_{sen}=\Phi(x_{0},v,t)$
and $v_{isen}=\Phi^{-1}(x_{0},v,t)$. This data is then used for training and
evaluation of the neural network to approximate the inverse sensitivity
function. We denote these networks as $\mathit{NN}_{\Phi}$ and
$\mathit{NN}_{\Phi^{-1}}$ respectively.

###  4.1 Evaluation on Standard Benchmarks

For approximating the sensitivity and inverse sensitivity functions, we pick a
standard set of benchmakrs consisting of nonlinear dynamical systems, hybrid
systems, and control systems with neural network feedback functions. Most of
the benchmarks are taken from standard hybrid systems benchmark suite [1, 29,
6]. The benchmarks Brussellator, Lotka, Jetengine, Buckling, Vanderpol,
Lacoperon, Roesseler, Steam, Lorentz, and Coupled vanderpol are continuous
nonlinear systems, where Lorentz and Roesseler are chaotic as well.
SmoothHybrid- -Oscillator and HybridOscillator are nonlinear hybrid systems.
The remaining benchmarks Mountain Car and Quadrotor are selected from [30],
where the state feedback controller is given in the from of neural network.

For each benchmark, we generated a given number ($N$ is typically 30 or 50) of
trajectories, where the step size for ODE solvers ($s$) and the time bound are
provided by the user. The data used for training the neural network is
collected as described in 4. We use 90% of the data for training and 10% for
testing. We used the Python Multilayer Perceptron implemented as Sequential
model in Keras [10] library with Tensorflow as the backend. The network has 8
layers with each layer having 512 neurons. The optimizer used is stochastic
gradient descent. The network is trained using Levenberg-Marquardt
backpropagation algorithm optimizing the mean absolute error loss function,
and the Nguyen-Widrow initialization.

The activation function used to train the network is relu for all benchmarks
except Mountain car for which sigmoid performs better because the NN
controller is sigmoid-based. Note that the choice of hyper-parameters such as
number of layers and neurons, the loss and activation functions is empirical,
and is motivated by some prior work [25]. We evaluate the network performance
using root mean square error (MSE) and mean relative error (MRE) metrics. The
training and evaluation are performed on a system running Ubuntu 18.04 with a
2.20GHz Intel Core i7-8750H CPU with 12 cores and 32 GB RAM. The network
training time, MSE and MRE are given in Tables 2 and 2 respectively.

Benchmark | Dims | Step size | Time bound | NN Training | MSE | MRE
---|---|---|---|---|---|---
|  |  | (sec) |  | Time (min) |  |
| Brussellator | 2 | 0.01 | 500 | 40.0 | 0.14 | 0.34
Buckling | 2 | 0.01 | 500 | 25.0 | 2.38 | 0.18
Lotka | 2 | 0.01 | 500 | 27.0 | 0.38 | 0.31
Jetengine | 2 | 0.02 | 300 | 35.0 | 0.086 | 0.63
Continuous | Vanderpol | 2 | 0.01 | 500 | 75.50 | 0.15 | 0.29
Nonlinear | Lacoperon | 2 | 0.1 | 500 | 66.0 | 0.12 | 0.33
Dynamics | Roesseler | 3 | 0.02 | 500 | 42.0 | 0.58 | 0.087
| Lorentz | 3 | 0.01 | 500 | 22.0 | 1.08 | 0.11
| Steam | 3 | 0.01 | 500 | 35.0 | 0.34 | 0.07
| C-Vanderpol | 4 | 0.01 | 500 | 70.0 | 0.18 | 0.15
| HybridOsc. | 2 | 0.01 | 500 | 80.0 | 0.35 | 0.11
Hybrid/ | SmoothOsc. | 2 | 0.01 | 500 | 38.5 | 0.40 | 0.096
NN Systems | Mountain Car | 2 | - | 100 | 12.5 | 0.015 | 0.79
| Quadrotor | 6 | 0.01 | 120 | 40.0 | 0.064 | 0.20

Table 1: Learning senstivity function. The set of benchmarks includes
nonlinear dynamical, hybrid and neural network based feedback control systems.
Time bound is number of steps for which the system simulation is computed.

Benchmark | Dims | Step size | Time bound | NN Training | MSE | MRE
---|---|---|---|---|---|---
|  |  | (sec) |  | Time (min) |  |
| Brussellator | 2 | 0.01 | 500 | 67.0 | 1.01 | 0.29
Buckling | 2 | 0.01 | 500 | 42.0 | 0.59 | 0.17
Lotka | 2 | 0.01 | 500 | 40.0 | 0.50 | 0.13
Jetengine | 2 | 0.01 | 300 | 34.0 | 1.002 | 0.26
Vanderpol | 2 | 0.01 | 500 | 45.50 | 0.23 | 0.23
Continuous | Lacoperon | 2 | 0.2 | 500 | 110.0 | 1.8 | 0.46
Nonlinear | Roesseler | 3 | 0.02 | 500 | 115.0 | 0.44 | 0.07
Dynamics | Lorentz | 3 | 0.01 | 500 | 67.0 | 0.48 | 0.08
| Steam | 3 | 0.01 | 500 | 58.0 | 0.13 | 0.057
| C-Vanderpol | 4 | 0.01 | 500 | 75.0 | 0.34 | 0.16
| HybridOsc. | 2 | 0.01 | 1000 | 77.0 | 0.31 | 0.077
Hybrid/ | SmoothOsc. | 2 | 0.01 | 1000 | 77.5 | 0.23 | 0.063
NN Systems | Mountain Car | 2 | - | 100 | 10.0 | 0.005 | 0.70
| Quadrotor | 6 | 0.01 | 120 | 25.0 | 0.0011 | 0.16

Table 2: Learning inverse senstivity function. Parameters and performance of
neural network tasked with learning inverse sensitivity function.

##  5 Space Space Exploration Using Neural Network Approximation

In this section, we present various applications in the domain of state space
exploration using the neural network approximations of sensitivity and inverse
sensitivity. The goal of state space exploration is to search for trajectories
that satisfy or violate a given specification. In this paper, we primarily
concern ourselves with a safety specification, that is, whether a specific
trajectory reaches a set of states labelled as _unsafe_. In order to search
for such trajectories, we present various empirical techniques that use both
forward and inverse sensitivity.

###  5.1 Reaching a Specified Destination Using Inverse Sensitivity
Approximation

In the course of state space exploration, after testing the behavior of the
system for a given set of test cases, the control designer might choose to
explore the behavior of a system that reaches a destination or approaches the
boundary condition for safe operation. Given a domain of operation $D$, we
assume that the designer provides a desired target state $z$ (with an error
threshold of $\epsilon$) that is reached by a trajectory at time $t$. Our goal
is to generate a trajectory $\xi$ such that $\xi(t)$ visits a state in the
$\epsilon$ neighborhood of the target $z$.

Our approach for generating the target trajectory is as follows. First, we
generate a random trajectory $\xi$ from the initial set $\Theta$, and compute
the difference vector of target state $z$ and $\xi(t)$. We now use the neural
network approximation of the inverse sensitivity function and estimate the
perturbtion required in the initial set such that the trajectory after time
$t$ goes through $z$. Since the neural network can only approximate the
inverse sensitivity function, the trajectory after the perturbation need not
visit $\epsilon$ neighborhood of the destination. However, we can repeat the
procedure until a threhold on the number of iterations is reached or the
$\epsilon$ threshold is satisfied. The pseudocode of this procedure, denoted
as reachTarget, is given in Algorithm 0.1

input : System simulator $\xi$, Time bound $T$, trained neural network
$\mathit{NN}_{\Phi^{-1}}$, time instance $t\leq T$, target state $z\in D$,
iteration count $I$, initial set $\theta$, and threshold $\epsilon$

output : State $x\in\theta$, $d_{r}$,
$d_{a}\mathrel{\stackrel{{\scriptstyle\scriptscriptstyle\Delta}}{{=}}}||\xi(x_{0},t)-z||_{2}$
such that $d_{a}\leq\epsilon$

1 $x\leftarrow x_{random}\in\theta$; $i\leftarrow 1$;

2 $x_{0}\leftarrow\xi(x,t)$; $d_{a}\leftarrow||x-z||_{2}$;

3 $d_{orig}\leftarrow d_{a}$; $d_{r}\leftarrow 1$;

4 while  _$(d_{a} >\epsilon$) & ($i\leq I)$ _ do

5 $v\leftarrow x_{0}-z$;

6 $v^{-1}\leftarrow\mathit{NN}_{\Phi^{-1}}(x_{0},v,t)$;

7 $x\leftarrow x+v^{-1}$; $x_{0}\leftarrow\xi(x,t)$;

8 $d_{a}\leftarrow||x_{0}-z||_{2}$; $i\leftarrow i+1$;

9

10 end while

11 $d_{r}\leftarrow\frac{d_{a}}{d_{orig}}$;

12 return $(x,d_{r},d_{a})$;

Algorithm 0.1 reachTarget. Finding an initial state from which the simulation
goes within $\epsilon$-neighborhood of destination $z$ at time $t$.

We initialize the variables in line 0.1. $\xi(x,.)$ is the simulation
generated by $\xi$ for the state $x$. The while loop (lines 0.1\- 0.1) runs
until either the threshold is reached or iteration count is exhausted.
$v^{-1}\mathrel{\stackrel{{\scriptstyle\scriptscriptstyle\Delta}}{{=}}}\Phi^{-1}(x_{0},v,t)$
is the inverse sensitivity function which is learned using a neural network
$\mathit{NN}_{{\Phi}^{-}1}$; $\left\lVert\cdot\right\rVert_{2}$ is the
$l2$-norm.

####  5.1.1 Evaluation of reachTarget on Standard Benchmarks

Benchmark | Dims | Iteration count = 1 | Iteration count = 5
---|---|---|---
|  | $d_{a}$ | $d_{r}$ | Time (ms) | $d_{a}$ | $d_{r}$ | Time (ms)
Brussellator | 2 | [0.19 - 1.87] | [0.23 - 0.74] | 11.38 | [0.003- 0.22] | [0.01 - 0.12] | 31.34
Buckling | 2 | [1.67 - 11.52 | [0.17 - 0.45] | 13.61 | [0.36- 2.09] | [0.06 - 0.31] | 34.51
Lotka | 2 | [0.08 - 0.24] | [0.21 - 0.45] | 12.38 | [0.02 - 0.07] | [0.09 - 0.22] | 34.28
Jetengine | 2 | [0.05 -0.20] | [0.19 - 0.28] | 15.96 | [0.0004 - 0.05] | [0.006 - 0.14] | 38.26
Vanderpol | 2 | [0.29 - 0.58] | [0.16 - 0.66] | 12.34 | [0.03 - 0.18] | [0.04 - 0.16] | 34.02
Lacoperon | 2 | [0.03 - 0.13] | [0.12 - 0.28] | 17.18 | [0.003 - 0.03] | [0.02 - 0.16] | 37.34
Roesseler | 3 | [0.72 - 2.02] | [0.20 - 0.34] | 16.08 | [0.21 - 0.63] | [0.06 - 0.14] | 38.26
Lorentz | 3 | [1.24 - 5.60] | [0.29 - 0.58] | 24.72 | [0.20 - 0.70] | [0.05 - 0.17] | 60.18
Steam | 3 | [1.59 - 5.21] | [0.31 - 0.67] | 8.68 | [0.41 - 1.8] | [0.08 - 0.30] | 69.80
C-Vanderpol | 4 | [0.87 - 1.72] | [0.34 - 0.60] | 17.44 | [0.20 - 0.40] | [0.07 - 0.18] | 44.86
HybridOsc. | 2 | [0.28 - 0.92] | [0.13 - 0.29] | 16.70 | [0.03 - 0.31] | [0.01 - 0.10] | 45.82
SmoothOsc. | 2 | [0.37 - 1.09] | [0.13- 0.23] | 52.22 | [0.04 - 0.42] | [0.02 - 0.18] | 136.72
Mountain Car | 2 | [0.004 - 0.24] | [0.08 - 0.22] | 138.90 | [0.0002 - 0.005] | [0.03 - 0.12] | 266.76
Quadrotor | 6 | [0.014 -1.09] | [0.10 - 0.67] | 284.96 | [0.004 - 0.04] | [0.02 - 0.13] | 668.78
Table 3: Evaluations. The results of reachTarget after iteration count 1 and
5. We compute average absolute distance $d_{a}$ and relative distance $d_{r}$
over 500 iterations of our algorithm for each benchmark. Additionally, a range
of values is obtained for $d_{a}$ and $d_{r}$ by performing the evaluation on
10 different targets.

We evaluate the performance of reachTarget algorithm as described in Section
5. We evaluate this by picking a random target state in the domain of interest
and let reachTarget generate a trajectory that goes through the neighborhood
of the target ($\epsilon=0.01$ or $0.001$) at a specified time $t$. Typically,
reachTarget executes the loop in lines 0.1-0.1 for 10 times before reaching
the target. In Table 3, we present the relative and absolute distance between
the target and the state reached by the trajectory generated by reachTarget
after one or five iterations of the main loop. The demonstration on benchmarks
Steam, Lotka and Coupled Vanderpol is shown in Figure 1. Experimental results
for the variations of the reachTarget are shown in in Figures 2 and 3

((a)) Steam

((b)) Coupled Vanderpol

((c)) Lotka

Figure 1: Illustration of reachTarget after iteration count 4. Iteration 0 is
the simulation from $x_{random}$. Subsequent 3 simulations are labeled as
Iteration 1, 2 and 3 respectively. As shown, with each iteration, the
proximity to the target increases.

((a)) Roesseler

((b)) Quadrotor

((c)) Lorentz

Figure 2: Illustration of fixed time instance based evaluation for single
target. Obtaining multiple simulations in the proximity of a given target at a
particular time instance.

((a)) Vanderpol

((b)) Brusselator

((c)) Hybrid Linear Oscillator

Figure 3: Illustration of fixed time instance based evaluation for multiple
targets. Obtaining multiple simulations in the proximity of each target at a
particular time instance.

We now discuss a few variations of our algorithm and approaches for evaluating
it.

  1. 1.

Uncertainty in time: The control designer might not be interested in reaching
the target at a precise time instance as long as it lies within a bounded
interval of time. In such cases, one can iterate reachTarget for every step in
this interval and generate a trajectory that approaches the target.

Consider the designer is interested in finding the maximum distance (or,
height) the car can go to on the left hill in Mountain Car. By providing an
ordered list of target states and a time interval, she can obtain the maximum
distance as well the time instance at which it achieves the maxima (Figure
4(a)). If there is no state in the given initial set from which the car can go
to a particular target, the approach, as a side effect, can also provide a
suitable initial candidate that takes the car as close as possible to that
target. In Quadrotor, one can easily find an initial configuration from which
the system can go to a particular location during a given time interval (Fig.
4(b)).

((a)) Mountain Car

((b)) Quadrotor

Figure 4: Illustration of time interval based evaluation. Finding an instance
in a given interval at which the distance to the target is minimum. The
initial set for Mountain Car is [-0.55, -0.45][0.0, 0.0]. Quadrotor has 5
discrete modes and each linear segment in 4(b) corresponds to a mode.

  2. 2.

Generalization: Based on our Mountain Car experiment, we observed that, for
the given initial set, the maximum distance the car can achieve on the left
hill is approx. 1.17. However, even after expanding the initial set to [-0.60,
-0.40][0.0, 0.0], our approach finds the maximum achievable distance (1.3019)
such that the car can still reach on the top of the right hill (shown in Fig.
5). This shows that our neural network is able to generalize the inverse
sensitivity over trajectories that go beyond the test cases considered during
the training process.

Figure 5: Illustration of Initial set expansion based evaluation. Computing
the maximum distance the car can achieve on the left hill after expanding the
initial set to [-0.60, -0.40][0.0, 0.0].

  3. 3.

Evaluating MRE for Random Targets: So far we have evaluated our technique with
respect to the target states that are reachable from given initial set. We
also perform evaluation of reachTarget by generating a random trajectory from
the domain and change its course at a provided time interval [25, 70] by a
randomly generated vector. The results are demonstrated in Figure 6.

((a)) Roesseler: Absolute error

((b)) Roesseler: Relative error

Figure 6: Illustration of Random vectors based evaluation. Plotting absolute
and relative errors for random vectors. The horizontal axis denotes the vector
norm. Each colored-plot corresponds to a particular time instance. The
absolute error increases linearly with vector size, whereas the relative error
appears to be converging. Interestingly, we observe similar behavior for each
system.

Discussion: It can be observed from Table 3 that our technique is capable of
achieving below 10% relative distance in almost all cases after 5 iterations.
That is, the trajectory generated by reachTarget algorithm after 5 iterations
is around 10% away from the target than the initial trajectory. This was the
case even for chaotic systems, hybrid sysems, and for control systems with
neural network components. While training the neural network might be time
taking process, the average time for generating new trajectories that approach
the target is very fast (less than a second for all cases). The high relative
distance in some cases might be due to high dimensionality or large distance
to the target which may be reduced further with more iterations.

###  5.2 Falsification of Safety Specification

One of the widely used methods for performing state space exploration are
falsification methods [43, 38]. Here, the specification is provided in some
temporal logic such as Signal or Metric Temporal Logic [35, 31]. The falsifier
then generates a set of test executions and computes the _robustness_ of
trajectory with respect to the specification. It then invokes heuristics such
as stochastic global optimization for discovering a trajectory that violates
the specification.

Given an unsafe set $U$, we provide a simple algorithm to falsify safety
specifications. We generate a fixed number ($m$) of random states in the
unsafe set $U$. Then, using the reachTarget sub-routine, generate trajectories
that reach a vicinity of the randomly generated states in $U$. We terminate
the procedure when we discover an execution that enters the unsafe set $U$.
For some standard benchmakrs, we compare the number of trajectories generated
by STL with the trajectories generated using inverse sensitivity in Fig. 7.
The box in each of the figures denotes the initial set and the red box
represents the unsafe set. Each of the points in the initial set represents a
sample trajectory generated by the falsification engine.

((a)) Brusselator

((b)) Brusselator

((c)) Lacoperon

((d)) Lacoperon

((e)) Simulated Annealing

((f)) Simulated Annealing

Figure 7: Illustration of falsification in S-Taliro and NeuralExplorer. When a
falsification tool fails to find a counterexample, NeuralExplorer may help
(using much less number of samples) the system designer by providing geometric
insight into the reason why the property is not satisfied.

((a)) Lorentz

((b)) Lorentz

((c)) Buckling

((d)) Buckling

((e)) Vanderpol

((f)) Vanderpol

Figure 8: Illustration of falsification in S-Taliro and NeuralExplorer. Some
scenarios where S-Taliro terminates with a falsification trajectory faster
than our approach.

Falsification using approximation of inverse sensitivity enjoys a few
advantages over other falsification methods. First, since our approach
approximates the inverse sensitivity, and we use the reachTarget sub-routine;
if the approximation is accurate to a certain degree, each subsequent
trajectory generated in reachTarget would make progress towards the
destination. Second, if the safety specification is changed slightly the
robustness of the trajectories with respect to new specification and the
internal representation for the stochastic optimization solver has to be
completely recomputed. However, since our trajectory generation does not rely
on computing the robustness for all the previously generated samples, our
algorithm is effective even when the safety specification is modified.

The third and crucial advantage of our approach lies when the falsification
tool does not yield a counterexample. In those cases, the typical
falsification tools cannot provide any geometric insight into the reason why
the property is not satisfied. However, using an approximation of inverse
sensitivity, the system designer can envison the required perturbation of the
reachable set in order to move the trajectory in a specific direction. This
geometric insight would be helpful in understanding why a specific trajectory
does not go into the unsafe set.

Considering these advantages, the results demonstrated in Fig. 7 should not be
surprising. We also would like to mention that these advantages come at the
computational price of training the neural networks to approximating the
inverse sensitivity. In addition to the examples shown above, we have included
some examples (see Fig. 8) where S-Taliro terminates with a falsification
trajectory faster than our approach. The reasons for such cases and methods to
improve falsification using NeuralExplore are a topic of future work.

###  5.3 Density Based Search Methods For State Space Exploration

One of the most commonly used technique for performing state space exploration
is generation of trajectories from a set of random points generated using an
apriori distribution. Based on the proximity of these trajectories to the
_unsafe_ set, this probability distribution can further be refined to obtain
trajectories that move closer to the unsafe set. However, one of the
computational bottlenecks for this is the generation of trajectories. Since
the numerical ODE solvers are sequential in nature, the refinement procedure
for probability distribution is hard to accelerate.

For this purpose, one can use the neural network approximation of sensitivity
to _predict_ many trajectories in an embarassingly parallel way. Here, a
specific set of initial states for the trajectories are generated using a pre-
determined distribution. Only a few of the corresponding trajectories for the
initial states are generated using numerical ODE solvers. These are called as
_anchor trajectories_. The remainder of trajectories are not generated, but
rather predicted using the neural network approximation of sensitivity and
anchor trajectories. That is, $\xi(x_{i},t)+\Phi_{NN}(x_{i},x_{j}-x_{i},t)$.
Additionally, the designer has the freedom to choose only a subset of the
initial states for only a specific time interval for prediction and refine the
probability distribution for generating new states. This would also allow us
to specifically focus on a time interval or a trajectory without generating
the prefix of it. A few examples of predictive trajectory generation for
performing _reachability analysis_ on Vanderpol oscillator and Jetengine are
provided in Fig. 9 and Fig. 10.

Figure 9: Illustration of sensitivity in reachability using Vanderpol. A
cluster of points is sampled in the neighborhood of a reference state. Actual
trajectories as well as predicted trajectories obtained by the neural network
which approximates sensitivity function are shown. Their corresponding
displacement from the reference trajectory at each time step are also
displayed in (d) and (e).

Figure 10: Illustration of sensitivity in reachability using Jetengine. Actual
trajectories as well as projected trajectories for a cluster of points sampled
in the neighborhood of a reference trajectory are shown. In addition, their
displacement from the reference trajectory at each time instance are
displayed.

####  5.3.1 Density Based Search for Falsification

Similar to the inverse sensitivity based falsification, one can use the
density based search space method for generating trajectories that reach a
destination and violate a safety specification. The forward density based
search procedure would work as follows. First, an anchor trajectory is
generated and time intervals of its trajectory that are closer to the unsafe
set are identified. Then a set of new initial states are generated according
to an apriori decided distribution. Instead of generating the trajectory from
these initial states, the _predicted trajectories_ using the anchor trajectory
and neural network approximation of sensitivity is generated specifically for
the time intervals of interest. Then, the initial state with the predicted
trajectory that is closest to the unsafe set is chosen and a new anchor
trajectory from the selected intial state is generated. This process of
generating anchor trajectory, new distribution of initial states, and moving
closer to the unsafe set is continued until you reach within the threhold that
is generated by the user.

Demonstration of this procedure for Vanderpol and Brusselator systems are
shown in Fig. 11 and 12 respectively. Notice that this approach gives an
underlying intuition about the geometric behavior of neighboring trajectories.
A similar method for density based estimation using inverse sensitivity
approximation can also be devised. Instead of sampling the initial set, the
density based method for inverse sensitivity generates random states around
the unsafe set to be reached and then generates a density map. An example of
such a density map generated is given in Fig. 13.

Figure 11: Illustration of sensitivity in falsification using Vanderpol. The
perturbation in the neighborhood of reference state are greedily chosen in an
iterative manner so as to minimize the distance to unsafe state. The projected
and actual displacements from reference trajectories are shown. The sampled
states are classified based on their euclidean distance to the unsafe state.

Figure 12: Illustration of sensitivity in falsification using Brusselator. For
a given unsafe state, a reference initial state and a set of perturbations in
its neighborhood are randomly chosen. The perturbation that minimizes the
distance to the unsafe state is greedily picked up and the new reference
trajectory is sampled from this perturbation. The procedure repeats until the
threshold is reached or the iteration count is exhausted. The states explored
in this iterative process are classified (in Fig. 12 and 12) based on their
respective distance to the unsafe state.

((a)) Brusselator: Distance profiles

((b)) Brusslator: Distance profiles for a different unsafe spec

((c)) Vanderpol: Distance profiles

((d)) Vanderpol: Distance profiles for a different time instance

Figure 13: Illustration of inverse sensitivity in falsification. Given a set
of states satisfying unsafe spec and a time instance, the routine, starting
from a reference trajectory, attempts to iteratively find a falsifying
trajectory. The initial states explored in the process are colored according
to their distance to unsafe states. These color densities help in identifying
regions in the initial set potentially useful for falsification. Notice the
difference in the distance profiles (color densities) as we change the time
instance (in Vanderpol) or select a difference unsafe spec (in Brusselator).

##  6 Conclusion and Future Work

We presented NeuralExplorer framework for state space exploration of closed
loop control systems using neural network. Our framework depends on computing
neural network approximations of two key properties of a dynamical system
called sensitivity and inverse sensitivity. We have demonstrated that for
standard benchmarks, these functions can be learned with less than $20\%$
relative error. We demonstrated that our method can not only be applied to
standard nonlinear dynamical systems but also for control systems with neural
network as feedback functions.

Using these approximations of sensitivity and inverse sensitivity, we
presented new ways to performing state space exploration. We also highlighted
the advantages of the falsification methods devised using the approximations.
Additionally, we demonstrated that our techniques give a geometric insight
into the behavior of the system and provide more intuitive information to the
user, unlike earlier black box methods. We believe that these techniques can
help the system designer in search of the desired executions. 111All the code
and examples for the evaluations performed in this section are available at
https://github.com/mag16154/NeuralExplorer

In future, we intend to extend this work to handle more generic systems such
as feedback systems with environmental inputs. We believe such a black-box
method for generating adversarial examples can be integrated into generative
adversarial training for training neural networks for control applications.

## References

  * [1] Benchmarks of continuous and hybrid systems. https://ths.rwth-aachen.de/research/projects/hypro/benchmarks-of-continuous-and-hybrid-systems/, benchmarks of continuous and hybrid systems
  * [2] Abbas, H., Fainekos, G.E.: Linear hybrid system falsification through local search. In: Automated Technology for Verification and Analysis, 9th International Symposium, ATVA 2011, Taipei, Taiwan, October 11-14, 2011. Proceedings. pp. 503–510 (2011)
  * [3] Althoff, M.: An introduction to cora 2015. In: Proc. of the Workshop on Applied Verification for Continuous and Hybrid Systems (2015)
  * [4] Alur, R., Dang, T., Ivančić, F.: Progress on reachability analysis of hybrid systems using predicate abstraction. In: Maler, O., Pnueli, A. (eds.) Hybrid Systems: Computation and Control (2003)
  * [5] Annpureddy, Y., Liu, C., Fainekos, G.E., Sankaranarayanan, S.: S-taliro: A tool for temporal logic falsification for hybrid systems (2011)
  * [6] Bak, S., Beg, O.A., Bogomolov, S., Johnson, T.T., Nguyen, L.V., Schilling, C.: Hybrid automata: from verification to implementation. International Journal on Software Tools for Technology Transfer (2019)
  * [7] Bak, S., Duggirala, P.S.: Hylaa: A tool for computing simulation-equivalent reachability for linear systems. In: Proceedings of the 20th International Conference on Hybrid Systems: Computation and Control (2017)
  * [8] Chen, R.T.Q., Rubanova, Y., Bettencourt, J., Duvenaud, D.: Neural ordinary differential equations (2018)
  * [9] Chen, X., Ábrahám, E., Sankaranarayanan, S.: Flow*: An analyzer for non-linear hybrid systems. In: Sharygina, N., Veith, H. (eds.) Computer Aided Verification (2013)
  * [10] Chollet, F., et al.: Keras. https://github.com/fchollet/keras (2015)
  * [11] Dang, T., Maler, O.: Reachability analysis via face lifting. In: Henzinger, T.A., Sastry, S. (eds.) Hybrid Systems: Computation and Control (1998)
  * [12] Deshmukh, J.V., Fainekos, G.E., Kapinski, J., Sankaranarayanan, S., Zutshi, A., Jin, X.: Beyond single shooting: Iterative approaches to falsification. In: American Control Conference, ACC (2015)
  * [13] Donzé, A.: Breach, A toolbox for verification and parameter synthesis of hybrid systems. In: Computer Aided Verification, 22nd International Conference, CAV 2010. Proceedings (2010)
  * [14] Donzé, A., Maler, O.: Systematic simulation using sensitivity analysis. In: International Workshop on Hybrid Systems: Computation and Control (2007)
  * [15] Donzé, A., Maler, O.: Robust satisfaction of temporal logic over real-valued signals. In: Formal Modeling and Analysis of Timed Systems - 8th International Conference, FORMATS 2010. Proceedings (2010)
  * [16] Duggirala, P.S., Mitra, S., Viswanathan, M.: Verification of annotated models from executions. In: Proceedings of the Eleventh ACM International Conference on Embedded Software (2013)
  * [17] Duggirala, P.S., Potok, M., Mitra, S., Viswanathan, M.: C2E2: a tool for verifying annotated hybrid systems. In: HSCC (2015)
  * [18] Dutta, S., Chen, X., Jha, S., Sankaranarayanan, S., Tiwari, A.: Sherlock - A tool for verification of neural network feedback systems: demo abstract. In: HSCC (2019)
  * [19] Fainekos, G.E., Pappas, G.J.: Robustness of temporal logic specifications for continuous-time signals. TCS (2009)
  * [20] Fan, C., Mitra, S.: Bounded verification with on-the-fly discrepancy computation. In: International Symposium on Automated Technology for Verification and Analysis. pp. 446–463 (2015)
  * [21] Fan, C., Qi, B., Mitra, S., Viswanathan, M.: Dryvr: Data-driven verification and compositional reasoning for automotive systems. In: CAV (1) (2017)
  * [22] Frehse, G., Le Guernic, C., Donzé, A., Cotton, S., Ray, R., Lebeltel, O., Ripado, R., Girard, A., Dang, T., Maler, O.: Spaceex: Scalable verification of hybrid systems. In: Proc. 23rd International Conference on Computer Aided Verification (CAV) (2011)
  * [23] Ghosh, S., Sadigh, D., Nuzzo, P., Raman, V., Donzé, A., Sangiovanni-Vincentelli, A.L., Sastry, S.S., Seshia, S.A.: Diagnosis and repair for synthesis from signal temporal logic specifications. In: Proceedings of the 19th International Conference on Hybrid Systems: Computation and Control (2016)
  * [24] Goyal, M., Duggirala, P.S.: On generating A variety of unsafe counterexamples for linear dynamical systems. In: ADHS (2018)
  * [25] Goyal, M., Duggirala, P.S.: Learning robustness of nonlinear systems using neural networks (2019), https://sites.google.com/view/dars2019/home
  * [26] Goyal, M., Duggirala, P.S.: Extracting counterexamples induced by safety violation in linear hybrid systems. Automatica 117, 109005 (2020). https://doi.org/https://doi.org/10.1016/j.automatica.2020.109005, http://www.sciencedirect.com/science/article/pii/S000510982030203X
  * [27] Huang, X., Kwiatkowska, M., Wang, S., Wu, M.: Safety verification of deep neural networks. In: CAV (1) (2017)
  * [28] Huang, Z., Mitra, S.: Proofs from simulations and modular annotations. In: Proceedings of the 17th international conference on Hybrid systems: computation and control (2014)
  * [29] Immler, F., Althoff, M., Chen, X., Fan, C., Frehse, G., Kochdumper, N., Li, Y., Mitra, S., Tomar, M.S., Zamani, M.: ARCH-COMP18 category report: Continuous and hybrid systems with nonlinear dynamics. In: ARCH@ADHS (2018)
  * [30] Ivanov, R., Weimer, J., Alur, R., Pappas, G.J., Lee, I.: Verisig: verifying safety properties of hybrid systems with neural network controllers (2018)
  * [31] Koymans, R.: Specifying real-time properties with metric temporal logic. Real-time systems (1990)
  * [32] Kyriakis, P., Deshmukh, J.V., Bogdan, P.: Specification mining and robust design under uncertainty: A stochastic temporal logic approach. ACM Trans. Embed. Comput. Syst. (2019)
  * [33] Lewis, F., Jagannathan, S., Yesildirak, A.: Neural network control of robot manipulators and non-linear systems. CRC Press (1998)
  * [34] Long, Z., Lu, Y., Ma, X., Dong, B.: PDE-net: Learning PDEs from data. In: Proceedings of the 35th International Conference on Machine Learning. pp. 3208–3216 (2018)
  * [35] Maler, O., Nickovic, D., Pnueli, A.: Checking temporal properties of discrete, timed and continuous behaviors. In: Pillars of computer science (2008)
  * [36] Miller, W., Werbos, P., Sutton, R.: Neural Networks for Control. A Bradford book, MIT Press (1995), https://books.google.com/books?id=prjMtIr\\_yT8C
  * [37] Moore, K.L.: Iterative learning control for deterministic systems. Springer Science & Business Media (2012)
  * [38] Nghiem, T., Sankaranarayanan, S., Fainekos, G., Ivancić, F., Gupta, A., Pappas, G.J.: Monte-carlo techniques for falsification of temporal properties of non-linear hybrid systems. In: Proceedings of the 13th ACM International Conference on Hybrid Systems: Computation and Control (HSCC 2010) (2010)
  * [39] Pan, S., Duraisamy, K.: Long-time predictive modeling of nonlinear dynamical systems using neural networks. Complexity 2018
  * [40] Pathak, J., Hunt, B., Girvan, M., Lu, Z., Ott, E.: Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach. Phys. Rev. Lett. (2018)
  * [41] Phan, D., Paoletti, N., Zhang, T., Grosu, R., Smolka, S.A., Stoller, S.D.: Neural state classification for hybrid systems. In: Proceedings of the Fifth International Workshop on Symbolic-Numeric Methods for Reasoning About CPS and IoT (2019)
  * [42] Raissi, M., Perdikaris, P., Karniadakis, G.E.: Multistep neural networks for data-driven discovery of nonlinear dynamical systems (2018)
  * [43] Sankaranarayanan, S., Fainekos, G.E.: Falsification of temporal properties of hybrid systems using the cross-entropy method. In: Hybrid Systems: Computation and Control, 2012 (2012)
  * [44] Sun, X., Khedr, H., Shoukry, Y.: Formal verification of neural network controlled autonomous systems. In: Proceedings of the 22Nd ACM International Conference on Hybrid Systems: Computation and Control (2019)
  * [45] Sun, Y., Huang, X., Kroening, D., Sharp, J., Hill, M., Ashmore, R.: Structural test coverage criteria for deep neural networks. ACM Trans. Embedded Comput. Syst. (2019)
  * [46] Tjeng, V., Xiao, K.Y., Tedrake, R.: Evaluating robustness of neural networks with mixed integer programming. In: ICLR (Poster). OpenReview.net (2019)
  * [47] Tran, H.D., Cai, F., Diego, M.L., Musau, P., Johnson, T.T., Koutsoukos, X.: Safety verification of cyber-physical systems with reinforcement learning control. ACM Trans. Embed. Comput. Syst. (2019)
  * [48] Zutshi, A., Deshmukh, J.V., Sankaranarayanan, S., Kapinski, J.: Multiple shooting, cegar-based falsification for hybrid systems. In: Proceedings of the 14th International Conference on Embedded Software (2014)