Exact and Efficient HamiltonJacobi Reachability for Decoupled Systems
Abstract
Reachability analysis is important for studying optimal control problems and differential games, which are powerful theoretical tools for analyzing and modeling many practical problems in robotics, aircraft control, among other application areas. In reachability analysis, one is interested in computing the reachable set, defined as the set of states from which there exists a control, despite the worst disturbance, that can drive the system into a set of target states. The target states can be used to model either unsafe or desirable configurations, depending on the application. Many HamiltonJacobi formulations allow the computation of reachable sets; however, due to the exponential complexity scaling in computation time and space, problems involving approximately 5 dimensions become intractable. A number of methods that compute an approximate solution exist in the literature, but these methods trade off complexity for optimality. In this paper, we eliminate complexityoptimality tradeoffs for timeinvariant decoupled systems using a decoupled HamiltonJacobi formulation that enables the exact reconstruction of high dimensional solutions via low dimensional solutions of the decoupled subsystems. Our formulation is compatible with existing numerical tools, and we show the accuracy, computation benefits, and an application of our novel approach using two numerical examples.
I Introduction
Optimal control problems and differential games have been extensively studied [1, 2, 3, 4], and have received growing interest in the recent past. These powerful theoretical tools allow us to analyze a variety of real world problems, including path planning, collision avoidance, safety verification, among other applications in robotics, aircraft control, security, and other domains [5, 6, 7, 8].
In an optimal control problem, one aims to drive a controlled dynamical system into a set of states called the target set; depending on the application, the target set can model the set of either desirable or undesirable configurations. In a reachability framework, one aims to determine the backwards reachable set, defined as the set of states from which a control exists to drive the system into the target set. Differential games involve two adversarial players (Player 1 and Player 2). Player 2 seeks to drive a system to a target set, while Player 1 seeks to prevent Player 2 from doing so. One again aims to determine the backwards reachable set, which in this case is defined as the set of states from which a control from Player 2 exists to drive the system into the target set, despite the optimal adversarial control from Player 1.
Reachability is an effective way to analyze optimal control problems and differential games because it provides guarantees on system performance and safety. Reachability problems involving one player can be posed as a minimum (maximum) cost game where the player minimizes the minimum value over time of some cost function representing the proximity to the target set. In the case of a differential game, Player 1 maximizes the minimum cost over time, while Player 2 minimizes it. [1] has shown that the backwards reachable set can be obtained by solving a HamiltonJacobi Partial Differential Equation (HJ PDE) with a terminal condition specifying the target set. Many similar formulations of the backwards reachability problems also exist [9, 10, 11]. HJ reachability has been successfully used to solve problems such as aircraft collision avoidance [1], automated inflight refueling [12], and reachavoid games [13, 14].
The techniques for computing backwards reachable sets via solving an HJ PDE are very flexible and can be applied to a large variety of system dynamics when the problem dimensionality is low. Furthermore, many numerical tools have been developed to solve these equations, making the HJ approach practically appealing [15, 16, 17, 18]. For higher dimensional problems, various techniques such as those involving projections [19, 20], approximate dynamic programming [21], and occupation measure theory [22] have been proposed. While these approximation techniques alleviate the computation complexity, they give up optimality and sometimes give overly conservative results.
This paper resolves the complexityoptimality tradeoff for timeinvariant systems with decoupled dynamics. We present a decoupled formulation of HJ reachability for decoupled systems, defined in (3). By considering the decoupled component separately and solving lower dimensional HJ PDEs for each subsystem, we reduce the computation time and space complexity substantially. Our approach also exactly recovers the solution to the original, high dimensional PDE.
Ii Problem Formulation
Consider a differential game between two players described by the timeinvariant system
(1) 
where is the system state, is the control of Player 1, and is the control of Player 2. We assume is uniformly continuous, bounded, and Lipschitz continuous in for fixed , and the control functions are drawn from the set of measurable functions
(2)  
We further assume that the system is a decoupled system.
Definition 1
Decoupled system. A system (1) is a decoupled system if it can be split into components, denoted where , that satisfy the following:
(3)  
where is th component of the full state, is th component of the control of Player 1, and is th component of the control of Player 2. Based on this assumption and the assumptions on , we have that is uniformly continuous, bounded, and Lipschitz continuous in for fixed and are measurable. Note that .
Denote system trajectories, which are solutions to (1), as
(4) 
satisfies initial conditions and the following differential equation almost everywhere
(5) 
In our differential game, the goal of Player 2 is to drive the system into some target set , and the goal of Player 1 is to drive the system away from it. The set is represented as the zero sublevel set of a bounded, Lipschitz continuous function , .
Such a function always exists, since we can choose to be a signed distance function; we call the implicit surface function representing the set . In accordance with our decoupled dynamics, we assume that can be represented as a maximum of bounded, Lipschitz continuous functions , where are implicit surface functions representing so that . Note that with the definition of and , we have that
(6) 
Given the decoupled system (3) and the target set in the form (6) represented by , our goal in this paper is to compute the backwards reachable set, , in the lowdimensional space of each of the decoupled components as opposed to in the full system state space . is defined as
(7)  
Remark 1
One may have noticed that if , one would be able to simply find by solving (8) with as the target set, and then obtain . However, it is crucial to observe that we are interested in the case where , in which a simple union of would not yield the correct reachable set .
Iii Solution
Iiia HJ Reachability: Full Formulation
In [1], the authors showed that the backwards reachable set can be obtained as the zero sublevel set of the viscosity solution [25] of the following terminal value HJ PDE:
(8)  
from which we obtain from the bounded, Lipschitz function that is also continuous in both and [23].
[1] and similar approaches, such as [9, 10, 11], are compatible with wellestablished numerical methods [15, 16, 17, 18]. However, these approaches become intractable quickly as the dimensionality of the problem increases. Numerically, the solution is computed on a grid, and the number of grid points increases exponentially with the number of dimensions.
Decoupled dynamics allow for tractable or faster computation of reachable sets in the individual decoupled components. Some authors [19, 20] have proposed methods for combining or stitching together these reachable set components into the full reachable set. These methods work reasonably well, but introduce conservatism in various ways. In the next subsections, we will provide a method for combining solutions to the lower dimensional HJ PDEs to construct the exact full solution in the original HJ PDE.
IiiB HJ Reachability: Decoupled Formulation
Observe that (8) can be viewed as an equation involving two cases. Depending on which of the arguments in the outermost minimum is active, becomes one of (9) or (10):
(9) 
(10) 
This motivates us to define the following sets which characterize which of the outermost minimum operation is active in (8).
(11)  
Note that is the complement of , , for all time. At a given , in , satisfies (9); in , satisfies (10). We now show an important property of and in the Lemma and Corollary below. These will be used to show that our proposed decoupled formulation allows exact computation of , by computation of , value functions of lower dimensional spaces.
Lemma 1
for some such that .
Suppose , then by (9) we have the following:

