Set-boundary based Reachability Analysis Toolbox in Python
Reachability analysis, which involves computing reachable state sets, plays a fundamental role in the temporal verification of nonlinear systems. Overly pessimistic over-approximations, however, render many temporal properties unverifiable in practice. This pessimism mainly arises due to the wrapping effect, which is the propagation and accumulation of over-approximation error through the iterative computation in the construction of reachable sets. As the extent of the wrapping effect correlates strongly with the volume of the initial set, techniques that partition the initial state space and independently compute reachable sets of those partitions are often used to reduce the wrapping effect, especially for large initial sets or/and large time horizons. Such partitioning may, however, induce extensive demand on computation time and memory, often rendering the existing reachability analysis techniques not suitable for complex real-world applications. Not being forced to explore the full, i.g. exponential in the dimensionality, number of partitions could help such procedures tremendously. This is the theme of this tool, which implements the so-called 'set-boundary based method' that explores means of computing the full reachable state space based on state-exploratory analysis of just a small sub-volume of the initial state set, namely a set enclosing its boundary. For theoretical analysis, please refer to 'Bai Xue, Arvind Easwaran, Nam-Joon Cho and Martin Fränzle.Reach-Avoid Verification for Nonlinear Systems Based on Boundary Analysis. IEEE Transactions on Automatic Control (IEEE TAC), vol. 62: 3518--3523, 2017.' and 'Bai Xue, Qiuye Wang, Shenghua Feng, and Naijun Zhan. Over-and underapproximating reach sets for perturbed delay differential equations. IEEE Transactions on Automatic Control (IEEE TAC), vol.66: 283--290,2020.'
The set-boundary based method can be used to perform reachability analysis for systems modelled by
- ordinary differential equations (subject to Lispchitz continuous perturbations)
- delay differential equations (subject to Lispchitz continuous perturbations)
- neural ordinary differential equations
PyBDR is now available as an installable Python package. We strongly recommend using Conda for installation, as it automatically handles all system dependencies.
# Create conda environment (i.e. pybdr here) with all dependencies
conda env create -f environment.yml --name pybdr
# Activate the environment
conda activate pybdr
# Install PyBDR
pip install -e .That's it! The environment.yml file automatically installs the required cddlib C library and all Python dependencies.
If you already have a conda environment:
# Activate your environment
conda activate your_env_name
# Install cddlib system dependency
conda install -c conda-forge cddlib
# Install PyBDR
pip install -e .For detailed installation instructions and troubleshooting, please see INSTALL.md.
PyBDR depends on pypoman, which requires the cddlib C library. Conda makes this dependency installation seamless. While pip-only installation is possible, it requires manual installation of system libraries and is not recommended.
We recommend using PyCharm or VS Code as the IDE for development. To ensure a smoother installation and running of third-party libraries, we advise users to use miniconda and create a virtual environment. The steps for this process are as follows:
First, open the user's current working directory, and use the command
conda create -n pybdr_labto initialize a virtual test environment called "pybdr_lab".
After the virtual environment has been created, the user needs to activate it before running any third-party libraries. This can be done using the command
conda activate pybdr_labBy activating the virtual environment, the user ensures that any package installations and other commands will run within the virtual environment, rather than the system environment.
Now, the user can install the necessary third party libraries in this virtual environment using the following commands.
conda install matplotlib
conda install -c conda-forge numpy
conda install -c conda-forge cvxpy
conda install scipy
conda install -c mosek mosek
pip install pypoman
conda install -c anaconda sympy
pip install open3dFor the reason we may use Mosek as a solver for optimisation, we highly recommend you to apply for an official personal licence, the steps for which can be found at this link.
The tool provides sample files which serve as demonstrations of the proper utilization for computing reachable sets. These sample files serve as a reference point for users to grasp the process of modifying the dynamics and parameters necessary for reachability analysis. This feature aids users in experimenting with their analyses, allowing them to assess the impact of different settings on the overall computation of the reachable sets
For example, consider the following dynamic system:
import numpy as np
from pybdr.algorithm import ASB2008CDC, ASB2008CDC
from pybdr.util.functional import performance_counter, performance_counter_start
from pybdr.geometry import Zonotope, Interval, Geometry
from pybdr.geometry.operation import boundary, cvt2
from pybdr.model import *
from pybdr.util.visualization import plot
# settings for the computation
options = ASB2008CDC.Options()
options.t_end = 6.74
options.step = 0.005
options.tensor_order = 3
options.taylor_terms = 4
options.u = Zonotope.zero(1, 1)
options.u_trans = np.zeros(1)
# settings for the using geometry
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
Zonotope.ORDER = 50
z = Interval([1.23, 2.34], [1.57, 2.46])
x0 = cvt2(z, Geometry.TYPE.ZONOTOPE)
xs = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
this_time = performance_counter_start()
ri_without_bound, rp_without_bound = ASB2008CDC.reach(vanderpol, [2, 1], options, x0)
this_time = performance_counter(this_time, 'reach_without_bound')
ri_with_bound, rp_with_bound = ASB2008CDC.reach_parallel(vanderpol, [2, 1], options, xs)
this_time = performance_counter(this_time, 'reach_with_bound')
# visualize the results
plot(ri_without_bound, [0, 1])
plot(ri_with_bound, [0, 1])| With Boundary Analysis (BA) | No Boundary Analysis (NBA) |
|---|---|
 |
 |
