Sparsity-Sensitive Finite Abstraction

TL;DR

This paper introduces a sparsity-sensitive abstraction algorithm that significantly accelerates the process of creating finite models for high-dimensional systems with sparse interconnections, enabling applications previously infeasible due to computational constraints.

Contribution

The authors propose a simple modification to existing abstraction algorithms that exploits system sparsity, achieving substantial speed-ups while maintaining accuracy.

Findings

Speed-up in abstraction computation for high-dimensional systems

Successful synthesis of a safety controller for a 12D system

Efficient abstraction of a 51D vehicular traffic network

Abstract

Abstraction of a continuous-space model into a finite state and input dynamical model is a key step in formal controller synthesis tools. To date, these software tools have been limited to systems of modest size (typically $\leq$ 6 dimensions) because the abstraction procedure suffers from an exponential runtime with respect to the sum of state and input dimensions. We present a simple modification to the abstraction algorithm that dramatically reduces the computation time for systems exhibiting a sparse interconnection structure. This modified procedure recovers the same abstraction as the one computed by a brute force algorithm that disregards the sparsity. Examples highlight speed-ups from existing benchmarks in the literature, synthesis of a safety supervisory controller for a 12-dimensional and abstraction of a 51-dimensional vehicular traffic network.

Equations

$${\mathbf{x}}^{+}=\mathbf{f}(\mathbf{x},\mathbf{u})=\left(\begin{array}{c}{f}_1(\mathbf{x},\mathbf{u})\\ ⋮\\ f_n(\mathbf{x},\mathbf{u})\end{array}\right),$$

$$\mathsf{T}({\mathbf{d}}^x,{\mathbf{d}}^u,{\mathbf{d}}^{x^{+}})=\underset{i=1}{\overset{k}{\bigwedge }}\left(d_i^x==d_i^{x^{+}}\right),$$

$${\mathsf{O}^{\bm{f}}}=\left\{{[x^{+}]}\in[\mathcal{X}^{+}]\Big{|}{[x^{+}]}\cap\llbracket\bm{f}(\bm{a}^{x},\bm{a}^{u}),\bm{f}(\bm{b}^{x},\bm{b}^{u})\rrbracket\neq\emptyset\right\}.$$

$${\mathsf{O}^{\bm{f}}}=\left\{{[x^{+}]}\in[\mathcal{X}^{+}]\Big{|}{[x^{+}]}\cap\left(\bm{f}(\bm{c}^{x},\bm{c}^{u})+\llbracket-\bm{\beta},\bm{\beta}\rrbracket\right)\neq\emptyset\right\}$$

$$\left(\begin{array}{c}{x}_1^{+}\\ x_2^{+}\\ x_3^{+}\end{array}\right)=\left(\begin{array}{c}{f}_1(x_1,x_3,u_1,u_2)\\ f_2(x_2,x_3,u_1,u_2)\\ f_3(x_3,u_1,u_2)\end{array}\right),$$

$$\left(\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\end{array}\right)=\left(\begin{array}{c}{f}_1(x_3,u_1)\\ f_2(x_3,u_2)\\ f_3(u_1,u_2)\end{array}\right).$$

$$\left(\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\end{array}\right)=\left(\begin{array}{c}{u}_1\cos (\alpha +x_3)/\cos (\alpha )\\ u_1\sin (\alpha +x_3)/\cos (\alpha )\\ u_1\tan (u_2)\end{array}\right)$$

$${\mathcal{E}}_v^{\text{in}}$$

$${\mathcal{E}}_v^{\text{out}}$$

$${\mathcal{E}}^{\text{local}}(e)={\mathcal{E}}_{\tau (e)}^{\text{in}}\cup {\mathcal{E}}_{\tau (e)}^{\text{out}}\cup {\mathcal{E}}_{\eta (e)}^{\text{out}}.$$

$$\mathcal{X}=\prod _{e\in \mathcal{E}}[0,{\bar{x}}_e]\subset {\mathbb{R}}^n,$$

