Efficient and Guaranteed Hamilton–Jacobi Reachability via Self-Contained Subsystem Decomposition and Admissible Control Sets

Hamilton-Jacobi reachability analysis is a useful tool for generating reachable sets and corresponding optimal control policies, but its use in high-dimensional systems is hindered by the “curse of dimensionality.” Self-contained subsystem decomposition is a proposed solution, but it can produce conservative or incorrect results due to the “leaking corner issue.” This issue arises from the inexact decomposition of the target set and inconsistencies across the computed control policies for each coupled subsystem. In this letter, we define and resolve this issue by introducing the notion of an admissible control set that enforces consistent control actions across the coupled subsystems. Our method efficiently computes exact reachable sets and the corresponding optimal control policy for self-contained subsystems with a decomposable goal (or failure) set. We also provide conservative under-approximations for goal (or failure) sets with inexact decomposition. In this conservative case, a local update method in the full dimensional space can be applied to recover exact results. We validate our approach on a 3D system and demonstrate its scalability on a 6D system.

Fig. 1.BRSs for a liveness problem with Dubins car model (7), where the goal cannot be decomposed exactly.Left: the BRS R * (t) (blue) from the original SCSD method incorrectly over-approximates the true BRS R(t) (brown).Right: the BRS R u (t) (yellow) from our method provides a conservative under-approximation while maintaining computational efficiency.
set (BRS) or backward reachable tube (BRT), which describes the states from which a system can be driven to reach a goal set (or avoid a failure set) despite disturbances.Among the various methods for computing a BRS/BRT [3], [4], Hamilton-Jacobi (HJ) formulations [5], [6] are particularly adept at handling general nonlinear dynamics, with flexible set representations facilitated through a grid-based approach.The most formidable challenge is the computational limits imposed by high-dimensional state spaces: the computational complexity grows exponentially as the dimensionality increases.This is referred to as the "curse of dimensionality" [7].
The SCSD methods provides exact results under certain conditions on both the system dynamics and problem formulation, with a significant reduction in computation time.It projects the goal or failure set to low-dimensional subspaces (of the original state space) and computes the corresponding BRSs in these subspaces.The high-dimensional BRS is reconstructed by back-projecting the low-dimensional BRSs and taking the intersection or union of them.
However, for some problem formulations, the BRSs from SCSD can be inconsistent with the BRS directly computed from the original higher-dimensional system, therefore loses the liveness or safety guarantees.Fig. 1, left, shows the BRS from the SCSD method (blue) incorrectly over-approximates the true BRS (brown).This issue is called the "leaking corner problem" and is primarily caused by: (1) the shared control between subsystems, and (2) the inexact decomposition of the goal (or failure) set.The same issue also appears in the BRTs.
The key insight of this letter is the admissible control set (ACS).Applying the ACS to subsystems drives all subsystems to reach their goal (or avoid the failure) set.In this letter, We focus on the single control case and leave the multiple control case as future work.Our main contributions are: • The formal definition of the "leaking corner" is provided, along with a method for solving this issue using the ACS.• A novel approach for identifying the ACS is presented, which we prove can be used to compute the exact BRS for the full system when the goal or failure set can be exactly decomposed.This solves the issue caused by shared control and produces exact results at greatly reduced computational cost.• For goal or failure sets that cannot be exactly decomposed onto the subspaces, our method can efficiently compute guaranteed under-approximations (for liveness formulations) and over-approximations (for safety formulations).For these conservative results, we provide an additional refinement step in the full-dimensional space that recovers the exact BRS.• The method is validated on a 3D 1 Dubins car example (as shown in yellow in Fig. 1, right).Scalability is shown by the 6D 1 Planar Vertical Take-Off and Landing (PVTOL) vehicle system [17].

II. BACKGROUND A. System Dynamics
We assume the dynamics are control-affine and are given by the following ordinary differential equation Assume t ≤ 0 is the initial time, s is the current time, z ∈ Z ⊂ R n is the state, and u ∈ U ⊂ R is the scalar control input, where U = {u : u min ≤ u ≤ u max }.Further assume the control signal u(•) : [t, 0] → U is drawn from the set of measurable functions U, i.e., u(•) is measurable.
Assume f : Z → R n and g : Z → R n are uniformly continuous, bounded, and Lipschitz continuous in z.As a result, given an initial state z, control u(•) ∈ U, we can solve a unique solution of (1), denoted as ζ(s; z, t, u(•)).