. Thus, becomes independent of at .

since is an open set, there exists a neighborhood around that is contained in . Thus is also independent of in a neighborhood of .

By (11), we have .
Let and suppose at . Then, by (11), . Since is independent of , this necessarily means that .
This implies in a neighborhood of such that .
(12) 
However, by (9), . In particular, then, we have for any
(13) 
This means that , such that
(14) 
which is a contradiction since is Lipschitz continuous. Therefore, since , we have that by (11).
Corollary 1
.
Suppose but . Since is the complement of , this implies . By Lemma 1, we must have that since , a contradiction.
We can now state our main theorem.
Theorem 1
(15)  
where are the viscosity solutions to
(16)  
and is the smallest time such that , i.e.
(17) 
Case 1: By Lemma 1, we have . Therefore, satisfies (9) , so . Case 2 of this proof would then imply .
Case 2: Consider a target set represented by the zero sublevel set of the function , where . By (8), we have that . Define functions such that , then at , we have , and
(18)  
where is the indicator function that is when its argument is true and otherwise, and is an matrix in of all zeros except for in the rows corresponding to the component where it is the identity matrix in .
(19)  
Equation (19) states that in the region where is the maximum among , we have , where satisfies (16).
Consider auxiliary functions which satisfy, for all and all ,
(20)  
By Corollary 1, we have that and both satisfy the same PDE with the same terminal conditions, . Therefore, .
IiiC Decoupled Formulation Algorithm
This gives us an efficient way to computed by computing , transforming the original dimensional problem of computing into the dimensional problems of computing . Based on the conclusions we drew, the following algorithm exactly computes , which satisfies (8), with the abovementioned computation benefits:

Initialize .

Compute , by solving (16).

Initialize .

Decrement from 0 to ; for each time step :

Set the auxiliary variable
This step correctly computes to be for all .

Update the value function
This step correctly computes to satisfy for all .

IiiD Computation Time and Space Complexity Comparison
For a state space discretization of grid points in each dimension, the computation time complexity decreases from for the original problem in , to for the subproblems in . This is a computation speed improvement of many orders of magnitude.
Directly solving (8) on a computational domain has a space complexity of , where is the number of time steps of being stored, since we need to store an dimensional grid for each of the time steps. The algorithm presented in IIIC involves computing on computation domains . Each thus has a space complexity of , making the overall space complexity .
From , we can then reconstruct in any domain . Thus, by choosing to be a small subset of , we can always avoid additional space complexity. This allows us to access . In practice, one would choose to be in a small region around a state of interest (eg. the current system state), and access the value function as well as its gradient at ; this allows one to determine of whether is in the reachable set based on the sign of , and compute the optimal controls for both Player 1 and Player 2 respectively based on as follows:
(21)  
Note that we do not need to store for all . In fact, in many situations, storing in the entire is infeasible since the space complexity is exponential with the dimension of . With that caveat, we restate our algorithm for computing from , explicitly noting memory allocation, to show that one only needs to store for , where is a very small subset of :

Initialize for .

Compute in by solving 16.

Initialize in a small computation domain .

Decrement from 0 to ; for each time step , perform the following computations in :

IiiE Numerical Implementation
Our proposed decoupled formulation involves solving (16) for each of the subsystems. As already mentioned, many numerical tools already exist for solving (16); we will use the implementation in [15]. For the examples in this paper, we used the numerical schemes below.
For the numerical Hamiltonian in (16), we used the LaxFriedrich approximation [26]. For the numerical spatial derivatives , we used a fifthorder accurate weighted essentially nonoscillatory scheme [26, 27]. For numerical time derivatives , we used a thirdorder accurate total variation diminishing RungeKutta scheme [27, 28].
Computations were done on a computer with an Intel Core i72600K CPU running at 3.4 GHz, with 16 GB of memory.
Iv 4D Quadrotor Collision Avoidance
Consider a simple quadrotor model consisting of two decoupled doubleintegrators:
(22)  
denote the  and position of the quadrotor, and denote the  and velocity. The control signals are the  and acceleration of the quadrotor, constrained to be between and .
Now, consider two quadrotors in a pursuitevasion game, in which the evader (Player 1) aims to avoid collision, while the pursuer (Player 2) aims to cause a collision. The relative coordinates of the two quadrotors are given by the following state variables:
(23)  
Given the above relative state variables, the relative dynamics of the two quadrotors are given by
(24)  
Note that this system is decoupled, with as the first decoupled component, and as the second decoupled component. In the relative coordinates of the two quadrotors, we define the collision set of size , representing the configurations in which the two quadrotors are considered to have collided, as the following set:
(25) 
with the corresponding implicit surface function where . Since we have a decoupled system, let be the following sets:
(26)  
with corresponding implicit surface functions . Then, we have and .
We will set as the target set in our reachability problem, and compute the backwards reachable set from using three methods:

Solve (8) directly in to obtain , whose zero sublevel set represents .