$$u_e=\{\begin{array}{ll}1& \text{if }\eta (e)\text{ actuates }e\text{ at time stepleavevmodenobreak t}\\ 0& \text{otherwise}\end{array}$$

$$\sum _{e\in {\mathcal{E}}_{\tau (k)}^{\text{in}}}{\alpha }_{ek}{u}_e=1.$$

$$f_e^{\text{out}}=u_e\min \left(x_e,c_e,\underset{k\in {\mathcal{E}}_{\eta (e)}^{\text{out}}}{\min }\left\{\frac{{\alpha }_{ek}}{{\beta }_{ek}}({\bar{x}}_k-x_k)\right\}\right).$$

$$f_e^{\text{in}}=\sum _{k\in {\mathcal{E}}_{\tau (e)}^{\text{in}}}{\beta }_{ke}{f}_k^{\text{out}}.$$

$$x_e^{+}=\min \left\{{\bar{x}}_e,x_e-f_e^{\text{out}}+f_e^{\text{in}}+d_e\right\},$$

$$\underset{\begin{array}{c}[x]\in [\mathcal{X}]\\ [u]\in [\mathcal{U}]\end{array}}{\bigvee }\left(\mathsf{V}([\mathcal{X}])==\mathsf{I}([x])\wedge \mathsf{V}([\mathcal{U}])==\mathsf{I}([u])\wedge \left(\underset{[x^{+}]\in {\mathsf{O}}^{\mathbf{f}}([x],[u])}{\bigvee }\mathsf{V}([{\mathcal{X}}^{+}])==\mathsf{I}([x^{+}])\right)\right).$$

$$\underset{\begin{array}{c}[x]\in [\mathcal{X}]\\ [u]\in [\mathcal{U}]\end{array}}{\bigvee }\left(\mathsf{V}([\mathcal{X}])==\mathsf{I}([x])\wedge \mathsf{V}([\mathcal{U}])==\mathsf{I}([u])\wedge \left(\underset{i=1}{\overset{n}{\bigwedge }}\left(\underset{[x^{+}]\in {\mathsf{O}}^{\mathbf{f}}([x],[u])}{\bigvee }\mathsf{V}([{\mathcal{X}}_i^{+}])==\mathsf{I}([x_i^{+}])\right)\right)\right).$$

$$\underset{i=1}{\overset{n}{\bigwedge }}\left(\underset{\begin{array}{c}[x]\in [\mathcal{X}]\\ [u]\in [\mathcal{U}]\end{array}}{\bigvee }\left(\mathsf{V}([\mathcal{X}])==\mathsf{I}([x])\wedge \mathsf{V}([\mathcal{U}])==\mathsf{I}([u])\wedge \left(\underset{[x^{+}]\in {\mathsf{O}}^{\mathbf{f}}([x],[u])}{\bigvee }\mathsf{V}([{\mathcal{X}}_i^{+}])==\mathsf{I}([x_i^{+}])\right)\right)\right).$$

$$\underset{i=1}{\overset{n}{\bigwedge }}\left(\underset{\begin{array}{c}[x]\in {\mathsf{P}}_i^{\mathbf{f}}([\mathcal{X}])\\ [u]\in {\mathsf{P}}_i^{\mathbf{f}}([\mathcal{U}])\end{array}}{\bigvee }\left(\mathsf{V}({\mathsf{P}}_i^{\mathbf{f}}([\mathcal{X}]))==\mathsf{I}([x])\wedge \mathsf{V}({\mathsf{P}}_i^{\mathbf{f}}([\mathcal{U}]))==\mathsf{I}([u])\wedge \left(\underset{[x_i^{+}]\in {\mathsf{O}}_i^{\mathbf{f}}([x],[u])}{\bigvee }\mathsf{V}([{\mathcal{X}}_i^{+}])==\mathsf{I}([x_i^{+}])\right)\right)\right),$$

Full text

I Introduction

A variety of software tools have been developed in the last decade to automatically synthesize controllers that enforce a formal specification. A common first step for symbolic controller synthesis requires abstracting the control system with a continuous state and input space into one with discrete states and inputs. However, this abstraction step has been cited as a key bottleneck in benchmarks for higher-dimensional systems [1][2]. Many abstraction procedures implemented in current tools traverse the state space in a brute force manner and suffer from an exponential runtime with respect to the sum of state and input dimensions. The majority of continuous-space systems however contain a coordinate structure, where the governing equation for one state variable depends only on a small subset of other variables.

We present a simple modification to the brute force algorithm that takes advantage of sparsely interconnected structures. The modified abstraction algorithm judiciously looks at lower dimensional sub-spaces while ignoring those sub-spaces that do not affect the corresponding coordinate. The proposed procedure returns the same abstraction of the control system as the brute force algorithm while exploiting the structure. Our method requires only a linear number of operations with respect to the number of state dimensions, but this number grows exponentially with respect to a new sparsity parameter. For control systems that exhibit a sparse interconnection structure, this introduced parameter is small compared to the sum of state and input dimensions.

We provide two examples to highlight the practical speed-ups that can be attained with this modified abstraction procedure. The first is a bicycle example that is a common benchmark among controller synthesis tools. The second example is an arterial traffic corridor that is actuated with signalized intersections. These benchmarks confirm that the abstraction time and size of the abstraction encoding grow roughly linearly with the traffic network length.

I-A Literature Review

Abstraction bottlenecks have motivated research into compositional synthesis [3][4] and construction of abstractions [5][6]. These techniques leverage a system decomposition and reason over interconnecting variables among sub-systems. Therefore, our abstraction procedure can complement compositional techniques by also taking advantage of sparsity within a component.

Similar to our sparsity-sensitive finite abstraction algorithm, multiple procedures have appeared in the model checking and verification literature over the last decades. Finite abstractions of Markov processes have been constructed by leveraging independence assumptions on the underlying process noise [7]. In an application to symbolic model checking of sequential Boolean circuits, the abstraction encoding of a monolithic circuit has been decomposed implicitly in abstractions for each logical component independently [8]. Whereas the literature above focuses on applications to model checking and verification, in this paper we construct finite abstractions with controller synthesis as end goal.

I-B Software Tools

In recent years, the breadth of abstraction and corresponding controller synthesis software tools has grown to include a greater variety of dynamical systems and classes of specifications. Examples of recent tools include CoSyMA [9], PESSOA [10], and SCOTS [11].

At their core, each of these tools relies on a system relation between the original and the abstract control system. Efficient algorithms and data structures from the formal methods community are then leveraged to enforce the specifications for the abstracted system. The obtained controller of the finite input and state representation is then refined to the original system with continuous input and state spaces. The correctness of this controller design is guaranteed by relating the finite abstraction with the original control system via a system relation [12][13]. In this paper, we utilize feedback refinement relations to construct controllers which require quantized state information only [13]. However, the content in our paper can also accommodate other abstraction formalisms.

II Preliminaries

Given two sets $\mathcal{A}$ and $\mathcal{B}$, let $|\mathcal{A}|$, $2^{\mathcal{A}}$ and $\mathcal{A}\times\mathcal{B}$ respectively represent $\mathcal{A}$’s cardinality, power set (set of all subsets), and Cartesian product with $\mathcal{B}$.

Let $\mathbb{R}^{n}$ represent a $n$-dimensional Euclidean space. Given a vector element $\bm{x}\in\mathbb{R}^{n}$, $x_{i}$ represents its $i^{\text{th}}$ component for $i\in\{1,\ldots,n\}$. When convenient, we explicitly express $\bm{x}$’s coordinates as an ordered tuple $(x_{1},\ldots,x_{n})$. Given a pair of vectors $\bm{x},\bm{y}\in\mathbb{R}^{n}$, let $\llbracket\bm{x},\bm{y}\rrbracket$ denote the hyper-rectangle $[x_{1},y_{1}]\times\ldots\times[x_{n},y_{n}]$ and $\bm{z}+\llbracket\bm{x},\bm{y}\rrbracket$ is the same except with an offset $\bm{z}\in\mathbb{R}^{n}$. A set $[a]\subset\mathbb{R}^{n}$ is a box if there exists $\bm{x},\bm{y}\in\mathbb{R}^{n}$ such that $[a]=\llbracket\bm{x},\bm{y}\rrbracket$. The vectors $\bm{x}$ and $\bm{y}$ are the lower and upper bound of $[a]$, respectively.

A set $[\mathcal{A}]=\prod_{i=1}^{n}[\mathcal{A}_{i}]$ is a box partition of $\mathcal{A}\subset\mathbb{R}^{n}$, if each $[\mathcal{A}_{i}]$ is an interval partition, i.e., $[\mathcal{A}]$ inherits a coordinate structure and the union of all $[a]\in[\mathcal{A}]$ forms a cover of $\mathcal{A}$. A non-empty intersection between two elements in $[\mathcal{A}]$ may only occur if they are adjacent boxes, and such an intersection has zero measure. Therefore, a box partition is a generalization of uniform space partitions, where the boxes do not have to be all congruent to each other.

Every box $[a]\in[\mathcal{A}]$ can be uniquely identified with a vector of indices $\bm{d}^{a}\in\mathbb{Z}_{>0}^{n}$, where $d^{a}_{i}$ is the corresponding index of the $i^{\text{th}}$ coordinate dimension. Let $\mathsf{I}\colon[\mathcal{A}]\to\mathbb{Z}_{>0}^{n}$ be an invertible map from an element of a box partition to its corresponding vector of indices.

The set of Booleans $\mathbb{B}=\{\top,\bot\}$ consists of two elements, Boolean true $\top$ and false $\bot$. Let $\neg$, $\land$ and $\lor$ denote the logical negation, conjunction and disjunction, respectively.

III Construction of Abstractions

In this section, we present the brute force algorithm which constructs an abstraction by traversing the whole state and input space. Subsequently, an overview of existing methods to compute a reachable set over-approximation from a continuous-space control system is given.

III-A Control System

A discrete-time, continuous-space control system is given by a system of update equations

[TABLE]

where $\bm{x}^{+}\in\mathbb{R}^{n}$ represents the state at the next time step. The function $\bm{f}\colon\mathcal{X}\times\mathcal{U}\to\mathcal{X}^{+}$ is defined over a state space $\mathcal{X}\subset\mathbb{R}^{n}$ and input space $\mathcal{U}\subset\mathbb{R}^{m}$ and maps to the space of next states $\mathcal{X}^{+}\equiv\mathcal{X}$, which are subsets of a Euclidean space of appropriate dimensions.

Assumption 1**.**

The state space $\mathcal{X}$ and input space $\mathcal{U}$ are bounded hyper-rectangles.

By construction, $[\mathcal{X}]$, $[\mathcal{U}]$, $[\mathcal{X}^{+}]$ are finite and all corresponding box elements ${[x]}$, ${[u]}$, ${[x^{+}]}$ can be represented as index vectors. Moreover, all three box partitions inherit a coordinate structure, for instance $[\mathcal{X}^{+}]=\prod_{i=1}^{n}[\mathcal{X}^{+}_{i}]$.

III-B Finite Abstractions

Definition 1**.**

A finite abstraction of a control system is a quadruple $([\mathcal{X}],[\mathcal{U}],[\mathcal{X}^{+}],\mathsf{T})$, where $[\mathcal{X}]$, $[\mathcal{U}]$, $[\mathcal{X}^{+}]$ are associated box partitions of $\mathcal{X}$, $\mathcal{U}$, $\mathcal{X}^{+}$, and $\mathsf{T}\colon\mathsf{I}([\mathcal{X}])\times\mathsf{I}([\mathcal{U}])\times\mathsf{I}([\mathcal{X}^{+}])\to\mathbb{B}$ is a Boolean function that encodes valid transitions.

A transition from ${\bm{d}^{x}}\in\mathsf{I}([\mathcal{X}])$ to $\bm{d}^{x^{+}}\in\mathsf{I}([\mathcal{X}^{+}])$ under input $\bm{d}^{u}\in\mathsf{I}([\mathcal{U}])$ is valid if and only if $\mathsf{T}({\bm{d}^{x}},\bm{d}^{u},\bm{d}^{x^{+}})=\top$. Thus, the pre-image $\mathsf{T}^{-1}(\top)\subseteq\mathsf{I}([\mathcal{X}])\times\mathsf{I}([\mathcal{U}])\times\mathsf{I}([\mathcal{X}^{+}])$ is a set that encodes all valid transitions. Although the information stored in the transition function $\mathsf{T}$ and the set of valid transitions $\mathsf{T}^{-1}(\top)$ are equivalent, we opt to use the symbolic function representation because it can easily be compressed, as shown by the following example.

Example 1** (Valid Transition Representations).**

Consider the transition function $\mathsf{T}({\bm{d}^{x}},\bm{d}^{u},\bm{d}^{x^{+}})=\top$ if and only if ${\bm{d}^{x}}=\bm{d}^{x^{+}}$, which encodes stationary dynamics and no effect from control inputs. If $k$ and $l$ are the total number of bits used to represent an element of $[\mathcal{X}]$ and $[\mathcal{U}]$, respectively, then enumerating all elements of $\mathsf{T}^{-1}(\top)$ in a sparse matrix lookup table requires $2^{k+l}$ non-zero entries. On the other hand, $\mathsf{T}$ can also be represented as Boolean function with $k$ clauses

[TABLE]

where the $i^{\text{th}}$ clause encodes equivalence between the $i^{\text{th}}$ index of ${\bm{d}^{x}}$ and $\bm{d}^{x^{+}}$ for $i\in\{1,\ldots,k\}$.

Section VI presents further examples where Boolean function representations allow us to construct abstractions with more than ${10}^{60}$ valid transitions with as little as hundreds of kilobytes.

III-C Brute Force Abstraction

A core subroutine in each of the controller synthesis tools presented in Section I-B is the computation of a finite abstraction from the continuous-space dynamics. Discrete states and inputs correspond with non-empty boxes in the continuous space. Transitions among discrete states ${[x]}$ under inputs ${[u]}$ are computed by first executing 1 from a representative state $\bm{x}\in{[x]}$ and input $\bm{u}\in{[u]}$, then adding additional transitions to account for the quantization error. Since the computation of exact reachable sets of boxes ${[x]}$ under input ${[u]}$ is too difficult, the finite abstraction over-approximates the behavior of the original continuous-space control system.

Let the reachable set over-approximation function of the control system be given by ${\mathsf{O}^{\bm{f}}}\colon[\mathcal{X}]\times[\mathcal{U}]\to 2^{[\mathcal{X}^{+}]}$ . It maps a state and input box to a set of next state boxes. We will provide further details about efficient over-approximation techniques for different classes of dynamical systems in Section III-D.

We assume that an over-approximation function ${\mathsf{O}^{\bm{f}}}$ is given to the brute force abstraction algorithm, which computes the transition function $\mathsf{T}$ iteratively. Since $\mathsf{T}$ is a function with inputs, we introduce as placeholders some free variables $\mathsf{V}([\mathcal{X}])$, $\mathsf{V}([\mathcal{U}])$ and $\mathsf{V}([\mathcal{X}^{+}])$. For instance, let $\mathsf{V}([\mathcal{X}])$ be a map from a state box partition $[\mathcal{X}]$ to a vector of corresponding free variables. The brute force abstraction procedure is given by Algorithm 1, where the double equal sign denotes component-wise equality.

This algorithm iterates over all elements of $[\mathcal{X}]$ and $[\mathcal{U}]$, which can be seen as a nested for loop of depth equal to the sum of the state and input dimensions. Thus, the number of discrete states and inputs grows exponentially with each additional dimension, i.e., obtaining a finite abstraction becomes an intractable problem for higher-dimensional systems.

III-D Reachable Set Over-approximations

In Algorithm 1, we assumed that a method ${\mathsf{O}^{\bm{f}}}$ to over-approximate reachable sets is given. There exist multiple techniques to over-approximate the system dynamics, tailored for different classes of systems. In the following, we briefly present a few of them.

Let the two boxes ${[x]}=\llbracket\bm{a}^{x},\bm{b}^{x}\rrbracket$ and ${[u]}=\llbracket\bm{a}^{u},\bm{b}^{u}\rrbracket$ be the inputs of ${\mathsf{O}^{\bm{f}}}$. An efficient estimate of the reachable set of nonlinear, monotone dynamical systems is presented in [14]. The authors consider discrete-time systems that are monotone with respect to the non-negative convex cone $\mathbb{R}_{\geq 0}^{n}$. The over-approximation function for this class is

[TABLE]

Monotone systems are a subclass of mixed monotone systems, for which there exists a generalized procedure that accounts for increasing or decreasing components of $\bm{f}$ [15].

An over-approximation method for systems with perturbed measurements is presented in [13]. A local error bound $\bm{\beta}$ is introduced, which bounds the uncertainty in each state coordinate independently based on the current ${[x]}$ and ${[u]}$. Let the centers of the state and input box be $\bm{c}^{x}=\frac{1}{2}(\bm{b}^{x}-\bm{a}^{x})$ and $\bm{c}^{u}=\frac{1}{2}(\bm{b}^{u}-\bm{a}^{u})$, respectively. Then

[TABLE]

is the corresponding over-approximation function for this class of dynamical systems.

IV Sparsity-Sensitive Abstraction for Discrete-Time Systems

In this section, we propose modifications of Algorithm 1 for discrete-time, sparsely interconnected control systems. The brute force algorithm is suited for systems where a change in each coordinate affects the next states in all other coordinates. However, when a coordinate update is independent of many other coordinates, then redundant computations are contained in Algorithm 1’s nested for loops. Thus, we can compute finite representations of systems more efficiently by exploiting this sparse interconnection structure.

Example 2** (Discrete-Time Redundant Iterations).**

Consider the following system consisting of three state and two input variables

[TABLE]

where each coordinate update depends on at most four instead of all five possible variables. Then Algorithm 1’s for loop iteration over independent variables provides no further information on the next state’s component. Thus, instead of dealing with three 5-dimensional update equations, we can consider only two 4-dimensional and one 3-dimensional problems.

IV-A Dependency Graph

Before we propose modifications to Algorithm 1, we illustrate system dependencies as a directed graph, called dependency graph subsequently.

Each vertex of the dependency graph $v_{i}$ for $i\in\{1,\ldots,n+m\}$ represents a state or input variable. If the state update equation associated with vertex $v_{i}$ depends on $v_{j}$, a directed edge from $v_{j}$ to $v_{i}$ is added. Consequently, the indegree of a vertex is equal to the number of directly dependent states and inputs. The dependency graph for the control system in Example 2 is illustrated in Fig. 1.

If the indegree of a next state vertex is smaller than the sum of state and input variables, i.e. if the next state variable does not depend on all possible variables, it is advantageous to exploit the sparse interconnection structure.

IV-B Coordinate-Wise Decomposition

In order to exploit sparsely interconnected systems, we first define a vector of projections $\mathsf{P}^{\bm{f}}$ associated with $\bm{f}$. Let the $i^{\text{th}}$ component $\mathsf{P}^{\bm{f}}_{i}$ of $\mathsf{P}^{\bm{f}}$ be a projection of a discrete box or corresponding box partition onto the lower dimensional domain of the $f_{i}$ component for $i\in\{1,\ldots,n\}$. For instance, the associated box partition projections of the first component of $\bm{f}$ in Example 2 are given by $\mathsf{P}^{\bm{f}}_{1}([\mathcal{X}])=[\mathcal{X}_{1}]\times[\mathcal{X}_{3}]$ and $\mathsf{P}^{\bm{f}}_{1}([\mathcal{U}])=[\mathcal{U}_{1}]\times[\mathcal{U}_{2}]$, respectively.

Before we present our algorithm modifications, we need to make the following assumption, which is fulfilled for all reachable set over-approximation functions in Section III-D.

Assumption 2**.**

The over-approximation function ${\mathsf{O}^{\bm{f}}}$ can be decomposed coordinate-wise. This means that ${\mathsf{O}^{\bm{f}}}({[x]},{[u]})=\prod_{i=1}^{n}{\mathsf{O}^{\bm{f}}_{i}}\left(\mathsf{P}^{\bm{f}}_{i}({[x]}),\mathsf{P}^{\bm{f}}_{i}({[u]})\right)$, where the over-approximation method associated with the $i^{\text{th}}$ coordinate dimension is given by ${\mathsf{O}^{\bm{f}}_{i}}\colon\mathsf{P}^{\bm{f}}_{i}([\mathcal{X}])\times\mathsf{P}^{\bm{f}}_{i}([\mathcal{U}])\to 2^{[\mathcal{X}^{+}_{i}]}$ for $i\in\{1,\ldots,n\}$.

IV-C Sparsity-Sensitive Abstraction

We propose Algorithm 2 below, which takes advantage of dependencies among state, input and next state variables by executing a lower dimensional brute force abstraction for each coordinate independently.

Claim 1**.**

(Equivalence between Abstraction Algorithms) Algorithm 1 and Algorithm 2 return the same abstraction.

A proof is provided in the appendix.

Let g=\max\big{(}|[\mathcal{X}_{1}]|,\ldots,|[\mathcal{X}_{n}]|,|[\mathcal{U}_{1}]|,\ldots,|[\mathcal{U}_{m}]|\big{)} be the maximum number of intervals along all dimensions. Furthermore, let the new sparsity parameter $h\leq n+m$ be the maximum indegree of all vertices in the dependency graph. The sparsity parameter $h$ can also be computed based on the adjacency matrix $A\in\mathbb{R}^{(n+m)\times(n+m)}$ of the dependency graph. An element $A_{ij}$ of the adjacency matrix is one if there exists an edge from vertex $v_{i}$ to $v_{j}$ and zero otherwise. The indegree of vertex $v_{i}$ is simply given by the sum of all entries in column $i$. For instance, the maximum indegree of the graph in Fig. 1 is $4$. Then, the brute force Algorithm 1’s $O(g^{n+m})$ complexity is reduced to $O(ng^{h})$. While Algorithm 2 is still exponential, the value of the exponent is smaller than the sum of states and inputs for control systems with sparse dependency structures, i.e. $h<n+m$.

V Sparsity-Sensitive Abstraction for Continuous-Time Systems

For continuous-time systems, the dynamics are given by a system of differential equations. An abstraction, i.e. a finite input and state representation, is associated with the corresponding sampled system with time quantization parameter $\tau\in\mathbb{R}_{>0}$. The system inputs are assumed to be constant for each sampling interval of length $\tau$. In Section IV, we have proposed Algorithm 2 as a modification of Algorithm 1, which is valid for discrete-time systems. In discrete-time, a state at the next time step depends solely on the variables appearing in its update equation. This idea can be generalized to continuous-time systems, where a change in one variable can have an indirect effect on the time evolution of another one in infinitesimal amount of time. When the dynamics are given in continuous time, the notion of dependencies requires a slight modification, which we present in the following.

The dependency graph for the continuous-time case, similar to the discrete one in Section IV-A, consists of vertices which represent a state or input variable. If the differential equation associated with vertex $v_{i}$ depends on $v_{j}$, a directed edge from $v_{j}$ to $v_{i}$ is added, similar to the update equation in discrete-time. In contrast to the state update equations, the time evolution of a state in continuous-time always depends on the initial value. Therefore, loops for all vertices which correspond to a state need to be added additionally in the dependency graph.

Example 3** (Continuous-Time Dependency Graph).**

Consider the following system of ordinary differential equations

[TABLE]

The associated dependency graph can be obtained from the graph in Fig. 1 by deleting the directed edges $(u_{1},x_{2})$ and $(u_{2},x_{1})$.

The number of directly dependent variables of a vertex of the dependency graph is equal to the indegree. Furthermore, the maximum indegree over all vertices corresponds to the sparsity parameter $h$ in the discrete-time case. Since there exist solely direct dependencies, we only need to consider walks of length one. However, for continuous-time systems there can be indirect influences that are not immediately reflected in the system equations. For instance, consider a double integrator system where the position state depends indirectly on the incoming acceleration state via the velocity state. In general, state vertex $v_{i}$ depends on any vertex $v_{j}$ if there exists a walk from $v_{j}$ to $v_{i}$. The number of dependencies for each vertex can be computed via a breadth-first search for all state vertices. Alternatively, we can compute these number for all vertices based on the adjacency matrix $A$. The element $\left(A^{l}\right)_{ij}$ corresponds to the number of walks of length $l$ from vertex $v_{i}$ to $v_{j}$. Therefore, the number of vertices from which $v_{i}$ can be reached is equal to the number of non-zero elements of the $i^{\text{th}}$ column of $\tilde{A}=\sum_{j=0}^{n}A^{j}$. We need to sum up to $n$, since the longest path from any vertex to another one is of length $n$. A path is defined as a walk which does not traverse any vertex more than once. Then, the sparsity parameter $h$ for continuous-time systems is given by the maximum number of non-zero elements of all columns of $\tilde{A}$.

VI Examples

We illustrate our sparsity-sensitive algorithm on two examples taken from the literature. All computations are run on a single thread of an Intel® Core™ i$7$-$5500$U ($3.00\text{\,}\frac{\mathrm{GHz}}{}$) with $8\text{\,}\frac{\mathrm{GB}}{}$ RAM. We use a modified version of the open-source tool SCOTS [11], which represents transition functions, specifications and controllers as single binary decision diagrams (BDD) [16].

VI-A Bicycle Model

As a first example, we want to obtain an abstraction for the continuous-time, continuous-space bicycle model with corresponding parameters presented in [13]. The position of the modeled vehicle in the $2$-dimensional plane and the orientation is given by the states $(x_{1},x_{2})$ and $x_{3}$, respectively. The two controllable inputs $u_{1}$ and $u_{2}$ are the velocity and steering angle. The vehicle dynamics are given by the following system of ordinary differential equations (ODEs)

[TABLE]

with $\alpha=\arctan(\tan(u_{2})/2)$.

We can solve 7 analytically to obtain a discrete-time control system with sampling time $\tau$, because the two inputs $u_{1}$ and $u_{2}$ are assumed to be constant for each sampling interval. The corresponding dependency graph of the exact, discrete-time representation of the sampled vehicle system is illustrated in Fig. 1. Since the software tools presented in Section I-B solve 7 numerically to compute an abstraction, we can also use the same update rules of the employed ODE solver to get a discrete-time system. Consequently, we obtain the identical transition function for both algorithms and can determine the computational speed-up.

The construction of the abstraction takes $2.3\text{\,}\frac{\mathrm{s}}{}$ based on our proposed Algorithm 2. In contrast, it requires $104.0\text{\,}\frac{\mathrm{s}}{}$ to compute the identical transition function based on the brute force procedure implemented in [11], while using the same over-approximation function 4. Therefore, we are able to speed-up the computation time by a factor of $45$, while the memory requirements are identical. Due to its partial implementation in MATLAB® for the nonlinear case, PESSOA [10] needs more than ${10}^{5}\text{\,}\frac{\mathrm{s}}{}$ on a $2.40\text{\,}\frac{\mathrm{GHz}}{}$ machine to compute the abstraction [13].

VI-B Traffic Model

The traffic model used in this example is taken from the modified, discrete-time cell transmission model presented in [3]. The traffic network can be represented as directed graph with a set $\mathcal{E}$ of edges or links and a set $\mathcal{V}$ of vertices or junctions. Every edge $e$ is a directed arc from tail vertex $\tau(e)\in\{\mathcal{V}\cup\emptyset\}$ to head vertex $\eta(e)\in\mathcal{V}$. By convention, an edge with $\tau(e)=\emptyset$ serves as entry point to the network and edges, which outflow the network, are not explicitly modeled. Let

[TABLE]

denote the incoming, respectively, outgoing edges of vertex $v$. The set of local edges to edge $e$ is given by

[TABLE]

Each edge $e$ $\in\mathcal{E}$ has an associated link occupancy $x_{e}\in[0,\bar{x}_{e}]$ at time step $t\in\mathbb{Z}_{\geq 0}$, where the constant $\bar{x}_{e}$ denotes the maximum vehicular capacity of edge $e$. Then, the continuous, $n$-dimensional state space of the traffic model is

[TABLE]

where $n=|\mathcal{E}|$. An edge $e$ is actuated if the traffic signal at vertex $\eta(e)$ permits outward flow from $e$ to $\mathcal{E}^{\text{out}}_{\eta(e)}$. The controllable input $u_{e}$ associated with edge $e$ is given by

[TABLE]

The constant split ratio $\beta_{ek}\in[0,1]$ denotes the fraction of vehicles flowing from edge $e$ to edge $k\in\mathcal{E}^{\text{out}}_{\eta(e)}$. However, edge $e$ can only send vehicles to $k$ if $k$ offers enough capacity. For this purpose, the constant supply ratio $\alpha_{ek}\in[0,1]$ denotes the capacity fraction of $k$ dedicated to $e$. Since the supply is split among actuated incoming edges, it follows that for all $k\in\mathcal{E}$ and time steps $t$

[TABLE]

In the following, we introduce the system dynamics based on the modified cell transmission model, which exhibits FIFO property. As we will see, the outflow and state update function of edge $e$ depends only on the corresponding states of edges in $\mathcal{E}^{\text{out}}_{\eta(e)}$ and $\mathcal{E}^{\text{local}}(e)$, respectively. Thus, the traffic network exhibits a locality property, which we can exploit to compute abstractions.

The flow out of edge $e\in\mathcal{E}$ is given by

[TABLE]

The minimization in 14 implies that the flow of exiting vehicles of edge $e$ cannot exceed the current link occupancy, the constant saturation flow $c_{e}$ and the resulting supply offered by edges in $\mathcal{E}^{\text{out}}_{\eta(e)}$. The saturation flow $c_{e}$ is a model parameter, which is determined by, e.g., number of lanes and speed limits. The inflow of edge $e$ is denoted by

[TABLE]

Based on the principle of mass conservation, the discrete-time update equation of state $x_{e}$ associated with edge $e$ is

[TABLE]

where $d_{e}\in[0,\bar{d}_{e}]$ represents a time-dependent, exogenous flow entering edge $e$ and the constant parameter $\bar{d}_{e}$ is an upper bound. It is clear based on 14 to 16, that the state update equation of edge $e$ depends only on the current states of edges in $\mathcal{E}^{\text{local}}(e)$. Therefore, we can reduce the computation time of the abstraction drastically by exploiting this sparse interconnection structure of the traffic network.

Arterial Traffic Corridor

We consider an easily extensible network to show the applicability of our proposed algorithm modifications to more than three dimensions. The considered $n$-dimensional traffic network is depicted in Fig. 2 on the left side of the big plus sign.

The dashed arrows represent road links, which exit the network and are not explicitly modeled. Since we want to evaluate the performance of our approach with respect to the dimensions of the state and input space, we consecutively add a new block to the existing, $n$-dimensional network as shown in Fig. 2. Each newly added block consists of three edges and one vertex, i.e., three road links and one signalized intersection. Thus, the dimension of the state space in 11 is increased by $3$. The input space grid doubles, since we assume that every vertex actuates either the incoming arterial road edge or the two intersecting collector edges at time step $t$.

The parameters of the $n$-dimensional network in Fig. 2 are as follows. The supply ratio $\alpha_{e_{i},e_{j}}$ is equal $1$ and the split ratio $\beta_{e_{i},e_{j}}$ is $0.8$ if $i+3=j$ for $j\in\{4,7,\ldots,n-2\}$. Otherwise, both ratios are equal $0.5$. The maximum vehicular capacity $\bar{x}_{e}$ is $10$ for all edges, and every dimension is quantized in $10$ equally spaced intervals. Therefore, the considered $n$-dimensional network has $10^{n}$ state space grid points. For all arterial road edges, the upper disturbance bound $\bar{d}_{e}$ is [math] and the saturation flow $c_{e}$ is $6$. For all intersecting collector edges, the bound is $1$ and the corresponding flow is $2$.

Results

Since the considered $n$-dimensional traffic network in Fig. 2 is a monotone system, we can use 3 as reachable set over-approximation function. The time needed to compute the finite abstraction with respect to the dimension of the state space is illustrated in Fig. 3.

In contrast to the exponential growth of the brute force algorithm, our proposed Algorithm 2 performs roughly linearly with the network length, although the total number of transitions grows exponentially. For instance, there are $3.3\text{\times}{10}^{62}$ valid transitions in the $51$-dimensional traffic network.

The file size of the resulting transition function BDD with respect to the dimension of the state space is shown in Fig. 4.

Due to the memory-efficient data structure, the size of the binary decision diagram which encodes the transition function grows also roughly linearly.

Based on the computed transition function, we can use the fixed point algorithm implemented in [11] to synthesize a memoryless, safety supervisory controller for the vehicular traffic network. For instance, assume the occupancy of every edge has to be smaller than $\frac{4}{5}$ of the maximum vehicular capacity for all time steps, i.e. in our example less than $8$. After representing this safety specification as BDD, we obtain a controller for the $12$-dimensional network after $46$ iterations of the fixed point algorithm. The synthesis takes $2.5\text{\,}\frac{\mathrm{min}}{}$ and $0.9\text{\,}\frac{\mathrm{MB}}{}$ are needed to store the obtained controller as binary decision diagram. This BDD represents all allowed inputs for each discrete state such that the safety specification is fulfilled. In this example, there exist more than $2.1\text{\times}{10}^{011}$ allowed state-input combinations.

Instead of using the memory-efficient BDD data structure, it is also possible to store the controller as sparse matrix lookup table by explicitly enumerating the state and input space. However, based on the given dimensions $44$ bits are needed to store a single allowed state-input combination. This results in $1.2\text{\,}\frac{\mathrm{TB}}{}$ of memory usage in contrast to $0.9\text{\,}\frac{\mathrm{MB}}{}$ based on the BDD approach.

In contrast to the $12$-dimensional network, it takes $4.3\text{\,}\frac{\mathrm{h}}{}$ until convergence in the $15$-dimensional case for the same safety specification. This increase of the controller synthesis computation time has multiple reasons. In general, contrary to the coordinate-wise, conjunctive representation of the transition function, a controller cannot be represented in such a compressed form. Furthermore, the fixed point algorithm determines all allowed state-input combinations. Since the number of possible control actions for every state grows exponentially with respect to the number of vertices, the controller synthesis problem becomes intractable.

VII Conclusion

We have proposed a modification to the abstraction algorithm, which dramatically reduces the computation time for a large class of systems. Our traffic example shows that the time and memory requirements for encoding the transition function grows roughly linearly with respect to the network length. However, synthesizing controllers for high-dimensional systems still suffers from an exponential runtime with respect to the number of input dimensions.

Future research will focus on speeding up the synthesis procedure in order to leverage abstraction-based controller synthesis methods. Additionally, we want to generalize our proposed modifications to continuous-time systems.

Appendix

1.

In the following, we show that Algorithm 1 and Algorithm 2 return the same abstraction. Without loss of generality, let $n$ be the dimension of the state space $\mathcal{X}$. Algorithm 1 returns a Boolean function that is of the form

[TABLE]

According to 2, the over-approximation function ${\mathsf{O}^{\bm{f}}}$ can be decomposed coordinate-wise. As a result, the next state term $\bigvee_{{[x^{+}]}\in{\mathsf{O}^{\bm{f}}}({[x]},{[u]})}\mathsf{V}([\mathcal{X}^{+}])==\mathsf{I}({[x^{+}]})$ can be split into $n$ components, i.e. 17 can be written equivalently as

[TABLE]

The Boolean expression 18 is equivalent to

[TABLE]

Equivalence of 18 and 19 can easily be seen when 19’s conjunction over $i$ is viewed as multiplication and the disjunction over ${[x]}\in[\mathcal{X}]$ and ${[u]}\in[\mathcal{U}]$ is viewed as addition. “Cross-terms” which correspond to different pairs $({[x]},{[u]})\neq({[x]}^{\prime},{[u]}^{\prime})$ must equate to Boolean false.

Observe that the output of Algorithm 2 and Boolean function 19 are both conjunctions over coordinates $i\in\{1,\ldots,n\}$. If for each coordinate the clauses are equivalent then Algorithm 1 and Algorithm 2’s outputs are equivalent, i.e. both algorithms return the same abstraction.

Consider the $i^{\text{th}}$ component of the system of update equations $\bm{f}$ which depends only on a subset of variables. The associated box partition projections of the update equation $f_{i}$ onto the lower dimensional domains are given by $\mathsf{P}^{\bm{f}}_{i}([\mathcal{X}])$ and $\mathsf{P}^{\bm{f}}_{i}([\mathcal{U}])$, respectively. Any state and input variables that do not lie in these two domains do not change the over-approximation of $f_{i}$. Therefore, any of these non-dependent variables with corresponding indices can removed from 19 without changing the Boolean expression. This removal results in

[TABLE]

which is equivalent to Algorithm 2’s output. Therefore, Algorithm 1 and Algorithm 2 return the same abstraction.

∎