B. Backward Reachable Sets and Tubes
HJ reachability analysis is used to solve for reachable sets of various forms.This letter focuses on BRSs, though the methods introduced are generalizable to BRTs with [15,Proposition 4,Th. 3].Assume T is the target set.The maximal BRS R(t) (for liveness formulation) is the set of 1 "D" refers to the "state space dimension".initial states that can be driven into the goal set at the end of the time horizon: The minimal BRS A(t) (for safety formulation) is the set of initial states that must enter the failure set at the end of the time horizon, regardless of the control signal applied: It has been shown that HJ reachability can be solved as an optimal control problem.The loss function : R n → R of the optimal control problem is designed such that its zero sub-level set is the goal (or failure) set, generally called the target set: T = {z : (z) ≤ 0}.A common choice of the loss function is the signed distance function to the set T .
The value function associated with the maximal BRS for reaching a goal set is defined as If V R (z, t) ≤ 0, among all the control signals, there exists one control signal such that the final state of z along this signal is in the goal set.This shows V R (z, t) ≤ 0 ⇐⇒ z ∈ R(t), i.e., the zero sub-level set of V R (z, t) is the maximal BRS.
The value function associated with the minimal BRS for avoiding a failure set is ( The zero sub-level set of V A (z, t) is the minimal BRS, i.e., the set from which for all control signals the system will inevitably enter the failure set at the end of the time horizon.These value functions are viscosity solutions to the Hamilton-Jacobi-Bellman Partial Differential Equations (HJB-PDEs) [7], (6) where D t and D z represent the derivative with respect to t and z.The initial condition is V(z, 0) = (z).

C. Self-Contained Subsystem Decomposition
Numerically, the HJB-PDE ( 6) is solved on a discrete grid, therefore the computation scales exponentially with state dimension.This motivates the use of SCSD method to reduce computation cost.Decomposing coupled systems while maintaining exact global results is challenging.SCSD allows for the decomposition of a certain class of coupled nonlinear systems, defined below.
Definition 1 (SCSD): Consider the following special case z c are called "state partitions" of the system.Given system (1), the two subsystems of it are The two subsystems are coupled through the common state partition z c , and common control u, and we say the subsystems have "shared control."In this letter, we consider two subsystems, but the results generalize to an arbitrary number of subsystems.As an example, consider the 3D Dubins car: Self-contained subsystem trajectories ξ i (s; x i , t, u(•)) can evolve independently, and relate to the original system trajectory with the projection and back projection operator [15].
Definition 2 (Projection Operator): A projection operator proj i : R n → R n i +n c maps a state z in high-dimensional space R n to a state x i in low-dimensional space R n i +n c : Definition 3 (Back Projection Operator): The back projection operator proj With an abuse of notation, one can apply the projection and back projection operator to a set z ∈ T and the trajectory In the Dubins car example, consider a target given by T = {z : |p x | ≤ 1, |p y | ≤ 1}.Its projections onto the subspaces are Suppose the original system can be separated into p subsystems using SCSD.The BRS in full-dimensional system is: where S * (t) represents either the maximal or minimal BRS, without making an explicit distinction.The same process in (11) can also be applied to target set, which will later help us understand exact and inexact decomposition: T * contains all the states z with all of their projected states x i inside T i .Note that T * is not always equal to T , as illustrated in Fig. 2. In the case when T * = T , we say the goal (failure) set can be exactly decomposed.Additionally, this letter considers only the case where T * is the intersection (rather than union) of the back-projections; this is the standard problem formulation for decomposition.The difference between exact decomposition and inexact decomposition of a target set T .In both plots, the blue regions are target sets in the higher dimension (2D x-y plane), and the yellow regions are the sets of states in the lower dimensions (1D along x axis, and 1D along y axis.)

D. Notation of Reachable Sets
• The true BRSs for the full system dynamics are denoted as R(t) and A(t).• The subsystem BRSs are denoted as R i (t) and A i (t).
They are computed based on the projected targets T i .• The BRSs from the original SCSD method are denoted as R * (t) and A * (t), and are obtained through (11).• The BRSs for the method proposed in this letter (described in Section III) are denoted as R u (t) and A u (t).

III. LEAKING CORNERS IN EXACT DECOMPOSITION
We'll only consider the case of a full system being decomposed into two subsystems using the method in Definition 1.
The reconstructed BRS R * (t) from the original SCSD method may provide inexact approximations to the true BRS R(t).The inconsistent region is called the "leaking corner," whose elements are called "leaking corner states." Definition 4 (Leaking Corner): When the target set can be exactly decomposed, the "leaking corner issue" is caused only by the shared control.

A. Inconsistency With Shared Control
When the target set can be exactly decomposed, i.e., T = T * , the SCSD method generates an exact solution for the safety formulation, i.e., A * (t) = A(t): Lemma 1 [15, Th. 2]: In safety formulation, . Lemma 1 states that the SCSD for minimal BRS won't suffer from the "leaking corner issue" despite shared controls.
In the liveness formulation, if there exists no shared control between subsystems, the SCSD will also be free from the "leaking corner issue" [15].When the shared control exists, Lemma 2 shows that the SCSD provides over-approximation for the liveness formulation.
Lemma 2: In liveness formulation, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Randomly pick z ∈ R(t), and project it into the subsystems space: Combine Lemma 1 and 2, the BRS reconstructed from SCSD is an over-approximation of the true BRS.

B. Admissible Control
HJ reachability analysis provides the BRS and associated optimal control policy for the system and problem formulation.The "leaking corner issue" arises when the optimal control generated from each subsystem disagrees (e.g., one subsystem requires that the control takes on its minimum value and the other subsystem requires the maximum.)However, it is not always necessary to only use the optimal control.In many cases, a sub-optimal control signal can fulfill the tasks as well.This motivates us to define the set of "admissible controls" that accomplish the liveness (or safety) formulation in each subsystem.

Definition 5 (Admissible Control Signal Set) [ACSS]
: is the set of all the control signals that can drive the system to the goal/failure set To numerically compute the ACSS, we also define the admissible control set (ACS) for the discrete-time manner and denote it as U adms (z, t) := {u ∈ U : lim δt→0 + V(ζ (t; z, t − δt, u), t) ≤ 0}.ACS is the set of control inputs u adms ∈ U that drives the system to the "immediate" target (i.e., the zero sub-level set of the value function at the previous time step) at each time step.With first-order Taylor's expansion on the value function as we have scalar control input, ( 13) is a linear inequality and has the solution In the liveness formulation, only states inside the BRS have ACSS.But for the safety one, states outside of the BRS also have ACSS.

C. Computation of the Admissible Control Set and Associated BRS
The construction of the ACS for the full system requires the following steps: 1) Project the state to subsystem spaces using (8), compute the subsystems ACSs U adms,i (x i , t) using ( 14). 2) Take the intersection of ACSs in the high-dimensional system: U adms (z, t) = i∈{1,2} (U adms,i (proj i (z), t)).To construct the BRS from the ACS, trajectories from each initial state in the full-dimensional grid are generated with The corresponding value function is: Note that V u (z, t) is not unique because u adms (•) is chosen randomly within the U adms (z, t).Therefore, a state will have a different value for every run.Definition 6: The admissible BRS is the zero sub-level set of V u (z, t) S u is used here because this definition applies to both liveness and safety formulation.If the type of case needs to be emphasized, we use R u (t) for liveness formulation and use A u (t) for safety formulation.
The following theorem shows that Algorithm 1 solves the "leaking corner issue" caused by the shared control.
On the other hand, proves that z ∈ S(t) =⇒ z ∈ S u (t).With Definition 5, if z ∈ S(t), it must be true that U adms (z, t) = ∅.Then according to Eq. ( 15) and Definition 6, Fig. 3 shows the result on the Dubins car running example.The goal set can be exactly decomposed onto the subspaces: The computation in each subsystem is approximately 10 seconds.The updating process using Algorithm 1 is approximately 28 seconds.In comparison, the direct computation with the original fulldimensional system takes 600 seconds.

