Exact and Efficient HamiltonJacobibased Guaranteed Safety Analysis via System Decomposition
Abstract
HamiltonJacobi (HJ) reachability is a method that provides rigorous analyses of the safety properties of dynamical systems. This method has been successfully applied to many lowdimensional dynamical system models such as coarse models of aircraft and quadrotors in order to provide safety guarantees in potentially dangerous scenarios. These guarantees can be provided by the computation of a backward reachable set (BRS), which represents the set of states from which the system may be driven into violating safety properties despite the system’s best effort to remain safe. Unfortunately, HJ reachability is not practical for highdimensional systems because the complexity of the BRS computation scales exponentially with the number of state dimensions. Although numerous approximation techniques are able to tractably provide conservative estimates of the BRS, they often require restrictive assumptions about system dynamics without providing an exact solution. In this paper we propose a general method for decomposing dynamical systems. Even when the resulting subsystems are coupled, relatively highdimensional BRSs that were previously intractable or expensive to compute can now be quickly and exactly computed in lowerdimensional subspaces. As a result, the curse of dimensionality is alleviated to a large degree without sacrificing optimality. We demonstrate our theoretical results through two numerical examples: a 3D Dubins Car model and a 6D Acrobatic Quadrotor model.
I Introduction
As the presence of safetycritical systems in everyday life has grown, so has the importance for the verification of these systems. Within the next decade we expect to see a rapid increase in the use of safetycritical systems such as autonomous cars, unmanned aerial vehicles, and other robots. Given the number and density of autonomous systems expected in civilian space, higherfidelity models are needed to more accurately characterize these systems so that safety can be guaranteed. In addition, analysis of higherdimensional dynamical system models has the potential to provide valuable insight into the behavior of system states that are frequently ignored to keep the system dimensionality low. Thus, tractable verification tools that are not overly conservative are urgently needed.
Optimal control and differential game theory are powerful tools for the verification of nonlinear systems due to their flexibility with respect to system dynamics, treatment of unknown disturbances, and guaranteed optimality [1, 2, 3, 4]. Reachability analysis is core to these methods; here, the goal is to compute the backward reachable set (BRS), defined as the set of states from which the system can be driven into some unsafe set despite using the optimal control to avoid the unsafe set. HamiltonJacobi (HJ) reachability has been successfully used to guarantee safety for lowdimensional systems in application such as pairwise collision avoidance [2], automated aerial refueling [5], and many others [6, 7]. HJ reachability theory is also very convenient to use due to the many numerical tools available to obtain optimal solutions [8, 9, 10].
Despite these advantages, HJ reachability can be impractical for many highdimensional systems due to issues with scaling. HJ reachabilitybased methods involve solving a partial differential equation (PDE) or variational inequality on a grid representing a numerical discretization of the state space. As a result, the computation complexity scales exponentially with the system dimension. Application of current formulations of HJ reachability is limited to systems with approximately five dimensions or fewer, making the verification of most highdimensional system models intractable.
For the analysis of highdimensional systems, a number of approximation techniques exist. Unfortunately, these techniques usually place strong assumptions on system dynamics, such as requiring a polynomial form [11, 12], a linear form [13, 14], or a Hamiltonian that is only dependent on the control variable [15]. Other methods that are less restrictive in terms of system dynamics include [16], which works with projections, and [17], which involves treating system states as disturbances. In all of the methods mentioned so far, varying degrees of approximation or conservatism is introduced. Under some special scenarios such as those outlined in [18] or [19], a small dimensionality reduction may be possible when obtaining exact optimal solutions.
The previous methods either are forced to trade off between optimality and computation complexity or provide only a small dimensionality reduction. In contrast, this paper presents the selfcontained subsystem (SCS) formulation for computing exact, optimal solutions of systems with dynamics while drastically reducing dimensionality. Motivated by the need to provide safety guarantees, we compute BRSs in lowerdimensional subspaces of the full system state space, and then combine these lowdimensional BRSs to exactly construct the fulldimensional BRS. The fulldimensional BRS can be exactly constructed through back projections of the lowerdimensional BRSs even with coupling between the different subsystems. Furthermore, the theory we present in this paper is compatible with any other method such as [17] and [18]. When different methods are combined together, even more substantial dimensionality reduction can be achieved.
This paper will be presented as follows:

Next, in Section IV we present the SCS formulation, our main theoretical result. We describe how BRSs in lowerdimensional subspaces can be combined to construct the fulldimensional BRS exactly.

Finally, in Section V we present two numerical examples: a lowdimensional 3D Dubins Car example to validate our theory and a highdimensional 6D Acrobatic Quadrotor example that was previously intractable using standard methods.
Ii Background
There are several HJ formulations that can compute BRSs exactly when the system dimensionality is low. Although these methods have been successfully used for lowerdimensional systems, they become intractable when the system dimension is greater than approximately five. In this section, we give a brief overview to provide a starting point on which we build the new proposed theory.
Iia Full System Dynamics
Definition 1
Full system. Let be the state variable of the system under consideration. We call this system the “full system,” or just “system” for short. The evolution of the state of the full system satisfies the ordinary differential equation (ODE)
(1)  
For clarity, we assume that the state space is , but our theory also applies to systems with periodic state dimensions such as angles. The control is denoted by , with the control function being drawn from the set of measurable functions^{1}^{1}1A function between two measurable spaces and is said to be measurable if the preimage of a measurable set in is a measurable set in , that is: , with algebras on ,.:
(2) 
IiB Backward Reachable Set
In this paper, we consider a common definition of the BRS relevant for guaranteeing safety. Intuitively, the BRS represents the set of states from which the system can be driven into an unsafe set at a particular time. For our definition of BRS, we stipulate that the system be driven to for all control functions . In this case, the unsafe set can often be interpreted as a set of states to be avoided (such as an obstacle), and the BRS represents the set of states that leads to the system entering the unsafe set despite all possible control functions. We now formally define the BRS.
Definition 2
Backward reachable set. We denote the BRS , and define it as follows:
(4) 
IiC The Full Formulation for Computing the BRS
There are various similar HJ formulations such as [1, 2, 4], and [22] that cast the reachability problem as an optimal control problem and directly compute the BRS in the full state space of the system. These numerical solutions to the optimal control problem involve solving an HJ PDE on a grid that represents a discretization of the state space. Although these methods are not scalable beyond relatively lowdimensional systems, they form the foundation on which we will build our theory. We now briefly summarize the necessary details related to the HJ PDEs, and what their solutions represent in terms of the cost function and value function of the corresponding optimal control problem.
Let the unsafe set be represented by the implicit surface function such that the unsafe set is the zero sublevel set of the implicit surface function: . Such a function always exists since we can choose to be the signed distance function from . Examples of implicit surface functions are shown as colored surfaces in Fig. 1, with the boundary of the corresponding sets they represent shown in black.
Consider the optimal control problem given by
(5)  
with the optimal control being given by
(6) 
It is wellknown that the value function is the implicit surface function representing : .
The Hamiltonian in (7) is given by
(8) 
Fig. 1 shows an illustration of HJ reachability. , the implicit surface function representing , and the value function , the implicit surface function representing the BRS , are shown as the blue and light gray surfaces respectively. The unsafe set and the BRS are the zero sublevel sets of these two surface functions; the boundaries of and are shown in black. Once the value function is computed, the optimal control (6) can be obtained by the expression
(9) 
We state the following algorithm for clarity and convenience:
Iii Problem Formulation
In this paper, we seek to obtain the BRS in Definition 2 via computations in a lowerdimensional subspace under the assumption that the system (1) can be decomposed into SCSs. Such a decomposition can be commonly found, since many systems involve components that are loosely coupled. In particular, in the dynamics of many vehicles, the evolution of the position variables is often weakly coupled though other variables such as heading.
We now proceed with some essential definitions required to precisely state our main results.
Iiia Definitions
IiiA1 Subsystem Dynamics
Let the system be partitioned as follows:
(10)  
Note that could be zero, and . We call the variables the “state partitions”, or just “partitions”, of the system.
Define the SCS states as follows:
(11)  
It is important to note that and in general have overlapping states in the partition . Note that our theory is applicable to any finite number of subsystems defined in the analogous way; however, for clarity and without loss of generality, in this paper we will assume that there are two subsystems.
For convenience, we have assumed that , but as previously mentioned, our theory also applies to systems with periodic state dimensions.
Definition 3
Selfcontained subsystem. We call each of the systems with states evolving according to (12) a “selfcontained subsystem” (SCS), or just ”subsystem” for short.
(12)  
Intuitively (12) means that the evolution of states in each subsystem depend only on the states in that subsystem: for example, the evolution of depends only on the states in . However, the two subsystems are coupled through the state partition . Note that the subsystem controls and depend on how the control inputs appear in subsystem states and , and may not exist in some subsystems. For example, consider the dynamics of a Dubins Car:
(13) 
with state and control . The state partitions are . The subsystems and the subsystem controls are
(14)  
where the overlapping state is .
The subsystem control signal spaces and control function spaces are defined appropriately according to the full system control signal and function spaces and based on how the control enters the dynamics of the subsystems. For another example of a system decomposed into two selfcontained subsystems, see (34) and (35).
Although there may be common or overlapping states in and , the evolution of each subsystem does not depend on the other explicitly. In fact, if we for example entirely ignore the subsystem , the evolution of the subsystem is welldefined and can be considered a full system on its own; hence, each subsystem is selfcontained.
IiiA2 Projection Operators
Define the projection of a state onto a subsystem state space as
(15) 
For convenience, we will define the projection operator applied on sets :
(16) 
Since we will aim to relate the BRSs of the subsystems to the BRS of the full system, we also define the back projection operator as
(17) 
We will also apply the back projection operator on sets. In this case, we abuse notation and define the back projection operator on some set as
(18) 
IiiA3 Subsystem Trajectories
Since each subsystem in (12) is selfcontained, we can denote the subsystem trajectories . The subsystem trajectories satisfy the subsystem dynamics and initial condition:
(19)  
The full system trajectory and subsystem trajectories are simply related to each other via the projection operator:
(20) 
where .
IiiB Goals of This Paper
We assume that the full system unsafe set can be written in terms of the subsystem unsafe sets in the way depicted in Fig. 3:
(21) 
where the full unsafe set is the intersection of the back projections of subsystem unsafe sets. In practice, this is not a strong assumption since many obstacles can be accurately modeled as rectangular prisms in position space, or hyperrectangles in the full state space. In fact, the unsafe set described by (21) turns out to only be rectangular in the nonoverlapping states, and can be arbitrarily shaped in the overlapping states. In addition, such an assumption is reasonable since the fulldimensional unsafe set should at least be representable in some way in the lowerdimensional spaces. However, in the worst case, taking always leads to a conservative approximation of the constructed BRS that overapproximates the true BRS. Also note that with the definition in (21), we have that .
Next, we define the subsystem BRSs the same way as in (4), but with the subsystems in (12) and subsystem unsafe sets , respectively:
(22) 
Given a system in the form of (12) with unsafe set that can be represented by (21), our goal is to compute the fulldimensional BRS by performing computations in the lowerdimensional subspaces. Specifically, we would like to first compute the subsystem BRSs , and then construct the full system BRS exactly. This process dramatically reduces computation complexity by decomposing the higherdimensional system into two lowerdimensional subsystems. Specifically, we will show that if the unsafe set can be decomposed in the way described by (21), then the fulldimensional BRS is decomposable in the same way:
(23) 
It is important to note that if the subsystem states have no overlapping states (and are therefore decoupled), the above statement is relatively intuitive and easy to show; however, when the subsystems have the overlapping states in the partition , they are coupled to each other through these overlapping states. Our main result in this paper proves that despite this coupling, (23) still holds.
Iv SelfContained Subsystems
With the background and definitions established, we now show the main result in a theorem, which relates lowerdimensional BRSs to the fulldimensional BRS we would like to compute. The consequence of the theorem is that for systems of the form (12), one can obtain the exact fulldimensional BRS by first computing the lowerdimensional BRSs , and then constructing the fulldimensional BRS via (23). We first prove a lemma involving a key property of the projection operator.
Lemma 1
Let for some subsystem . Then,
(24) 
Forward direction: Suppose , then trivially (the that “exists” is just itself). By the definition of back projection in (18), we have .
Backward direction: Suppose , then by the definition of back projection in (18), we have .
Let such an be denoted , and suppose . Then, we must have , which is a contradiction, since .
Corollary 1
If , then
(25) 
Theorem 1
System decomposition for computing the BRS. Suppose that the full system in (1) can be decomposed into the form of (12), then
(26)  
Proof 1
We will prove Theorem 1 by proving the following equivalent statement:
(27) 
By the definition of BRS in (4), we have
(28) 
Consider the property (20), and let
(29)  
Noting that and using Corollary 1, we have the following equivalent statement in terms of the subsystem trajectories:
(30) 
which, by the definition of the subsystem BRS (22), is in turn equivalent to
(31) 
By Lemma 1, this is equivalent to
(32) 
With the above theorem, we now summarize our main theoretical result and its consequences with the following algorithm:
Algorithm 2
SCS formulation. Given an unsafe set that can be decomposed as and SCSs with dynamics in the form (12), the HJbased SCS formulation for computing the BRS is given in the following algorithm:
V Numerical Examples
We now present two numerical examples to illustrate our method. For each example, we present a common dynamical system that can be decomposed into the form of (12). The first example, the 3D Dubins Car, illustrates that our decomposition method produces the exact fulldimensional BRS at a substantially lower computation cost. The second example, the 6D Acrobatic Quadrotor, demonstrates that our technique enables the exact computation of a BRS that was previously intractable to compute with the full formulation.
Va Dubins Car
The Dubins Car is a wellknown system whose dynamics are given by (13). This system is only 3D, and its BRS can be tractably computed in the fulldimensional space, so we use it to compare the full formulation with the SCS formulation.
As previously mentioned, the Dubins Car dynamics can be decomposed according to (14).
For this example, we computed the BRS from the unsafe set representing the set of positions near the origin in both the and dimensions. More concretely, our unsafe set was defined to be . Such an unsafe set can be used to model an obstacle that the vehicle must avoid. Given the unsafe set, the interpretation of the BRS is the set of states from which a collision with the obstacle may occur after a duration of .
From , we computed the BRS of time horizon . The resulting full formulation BRS is shown in Fig. 4 as the red surface which appears in the two subplots on the right.
To compute the BRS using the SCS formulation, note that the unsafe set can be written as , with
(33)  
From these lowerdimensional unsafe sets, we computed the lowerdimensional BRSs and , and then constructed the fulldimensional BRS using Theorem 1: . The subsystem BRSs and their back projections are shown in magenta and green in the left subplot of Fig. 4. The constructed BRS is shown in the three left subplots of Fig. 4 as the black mesh.
In the middleright plot of Fig. 4, we superimpose the fulldimensional BRS computed using the two methods. We show the comparison of the computation results viewed from several different angles in Fig. 5. The results are indistinguishable. However, when using the SCS formulation, Theorem 1 allows the computation to be done significantly faster in lowerdimensional subspaces. An additional benefit of the SCS formulation is that in the numerical methods for solving the HJ PDE, the amount of numerical dissipation increases with the number of state dimensions. Thus, computations in lowerdimensional subspaces lead to a slightly more accurate numerical solution.
The computation benefits of using the SCS formulation can be seen from Fig. 6. Both subplots show the computation time in seconds versus the number of grid points per dimension in the numerical computation. From the top subplot, one can easily see that the direct computation of the full formulation BRS in 3D becomes very timeconsuming as the number of grid points per dimension is increased, while the computation using the SCS formulation hardly takes any time in comparison. The bottom subplot shows the same data, but on a loglog scale for more detail. Directly computing the BRS with 251 grid points per dimension using the full formulation took approximately 80 minutes, while computing the BRS using the SCS formulation is faster by several orders of magnitude: the computation only took approximately 30 seconds! The computations were timed on a desktop computer with an Intel Core i72600K processor and 16GB of randomaccess memory.
VB The 6D Acrobatic Quadrotor
This example illustrates the ability of the SCS formulation to produce BRSs for highdimensional systems that would be otherwise intractable to analyze by current HJbased methods. In [25], a 6D Acrobatic Quadrotor model used to perform backflips was simplified into a series of smaller hybrid models due to the intractability of computing a BRS over a 6D state space. Using the new SCS formulation we can accurately compute a BRS for the full 6D system.
The 6D Acrobatic Quadrotor’s state is ; the dynamics are given by [25]:
(34) 
where , , and represent the quadrotor’s horizontal, vertical, and rotational positions, respectively. Their derivatives represent the velocity with respect to each corresponding positional state. The inputs and represent the thrust exerted on either end of the quadrotor, and the constant system parameters are for mass, for translational drag, for rotational drag, for acceleration due to gravity, for the length from the quadrotor’s center to an edge, and for moment of inertia.
The state partitions of this system are . Using the SCS formulation, we decompose the full system into the following set of subsystems:
(35)  
For this example we will compute the BRS that describes the set of initial conditions from which the system may enter the unsafe set after a given time period despite best possible control. We define the unsafe set as a square of length 1m centered at described by . This can be interpreted as a positional box centered at the origin that must be avoided for all angles and velocities. From the unsafe set, we define such that . This unsafe set must be decomposed to provide a suitable unsafe set for each subsystem. This is done by letting be
(36)  
The BRS of each 4D subsystem is computed and then combined into the 6D BRS using the SCS formulation. To visually depict the 6D BRS, 3D slices of the BRS along the positional and velocity axes were computed. Fig. 7 shows a 3D slice in space at m/s, rad/s. The dark blue set represents the unsafe set , with the BRS in light blue.
In Fig. 8, 3D slices in space are visualized at m, rad. These colored sets represent the BRS at different points in time.
Vi Conclusions and Future Work
The SCS formulation that we proposed for computing BRSs significantly reduces computation burden, and makes many previously intractable computations possible. At the same time, the computation savings do not come at the cost of optimality: the fulldimensional BRS can be computed exactly in lowerdimensional subspaces. The construction of the fulldimensional BRS from lowerdimensional BRSs is exact even when the subsystem dynamics are coupled.
The SCS formulation will be the basis for future system decomposition methods, for which we already have several other preliminary theoretical results. These results include other definitions of BRSs, such as those used for reaching, instead of avoiding, a set; incorporation of disturbances into the problem formulation; and treatment of reachable tubes. In the future, we plan to apply the theory to a larger number of practical systems in these different settings.
References
 [1] E. N. Barron, “Differential Games with Maximum Cost,” Nonlinear analysis: Theory, methods & applications, pp. 971–989, 1990.
 [2] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin, “A timedependent HamiltonJacobi formulation of reachable sets for continuous dynamic games,” IEEE Transactions on Automatic Control, vol. 50, no. 7, pp. 947–957, July 2005.
 [3] K. Margellos and J. Lygeros, “HamiltonJacobi Formulation for ReachAvoid Differential Games,” IEEE Transactions on Automatic Control, vol. 56, no. 8, Aug 2011.
 [4] O. Bokanowski and H. Zidani, “Minimal time problems with moving targets and obstacles,” {IFAC} Proceedings Volumes, vol. 44, no. 1, pp. 2589 – 2593, 2011.
 [5] J. Ding, J. Sprinkle, S. S. Sastry, and C. J. Tomlin, “Reachability calculations for automated aerial refueling,” in IEEE Conference on Decision and Control, Cancun, Mexico, 2008.
 [6] E. M. Vaisbord and V. I. Zhukovskii, Introduction to Multiplayer Differential Games and Their Applications. Routledge, 1988.
 [7] I. Mitchell, “Application of level set methods to control and reachability problems in continuous and hybrid systems,” Ph.D. dissertation, Stanford University, 2002.
 [8] J. A. Sethian, “A fast marching level set method for monotonically advancing fronts,” Proceedings of the National Academy of Sciences, vol. 93, no. 4, pp. 1591–1595, 1996.
 [9] S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces. SpringerVerlag, 2002.
 [10] I. Mitchell, A Toolbox of Level Set Methods, 2009. [Online]. Available: {}{}}{http://people.cs.ubc.ca/~mitchell/ToolboxLS/index.html}{cmtt}
 [11] A.~Majumdar and R.~Tedrake, Algorithmic Foundations of Robotics X: Proceedings of the Tenth Workshop on the Algorithmic Foundations of Robotics. Springer Berlin Heidelberg, 2013, ch. Robust Online Motion Planning with Regions of Finite Time Invariance, pp. 543558.
 [12] T.~Dreossi, T.~Dang, and C.~Piazza, ``Parallelotope bundles for polynomial reachability,'' in 19th International Conference on Hybrid Systems: Computation and Controls, 2016.
 [13] A.~B. Kurzhanski and P.~Varaiya, ``Ellipsoidal techniques for reachability analysis: internal approximation,'' Systems & control letters, vol.~41, no.~3, pp. 201211, 2000.
 [14] J.~N. Maidens, S.~Kaynama, I.~M. Mitchell, M.~M. Oishi, and G.~A. Dumont, ``Lagrangian methods for approximating the viability kernel in highdimensional systems,'' Automatica, July 2013.
 [15] J.~Darbon and S.~Osher, ``Algorithms for overcoming the curse of dimensionality for certain hamiltonjacobi equations arising in control theory and elsewhere,'' arXiv preprint arXiv:1605.01799, 2016.
 [16] I.~M. Mitchell and C.~J. Tomlin, ``Overapproximating reachable sets by hamiltonjacobi projections,'' Journal of Scientific Computing, vol.~19, no. 13, pp. 323346, 2003.
 [17] M.~Chen*, S.~Herbert*, and C.~J. Tomlin, ``Fast reachable set approximations via state decoupling disturbances,'' in 55th IEEE Conference on Decision and Control, 2016.
 [18] I.~M. Mitchell, ``Scalable calculation of reach sets and tubes for nonlinear systems with terminal integrators: A mixed implicit explicit formulation,'' in Proceedings of the 14th International Conference on Hybrid Systems: Computation and Control, 2011, pp. 103112.
 [19] J.~F. Fisac, M.~Chen, C.~J. Tomlin, and S.~S. Shankar, ``Reachavoid problems with timevarying dynamics, targets and constraints,'' in 18th International Conference on Hybrid Systems: Computation and Controls, 2015.
 [20] L.~C. Evans and P.~E. Souganidis, ``Differential games and representation formulas for solutions of HamiltonJacobiIsaacs equations,'' Indiana Univ. Math. J., vol.~33, no.~5, pp. 773797, 1984.
 [21] E.~A. Coddington and N.~Levinson, Theory of ordinary differential equations. Tata McGrawHill Education, 1955.
 [22] P.~Varaiya, ``On the existence of solutions to a differential game,'' SIAM Journal on Control, vol.~5, no.~1, pp. 153162, 1967.
 [23] M.~G. Crandall and P.L. Lions, ``Viscosity solutions of HamiltonJacobi equations,'' Transactions of the American Mathematical Society, vol. 277, no.~1, pp. 142, 1983.
 [24] M.~G. Crandall, L.~C. Evans, and P.~L. Lions, ``Some properties of viscosity solutions of hamiltonjacobi equations,'' Transactions of the American Mathematical Society, vol. 282, no.~2, p. 487, April 1984.
 [25] J.~H. Gillula, G.~M. Hoffmann, H.~Huang, M.~P. Vitus, and C.~J. Tomlin, ``Applications of hybrid reachability analysis to robotic aerial vehicles,'' International Journal of Robotics Research, vol.~30, no.~3, pp. 335354, March 2011.