Compute the analytic boundary of the reachable set.
For comparison purposes, for the first two methods we will compute on the computation domain . However, it is important to recall that for our proposed method described in Section IIIB, we can significantly reduce space complexity by only storing a small part of .
Iva Reachable Set
Since the state space of our system is 4D, we visualize various slices of the reachable set, whose boundary is given by . Figures 1 shows these slices. The reachable set boundary computed using our proposed decoupled method is very close to the reachable set boundary computed by solving the full PDE (8) in and to analytic reachable set boundary. Figures 2 zooms in on the plots for a closer look.
IvB Performance
In order to quantify the computation error, we converted into a signed distance function from the boundary . This operation was first proposed in [29] and can be done by solving the reinitialization PDE formulated in [30]; for this operation, we use the implementation in [15]. We then evaluated approximately 24 million analyticallycomputed reachable set boundary points on the ; the resulting values represent how far each of the analyticallycomputed points are from the numericallycomputed boundary. The values of on analytic boundary points are defined as the computation error.
Figure 3 shows the error as a function of grid spacing. In terms of the maximum error (red curve), the decoupled formulation results in a numericallycomputed reachable set boundary that is accurate within the size of the grid spacing (black line). On average, the error is approximately an order of magnitude smaller than the size of grid spacing (blue curve). Furthermore, we can see numerical convergence to the analytic solution as the grid spacing size decreases.
Figure 4 shows the computation time as a function of the number of grid points in each dimension. Here, we can see that the decoupled formulation is orders of magnitude faster than the full formulation, and can be done with many more grid points in each dimension. Lastly, the slopes of curves in the loglog plot show an time complexity for the full formulation, and only for the decoupled formulation.
For the decoupled formulation, when we reconstruct the full value function in 4D (blue curve), the computation time hardly increases compared to when we do not perform full reconstruction (green curve). However, in general, we recommend that the value function in only a region near a state of interest should be computed. Without full reconstruction of the value function, we are able to obtain results with many more grid points (green curve), improving the accuracy of the numerical computation.
V 6D Quadrotor Collision Avoidance
Consider relative dynamics augmented by the velocity of evader quadrotor, given in Equation (27). These dynamics are needed to impose a velocity limit on the quadrotor.
(27)  
For this system, we consider a collision between the two quadrotors, as defined previously, to be unsafe configurations. Here, we denote this set
(28) 
with corresponding implicit surface function .
We also consider configurations in which Player 1 is exceeding a velocity limit of in the  or  directions to be unsafe. This set of configurations is denoted :
(29) 
We define the target set to be the union of the above configurations; represents all unsafe configurations. We will compute the reachable set from the target set as follows:

Compute , the reachable set from target set .

Compute , the reachable set from target set .

Take the union to obtain .
Va Reachable set
Taking the union , we obtain the 6D reachable set. To visualize , we compute 2D slices of the 6D reachable set at various values. This is done without computing the entire 6D reachable set by setting the computation domain to be in a large portion of the plane, at a small range of values. The 2D slices at shown in Figure 5. Each subplot shows two different pairs for a particular . The red boundary represents the slice with , while the blue boundary represents a slice with the evader velocity almost exceeding the limit of . The blue boundaries contain the red ones, because if the evader is already near the velocity limit, it would have more limited capability to avoid collisions.
Figure 6 shows a simulation of the collision avoidance maneuver resulting from the reachable set. The red and blue quadrotors are initially traveling in opposite directions. Consider the situation in which the red quadrotor insists on staying on its intended path. In this case, the blue quadrotor must perform the optimal avoidance control whenever it reaches the boundary of the reachable set (), shown as the blue dotted boundary. When the blue quadrotor is no longer at the boundary of the reachable set (), it is free to perform any control.
Vi Conclusions and Future Work
We have presented a decoupled formulation of the HJ reachability, enabling us to efficiently compute the viscosity solution of the HJ PDE (8), which gives the reachable set for optimal control problems and differential games. When the system dynamics are decoupled, our decoupled formulation allows the exact reconstruction of the solution to the full dimensional HJ PDE by solving the HJ PDEs corresponding to each of the decoupled components. Our novel approach enables the analysis of otherwise intractable problems. In addition, our formulation achieves this without sacrificing any optimality compared to the full formulation. We demonstrated the benefits of our approach using a 4D and 6D quadrotor system.
Footnotes
 A 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 ,.