IV. COMPUTATION OF THE ADMISSIBLE BRS FOR
INEXACT TARGET DECOMPOSITION In most real-world applications, the target set T cannot be exactly decomposed, i.e., T = T * .In such cases, the approximate decomposition of the target set also causes "leaking corners".Our method also provides theoretically guaranteed and efficient over/under approximations for the BRSs.To further refine these approximations (when time and Fig. 3.The first plot shows the inconsistency between the SCSD BRS R * (t) (blue) and the true BRS R(t) (brown).The second plot shows the comparison with our method R u (t) (yellow).With a clearer view of one slice of all the BRSs in the third plot, it shows that the "leaking corners" (with cyan shade in plot 3) are eliminated.The difference between R u (t) and R(t) is caused by numerical error accumulated during the Algorithm 1.The quantitative comparison is in Table I. computation are available), we provide a method to correctly locate the leaking corners for local update methods.
First, we show that admissible BRSs from Algorithm 1 always produce optimistic over-approximations (for the liveness formulation) and conservative under-approximations (for the safety formulation) of the true BRS.
Theorem 2: The admissible BRSs have the following relationship to the true BRS as 1) R u (t) ⊆ R(t).2) A(t) ⊆ A u (t).In order to concisely determine the "leaking corners", an over-approximated BRS (for the liveness formulation) and an under-approximated BRS (for the safety formulation) are needed from the SCSD method.Then we show that with a certain approximation of the goal (or failure) set, the SCSD method from Definition 1 produces wanted results.
Theorem 3: The SCSD method will provide a conservative under (over) approximation of the true BRS: Proof: The proof follows from the proof of Lemma 2. 1) Note the fact that reaching a subset of a set guarantees reaching that set itself, i.e., T sub ⊆ T =⇒ R sub (t) ⊂ R(t) From Lemma 2, we have: 2) The proof is similar to the reachable case with T * ⊆ T =⇒ A * (t) ⊆ A(t), and is omitted here.
We illustrate Theorem 2 and 3 with the Dubins car example.The goal (failure) set is T = {z : (p 2 x + p 2 y ) ≤ 1} and therefore cannot be decomposed exactly onto the subspaces.Fig. 1 shows the liveness formulation, with the optimistic overapproximation of the target set in the SCSD method T 1 = {z : |p x | ≤ 1} and T 2 = {z : |p y | ≤ 1}.Fig. 4 shows the safety  I.