For large initial sets,
| System | Code | Reachable Sets (Orange-NBA,Blue-BA) |
|---|---|---|
| synchronous machine | benchmark_synchronous_machine_cmp.py |  |
| Lotka Volterra model of 2 variables | benchmark_lotka_volterra_2d_cmp.py |  |
| Jet engine | benchmark_jet_engine_cmp.py |  |
For large time horizons, i.e. consider the system Brusselator
For more details about the following example, please refer to our code.
| Time instance | With Boundary Analysis | Without Boundary Analysi |
|---|---|---|
| t=5.4 |  |
 |
| t=5.7 |  |
 |
| t=6.0 |  |
 |
| t=6.1 |  |
Set Explosion Occurred! |
For example, consider a neural ODE with the following parameters and
import numpy as np
from pybdr.algorithm import ASB2008CDC
from pybdr.geometry import Zonotope, Interval, Geometry
from pybdr.model import *
from pybdr.util.visualization import plot, plot_cmp
from pybdr.geometry.operation import boundary, cvt2
from pybdr.util.functional import performance_counter_start, performance_counter
# settings for the computation
options = ASB2008CDC.Options()
options.t_end = 1
options.step = 0.01
options.tensor_order = 2
options.taylor_terms = 2
options.u = Zonotope([0], np.diag([0]))
options.u_trans = options.u.c
# settings for the using geometry
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
Zonotope.ORDER = 50
z = Interval([0, -0.5], [1, 0.5])
x0 = cvt2(z, Geometry.TYPE.ZONOTOPE)
xs = boundary(z, 2, Geometry.TYPE.ZONOTOPE)
print(len(xs))
this_time = performance_counter_start()
ri_without_bound, rp_without_bound = ASB2008CDC.reach(neural_ode_spiral1, [2, 1], options, x0)
this_time = performance_counter(this_time, "reach_without_bound")
ri_with_bound, rp_with_bound = ASB2008CDC.reach_parallel(neural_ode_spiral1, [2, 1], options, xs)
this_time = performance_counter(this_time, "reach_with_bound")
# visualize the results
plot_cmp([ri_without_bound, ri_with_bound], [0, 1], cs=["#FF5722", "#303F9F"])In the following table, we show the reachable computed with boundary analysis and without boundary analysis on different time instance cases.
| Time Instance | With Boundary Analysis | Without Boundary Analysis |
|---|---|---|
| t=0.5 |  |
 |
| t=1.0 |  |
 |
| t=1.5 |  |
Set Explosion Occured! |
Two modes of computation are supported by the tool for reachable sets. One mode is to compute the reachable set of evolved states using the entire initial set in a set propagation manner, while the other mode is to compute the reachable set of evolved states based on the boundary of the initial state set.
The computation may be slow for several reasons such as large computational time intervals, small steps, high Taylor expansion orders, or a large number of state variables.
To accelerate the computations, experiments can be performed with a smaller computational time horizon, a smaller order of expansion (such as 2), and a larger time step. Then gradually increase the computational time horizon and order of expansion based on the results of this setting to achieve the desired set of reachable states at an acceptable time consumption.
To enhance the precision of the reachable set computation, one can split the boundaries of initial sets or increase the order of the Taylor expansion while reducing the step size.
Feel free to contact dingjianqiang0x@gmail.com if you find any issues or bugs in this code, or you struggle to run it in any way.
This project is licensed under the GNU GPLv3 License - see the LICENSE file for details.
@inproceedings{ding2024pybdr,
title={PyBDR: Set-Boundary Based Reachability Analysis Toolkit in Python},
author={Ding, Jianqiang and Wu, Taoran and Liang, Zhen and Xue, Bai},
booktitle={International Symposium on Formal Methods},
pages={140--157},
year={2024},
organization={Springer}
}
When developing this tool, we drew upon models used in other tools for calculating reachable sets, including Flow*, CORA, and various others.