Reachability Analysis Toolbox Based on Boundary Propagation in Python

Online Documents »

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 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.

To incorporate the source code, there is no requirement for installation; simply copy it into the appropriate directory, similar to using a third-party library in Matlab. We strongly suggest utilizing `Pycharm' as an IDE for both development and testing purposes.

We recommend that users use miniconda to initialise a virtual environment for the subsequent installation and running of third party libraries. The steps are as follows.

In the user's current working directory, the user can initialise a virtual test environment called pybdr_lab using the following command.

conda create -n pybdr_lab

After the virtual environment has been initialized, the user activates the virtual test environment named pybdr_lab that has just been initialized with the following command.

conda activate pybdr_lab

Now, you can install the necessary third party libraries in this virtual environment using the following series of 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 open3d

For 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.

Example files are provided to show how to use the tool to calculate reachable sets. Users can refer to the example files provided and modify the dynamics and parameters required for the calculation to see the effect of using different settings for calculating system reachable sets.

For example, consider following dynamic system:

$$ \begin{align*} \dot{x} &= y \\ \dot{y} &= (1-x^2)y-x \end{align*} $$

import numpy as np

from pyrat.algorithm import ASB2008CDC
from pyrat.dynamic_system import NonLinSys
from pyrat.geometry import Zonotope, Interval, Geometry
from pyrat.model import *
from pyrat.util.visualization import plot

# init dynamic system
system = NonLinSys(Model(vanderpol, [2, 1]))

# 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.r0 = [Zonotope([1.4, 2.4], np.diag([0.17, 0.06]))]
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

# reachable sets computation
ti, tp, _, _ = ASB2008CDC.reach(system, options)

# visualize the results
plot(tp, [0, 1])

• The tool supports two modes of computation for reachable sets. The first mode calculates the reachable set of evolved states based on the entire set of states in a set propagation manner. The second mode calculates the reachable set of evolved states based on the initial state set boundary. There are several factors contributing to the slow computation of the tool: large computational time intervals, small steps, high Taylor expansion orders, and high dimensionality of the dynamical system. To address these issues, it is recommended to conduct experiments with a smaller computational time horizon, a lower expansion order (e.g., 2), and a larger time step. Based on the results obtained from this initial setting, gradually increase the computational time horizon and expansion order to obtain the desired set of reachable states while managing time consumption effectively.

• Due to the wrapping effect of set propagation-based algorithms, it is inevitable that the computed reachable state sets may be overly conservative under inappropriate settings, thereby preventing the attainment of reachable state sets that satisfy the desired requirements. To address these issues, several strategies can be employed to enhance the accuracy of the reachable set calculation. One technique involves dividing the boundary of the initial set into smaller cells. Furthermore, reducing the step size and increasing the order of the Taylor expansion can also contribute to improved accuracy.

Feel free to contact [email protected] 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.

This tool has been developed with reference to models used in tools such as Flow*, CORA and other reachable set calculation tools.