A. BRS Refinement
If an exact result is needed, we can apply a local recomputation method on uncertain states.The key problem is to find a region that contains all the leaking corners with the guarantee from Theorems 2 and 3.
In Theorem 2, the admissible BRS is a subset of the true BRS for the liveness formulation, thus we have With Definition 4, and Equation ( 16) and ( 17), The "leaking corner" that contains the boundary of the true BRS for the liveness formulation is: Following the same logic, the "leaking corner" containing the boundary of the true BRS for the safety formulation is: The region of computation in this space can be restricted to L(t) through our approach.The "leaking corners" can be locally updated with methods such as [18], which are beyond the scope of this letter.

V. NUMERICAL RESULTS
We first validate our method with the 3D Dubins car example, then demonstrate scalability on the 6D PVTOL system example.All simulations were executed on a platform with an Apple M1 chip with 4 performance and 4 efficiency cores, and 16 GB of RAM.The code can be found at https://github.com/ChongChongbst/Admissible-Control-Reconstruction.
For the 3D Dubins car, the result of the exact decomposition of the target set in the liveness formulation is shown in Fig. 3.The results for inexact decomposition of the target set are shown in Fig. 1 (for liveness) and Fig. 4 (for safety).
We use the Jaccard Index (JI) [19] to quantify the similarity of the sets over a common grid, and also quantify the percentage of falsely included (FI) and excluded (FE) states Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
These results are shown in Table I.For liveness formulation, an under-approximated result is more desirable to ensure the liveness of the system.However, for the safety formulation, only the over-approximated result can ensure safety reliably.
We test our method on the PVTOL system, which is computationally intractable with HelperOC because the dimension of its state space is too high.The scalability of our method is highlighted.The states are z = (p x , v x , p y , v y , θ, ω), and ṗx = v x , vx = −F sin θ, ṗy = v y , vy = F cos θ − g, θ = ω, ω = u, where the vectors p = (p x , p y ) and v = (v x , v y ) denote the position and velocity of the aircraft.g is the gravitational acceleration, θ represents the roll angle with respect to the horizon, and ω denotes the corresponding angular velocity.The control input u represents the rolling moment, which is bounded by [−2, 2].We specifically assume F is a constant thrust moment, and this system is being controlled over a relatively small time interval.In future work, we will reimplement this system where the F is an actual input.
With SCSD, the two subsystems states are x 1 = (p x , v x , θ, ω) and A 3D slice of the admissible BRS R u (t) at different times is shown in Fig. 5.The computation time is 2303 seconds.The result is restricted but reasonable to be analyzed under the general condition of this letter.

Fig. 2 .
Fig. 2.The difference between exact decomposition and inexact decomposition of a target set T .In both plots, the blue regions are target sets in the higher dimension (2D x-y plane), and the yellow regions are the sets of states in the lower dimensions (1D along x axis, and 1D along y axis.)

Fig. 4 .
Fig. 4. BRSs for the safety problem when the target set cannot be decomposed exactly.Left: the SCSD BRS R * (t) (blue) incorrectly under-approximates the true BRS R(t) (brown).Right: the admissible BRS R u (t) (yellow) is a conservative over-approximation.formulation with the same target, with the conservative underapproximation of the failure set in the SCSD method T 1 = {z : |p x | ≤ √ 0.5} and T 2 = {z : |p y | ≤ √ 0.5}.A quantitative comparison to the true BRSs is shown in TableI.

Fig. 5 .
Fig. 5.The 3D slice (in [p x , p y , θ] space) of the admissible BRS R u (t) at different time horizons for the 6D PVTOL system.

TABLE I QUANTITATIVE
MEASUREMENT OF (1) THE ADMISSIBLE BRS AND (2) THE BRS FROM SCSD, BOTH COMPARED TO A GLOBALLY COMPUTED BRS AS GROUND TRUTH