References
 I. Mitchell, A. Bayen, and C. 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.
 C. Tomlin, J. Lygeros, and S. Shankar Sastry, “A game theoretic approach to controller design for hybrid systems,” Proceedings of the IEEE, vol. 88, no. 7, pp. 949 –970, Jul 2000.
 T. Basar and G. Olsder, Dynamic Noncooperative Game Theory, 2nd ed. Philadelphia, PA: SIAM, 1999.
 J. Lygeros, C. Tomlin, and S. Sastry, “Controllers for reachability specifications for hybrid systems,” Automatica, vol. 35, no. 3, pp. 349 – 370, 1999.
 Department of the Air Force, “United states air force unmanned aircraft systems flight plan 20092047,” oai.dtic.mil, Jan 2009.
 A. Madrigal, “Autonomous robots invade retail warehouses,” http://www.wired.com/wiredscience/2009/01/retailrobots/.
 H. Erzberger, “Automated conflict resolution for air traffic control,” 25 th International Congress of the Aeronautical Sciences, Jan 2006.
 B. P. Tice, “Unmanned aerial vehicles – the force multiplier of the 1990s,” Airpower Journal, Spring 1991.
 J. F. Fisac, M. Chen, C. J. Tomlin, and S. S. Sastry, “ReachAvoid Problems with TimeVarying Dynamics, Targets and Constraints.” in 18th International Conference on Hybrid Systems: Computation and Controls, 2015.
 O. Bokanowski, N. Forcadel, and H. Zidani, “Reachability and minimal times for state constrained nonlinear problems without any controllability assumption,” SIAM Journal on Control and …, pp. 1–24, 2010.
 E. Barron, “Differential Games with Maximum Cost,” Nonlinear analysis: Theory, methods & applications, pp. 971–989, 1990.
 J. Ding, J. Sprinkle, S. S. Sastry, and C. J. Tomlin, “Reachability calculations for automated aerial refueling,” in Proceedings of the 47th IEEE Conference on Decision and Control, Cancun, Mexico, 2008.
 M. Chen, Z. Zhou, and C. Tomlin, “Multiplayer reachavoid games,” in American Control Conference 2014, 2014.
 H. Huang, J. Ding, W. Zhang, and C. Tomlin, “A differential game approach to planning in adversarial scenarios: A case study on capturetheflag,” in Robotics and Automation (ICRA), 2011 IEEE International Conference on, 2011, pp. 1451–1456.
 I. Mitchell, A Toolbox of Level Set Methods, 2009, http://people.cs.ubc.ca/~mitchell/ToolboxLS/index.html.
 I. M. Mitchell, “The flexible, extensible and efficient toolbox of level set methods,” Journal of Scientific Computing, vol. 35, no. 23.
 S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces. SpringerVerlag, 2002, ISBN: 9780387954820.
 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. [Online]. Available: http://www.pnas.org/content/93/4/1591.abstract
 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, ser. HSCC ’11. New York, NY, USA: ACM, 2011, pp. 103–112. [Online]. Available: http://doi.acm.org/10.1145/1967701.1967718
 I. M. Mitchell and C. J. Tomlin, “Overapproximating reachable sets by hamiltonjacobi projections,” Journal of Scientific Computing, vol. 19, no. 13, pp. 323–346, 2003.
 J. McGrew, J. How, L. Bush, B. Williams, and N. Roy, “Air combat strategy using approximate dynamic programming,” AIAA Guidance, Navigation, and Control Conference, Aug 2008.
 J. B. Lasserre, D. Henrion, C. Prieur, and E. Trélat, “Nonlinear optimal control via occupation measures and lmirelaxations,” SIAM J. Control Optim., vol. 47, no. 4, pp. 1643–1666, Jun. 2008. [Online]. Available: http://dx.doi.org/10.1137/070685051
 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. 773–797, 1984.
 P. Varaiya, “On the existence of solutions to a differential game,” SIAM Journal on Control, vol. 5, no. 1, pp. 153–162, 1967.
 M. G. Crandall and P.L. Lions, “Viscosity solutions of HamiltonJacobi equations,” Transactions of the American Mathematical Society, vol. 277, no. 1, pp. 1–42, 1983.
 S. Osher and C.W. Shu, “Highorder essentially nonoscillatory schemes for HamiltonJacobi equations,” SIAM Journal on Numerical Analysis, vol. 28, no. 4, pp. pp. 907–922, 1991. [Online]. Available: http://www.jstor.org/stable/2157779
 S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces. Springer Verlag, 2003.
 C.W. Shu and S. Osher, “Efficient implementation of essentially nonoscillatory shockcapturing schemes,” Journal of Computational Physics, vol. 77, no. 2, pp. 439 – 471, 1988. [Online]. Available: http://www.sciencedirect.com/science/article/pii/0021999188901775
 D. L. Chopp, “Computing minimal surfaces via level set curvature flow,” Journal of Computational Physics, vol. 106, no. 1, pp. 77 – 91, 1993.
 M. Sussman, E. Fatemi, P. Smereka, and S. Osher, “An improved level set method for incompressible twophase flows,” Computers & Fluids, vol. 27, no. 5â6, pp. 663 – 680, 1998.