April 02, 2024

Control Lyapunov functions (CLFs) play a vital role in modern control applications, but finding them remains a problem. Recently, the control Lyapunov-value function (CLVF) and robust CLVF have been proposed as solutions for nonlinear time-invariant systems with bounded control and disturbance. However, the CLVF suffers from the “curse of dimensionality,” which hinders its application to practical high-dimensional systems. In this paper, we propose a method to decompose systems of a particular coupled nonlinear structure, in order to solve for the CLVF in each low-dimensional subsystem. We then reconstruct the full-dimensional CLVF and provide sufficient conditions for when this reconstruction is exact. Moreover, a point-wise optimal controller can be obtained using a quadratic program. We also show that when the exact reconstruction is impossible, the subsystems’ CLVFs and their “admissible control sets” can be used to generate a Lipschitz continuous CLF. We provide several numerical examples to validate the theory and show computational efficiency.

Ensuring stability is one of the most important tasks for autonomous systems operating in the real world. Control Lyapunov functions (CLFs) are energy-like functions that stabilize systems to their equilibrium points [1], [2]. However, it is well known that there lacks a “universal” way of finding CLFs for general nonlinear systems, especially with state or input constraints. Hand-designed and task-specific CLFs have been proposed [3]–[5].

A method to construct control Lyapunov-like functions for nonlinear systems with control and disturbance bounds has been proposed [6], [7]. These control Lyapunov value functions (CLVFs) are based on Hamilton-Jacobi reachability analysis [8], [9], and are constructed through dynamic programming. The CLVF guarantees exponential stabilizability, works for general nonlinear dynamics, and deals with state and input constraints well. However, it requires solving a CLVF variational inequality (CLVF-VI) on a discrete grid iteratively, therefore suffering from the “curse of dimensionality.” In the HJ reachability community, many works are proposed to solve this issue, including incorporating reinforcement learning [10], deep learning [11], self-contained subsystems decomposition (SCSD) [12], [13], warming starting [14] and model reduction [15], [16].

The SCSD method provides guaranteed exact solutions for a Hamilton-Jacobi reachability problem under certain assumptions on the coupled nonlinear system and the problem formulation. This approach was further generalized recently through the use of an “admissible control signal set” (ACSS) and “admissible control set” (ACS) [17], which are used to forward propagate and refine the reconstructed value function. This line of work has only been applied to reachability problems, and only provides guarantees on the level sets (reachable sets) of the computed value function, rather than guarantees on the value function itself.

In this work, we generalize the SDSC and ACS decomposition work to numerically compute CLVFs for high dimensional nonlinear systems with input constraints. The proposed method first applies the standard SCSD and finds the CLVFs for the subsystems. A value function in the original state space is then reconstructed using the subsystems’ CLVFs, and the ACS is used to determine the domain where this reconstruction is exact. This is demonstrated in Fig. 1. For systems in which exact reconstruction is not possible, the subsystems’ CLVFs and ACSs can be used to generate a Lipschitz continuous CLF and stabilize the original system. The improved scalability of the proposed method is validated in several numerical examples.

The main contributions of this work are:

We extend the definition of the ACSS and ACS to multi-input systems’ CLVFs.

We propose a method that combines the SCSD, ACS, and CLVF, to reconstruct a value function for high-dimensional systems whose CLF is not trivial to obtain,

We provide sufficient conditions on when the reconstruction is exact, i.e., the reconstructed value function is the CLVF for the original system.

For the case when the exact reconstruction is impossible, we show that a Lipschitz continuous CLF can be reconstructed using the subsystems’ CLVFs.

Consider the following nonlinear control-affine system: \[\begin{align} \label{eqn:dyn95sys} \frac{d x}{d s} = \dot{x}= f(x)+ g(x) \cdot u,x(0) = x_0,s\in [t,0], \end{align}\tag{1}\] where \(s\) is the time, \(x\in \mathcal{X}\subseteq \mathbb{R}^n\) is the state, and \(u\in \mathcal{U}\subseteq \mathbb{R}^m\) is the control input, where \(\mathcal{U}\) is a compact set. Assume \(f: \mathcal{X}\mapsto \mathbb{R}^n\) and \(g:\mathcal{X}\mapsto \mathbb{R}^n\) are uniformly continuous, bounded, and Lipschitz continuous in \(x\), and the origin is one equilibrium point, i.e., \(f(0)+g(0)\cdot 0 = 0\). Further assume the control signal \(u(\cdot)\) is drawn from measurable control functions: \[\begin{align} u(\cdot) \in \mathbb{U}:= \{ u: [t,0] \mapsto \mathcal{U}, u(\cdot) \text{ is measurable}\} . \end{align}\] Under these assumptions, we can solve for a unique solution of 1 , denoted as \(\xi(s;t,x,u(\cdot))\) (in short \(\xi(s)\)).

**Definition 1**. *Given the following time-varying CLVF (TV-CLVF) \[\begin{align} \label{eqn:tv95clvf} V_\gamma(x,t) = \min _{u\in \mathbb{U}} \max_{ s\in [t, 0]} e^{\gamma( s-t)}
\ell ( \xi(s) ) ,
\end{align}\qquad{(1)}\] the CLVF \(V_\gamma^\infty: D_\gamma \mapsto \mathbb{R}\) of 1 is the limit function: \[\begin{align}
\label{eqn:clvf} V_\gamma^\infty(x) = \lim _{t\rightarrow -\infty} V_\gamma(x,t).
\end{align}\qquad{(2)}\] Here, \(\gamma\) is a user-specified parameter that represents the desired exponential decay rate. The domain of the CLVF is given as \(D_\gamma \subseteq
\mathbb{R}^n\).*

Proposition 2 of [6] shows there exists \(u\in \mathcal{U}\) s.t. \[\begin{align} \label{eqn:Vdot} \dot{V}_\gamma^\infty (x) \leq \gamma V_\gamma^\infty(x), \quad \forall x\in \mathcal{D}_\gamma. \end{align}\tag{2}\]

Note that since we assume the system has an equilibrium point, the definition of CLVF is simplified compared to [6]. The CLVF value of a state captures the largest exponentially amplified deviation between the origin and the trajectory starting from that state. If this is finite, then there must exist some control signal for which the system trajectory converges to the origin at the same rate \(\gamma\). The domain of the CLVF \(D_\gamma\) is also known as the region of exponential stabilizability (ROES). A larger \(\gamma\) results in a faster convergence to the origin, but a smaller ROES.

Under the scope of this paper, the CLVF is a Lipschitz continuous CLF, while the main difference is that using CLVF, the user can trade off between a faster convergence or a larger ROES.

It has been shown that the TV-CLVF is the unique viscosity solution to the following VI, \[\begin{align} &0 = \max\biggl\{\ell(x) - V_\gamma(x, t), \notag \\ & D_tV_\gamma+ \min_{ u\in \mathcal{U}} D_xV_\gamma\cdot (f(x)+ g(x) \cdot u)+ \gamma V_\gamma(x,t) \biggl\}, \end{align}\] with terminal condition \(V_\gamma(x,0) = \ell(x)\). The CLVF can be obtained by solving this VI backward in time until convergence. Here, \(D_t\) denotes the time derivative, and \(D_x\) denotes the gradient with respect to \(x\).

**Definition 2**. *(SCSD) Given system 1 and assume there exists state partitions \(z_1 = (x_1,x_c) \in \mathcal{Z}_1\), \(z_2 = (x_2,x_c) \in
\mathcal{Z}_2\), where \(x_1\in\mathbb{R}^{n_1}\), \(x_2\in\mathbb{R}^{n_2}\), \(x_c\in\mathbb{R}^{n_c}\), \(n_1,n_2>0\), \(n_c\ge 0\), \(n_1+n_2+n_c=n\). Assume also the control inputs can be partitioned similarly with \(v_1 = (u_1,u_c) \in
\mathcal{V}_1\), \(v_2 = (u_2,u_c) \in \mathcal{V}_2\), where \(u_1 \in \mathbb{R}^{m_1}\), \(u_2 \in \mathbb{R}^{m_2}\), \(u_c \in \mathbb{R}^{m_c}\) and \(m_1+m_2+m_c = m\). The two self-contained subsystems of 1 are \[\begin{align}
\label{eqn:subsys} \dot{z}_1 = f_1(z_1)+g_1(z_1)v_1, \quad \dot{z}_2 = f_2(z_2)+g_2(z_2)v_2,
\end{align}\qquad{(3)}\] with corresponding solution \(\phi_1 (s)\), \(\phi_2 (s)\).*

When \(x_c\) or/and \(u_c\) exist, we say the subsystems are coupled through the common states \(x_c\) or/and controls \(u_c\), and they have shared states or/and controls. The shared controls are one of the main causes of the inconsistency between the reconstructed value function from the subsystems and the original value function.

We introduce the projection and back projection operator on both state and control.

**Definition 3**. *(Projection and back projection operators.)
A state projection operator \(\mathcal{P}_{x, i}: \mathbb{R}^n \mapsto \mathbb{R}^{n_i+n_c}\) maps a state \(x\in \mathbb{R}^n\) to a state \(z_i
\in \mathbb{R}^{n_i+n_c}\): \[\begin{align} \mathcal{P}_{x,i} (x) := (x_i,x_c) = z_i.
\end{align}\] A control projection operator \(\mathcal{P}_{u, i}: \mathbb{R}^m \mapsto \mathbb{R}^{m_i+m_c}\) maps a control \(u\in \mathbb{R}^m\) to a control \(v_i \in \mathbb{R}^{m_i+m_c}\): \[\begin{align} \mathcal{P}_{u,i} (u) := (u_i,u_c) = v_i.
\end{align}\] A state back projection operator \(\mathcal{P}^{-1}_{x,i}: \mathbb{R}^{n_i+n_c} \mapsto \mathbb{R}^{n}\) maps a state \(z_i\in \mathbb{R}^{n_i+n_c}\) to
a set of states \(x\in \mathbb{R}^{n}\): \[\begin{align}
\mathcal{P}^{-1}_{x,i}(z_i) :=\{x\in \mathcal{X}: (x_i,x_c)=z_i\}.
\end{align}\] A control back projection operator \(\mathcal{P}^{-1}_{u,i}: \mathbb{R}^{m_i+m_c} \mapsto \mathbb{R}^{m}\) maps a control \(v_i\in
\mathbb{R}^{m_i+m_c}\) to a set of controls \(u\in \mathbb{R}^{m}\): \[\begin{align}
\mathcal{P}^{-1}_{u,i}(v_i) :=\{u\in \mathcal{U}: (u_i,u_c)=v_i\}.
\end{align}\]*

Following [17], we define the ACSS for the TV-CLVF.

**Definition 4**. *The ACSS of ?? is the set of all the control signals such that the corresponding trajectory achieves the same value \(V_\gamma(x,t)\):*

*\[\begin{align}
\mathbb{U}_\text{a}(x,t) = \{& u(\cdot)\in \mathbb{U}: \\
&\max_{ s\in [t, 0]} e^{\gamma( s-t)} \ell ( \xi(s;t,x,u(\cdot)) ) = V_\gamma(x,t) \}
\end{align}\]*

For numerical computation, define the ACS given a time step \(\delta t\) \[\begin{align} \label{eqn:32ACS} \mathcal{U}_\text{a}(x,t) = &\{ u\in \mathcal{U}: V_\gamma(x,t-\delta t) - V_\gamma(x,t) - \notag \\ & \bigl( D_xV_\gamma\cdot (f(x)+ g(x) \cdot u)+ \gamma V_\gamma(x,t) \bigl) \delta t\geq 0 \}. \end{align}\tag{3}\]

One admissible control signal \(u_a(\cdot)\) for \(x\) can be computed by concatenating the ACSs over time: \(u_a(\cdot) = [\mathcal{U}_\text{a}(x,t), \mathcal{U}_\text{a}(\xi(s_1),s_1),...,\mathcal{U}_\text{a}(\xi(0),0)]\). Both the ACSS and ACS are guaranteed to be non-empty [17], and they can be similarly defined for subsystems.

In this section, we begin with an idealized case, where the system can be decomposed into two subsystems with no shared control, and show the reconstructed value function is the CLVF for the original system. This is called *exact reconstruction*.
We then provide a sufficient condition of exact reconstruction for cases where subsystems are coupled through shared controls. Further, for systems where an exact reconstruction is impossible, we apply the concept of ACS, and show how the reconstructed
value function is a Lipschitz continuous CLF. All theorems follow naturally to the cases where the original system is decomposed into more than two subsystems.

Denote the TV-CLVF of the two subsystems as \[\begin{align} &V_{\gamma,i} (z_i,t) = \min _{v_i \in \mathbb{V}_i} \max_{ s\in [t, 0]} e^{\gamma( s-t)} \ell_i ( \phi_i (s) ) , \quad i = 1,2,
\end{align}\] and the CLVF for subsystems as \(V^\infty_{\gamma,i} : \mathcal{D}_{\gamma,i} \subseteq \mathcal{Z}_i \mapsto \mathbb{R}\), \[\begin{align} V^\infty_{\gamma,i}(z_i) =
\lim_{t\rightarrow -\infty } V_{\gamma,i} (z_i,t), \quad i = 1,2.
\end{align}\] For clarity, the CLVF computed directly for the high-dimensional system 1 is called the *original CLVF*, while the CLVF (TV-CLVF) reconstructed from subsystems is called the *reconstructed CLVF
(TV-CLVF)*, denoted as \(\bar{V}^\infty_\gamma(x)\) (\(\bar{V}_\gamma(x,t)\)). Using 3 , the ACS can be solved:
\[\begin{align} \label{eqn:32ACS95halfspace} D_xV_\gamma(x,t) g(x) \cdot u&\leq \frac{V_\gamma(x,t-\delta t) - V_\gamma(x,t)}{\delta t} \notag \\ &-\gamma V_\gamma(x,t) -
D_xV_\gamma(x,t) f(x) ,
\end{align}\tag{4}\] which is a linear inequality for single-input systems, and a half-space for multi-input systems. Though we cannot get an analytic solution of this half-space, it is numerically easy to obtain one element of it.

In this section, we assume \(\ell(x) = \max (\ell_1 (z_1),\ell_2 (z_2))\), meaning the loss function for the original system equals the reconstructed loss function from subsystems. This is possible, because [7] relaxes the choice of the loss to be any vector norm. One possible choice is \(\ell_i(z_i) = ||z_i||_\infty\).

**Lemma 1**. *Assume there are no common controls after decomposition. Then, \(\bar{V}^\infty_\gamma(x) = \max (V^\infty_{\gamma,1} (z_1) ,V^\infty_{\gamma,2} (z_2)) = V_\gamma^\infty(x)\) , where \(z_1 = \mathcal{P}_{x,1} (x)\), \(z_2 = \mathcal{P}_{x,2} (x)\), and \(\mathcal{D}_\gamma = \mathcal{P}^{-1}_{x,1}(\mathcal{D}_{\gamma,1}) \bigcap
\mathcal{P}^{-1}_{x,2}(\mathcal{D}_{\gamma,2}).\)*

**Proof.* Define \(\bar{V}_\gamma(x, t) = \max ( V_{\gamma,1} (z_1,t),V_{\gamma,2} (z_2,t) )\), we first show \(\bar{V}_\gamma(x, t) = V_\gamma(x,t)\) by contradiction.*

*Assume \(\bar{V}_\gamma(x, t) < V_\gamma(x,t)\), then there exists optimal control signals \(u^*(\cdot) \in \mathbb{U}\), \(v_1^*(\cdot) \in
\mathbb{V}_1\) and \(v_2^*(\cdot) \in \mathbb{V}_2\), s.t. \[\begin{align} \label{lemma1:eqn1} &\max_{ s\in [t, 0]} e^{\gamma( s-t)} \ell ( \xi(s;t,
x,u^*(\cdot) ) ) \notag \\ >& \max \bigl( \max_{ s\in [t, 0]} e^{\gamma( s-t)} \ell_1 ( \phi_1 (s;t, z_1,v_1^*(\cdot) ) ) ,\notag \\ &\max_{ s\in [t, 0]} e^{\gamma( s-t)} \ell_2 ( \phi_2 (s;t, z_2,v_2^*(\cdot) ) ) \bigl) \notag\\ = & \max_{
s\in [t, 0]} \max \bigl( e^{\gamma( s-t)} \ell_1 ( \phi_1 (s;t, z_1,v_1^*(\cdot) ) ) , \notag\\ &e^{\gamma( s-t)} \ell_2 ( \phi_2 (s;t, z_2,v_2^*(\cdot) ) ) \bigl)\notag \\ & = \max_{ s\in [t, 0]} e^{\gamma( s-t)} \max \bigl( \ell_1 ( \phi_1 (s;t,
z_1,v_1^*(\cdot) ) ) ,\notag \\ &\ell_2 ( \phi_2 (s;t, z_2,v_2^*(\cdot) ) ) \bigl).
\end{align}\qquad{(4)}\] Since the two subsystems have no shared controls, we could reconstruct one trajectory for 1 using the back projection: \(\bar \xi(s) = \mathcal{P}^{-1}_{x,1}
(\phi_1 (s)) \bigcap \mathcal{P}^{-1}_{x,2}(\phi_2 (s))\). Since the shared states evolve independent of the controls applied, ?? becomes \[\begin{align} &\max_{ s\in [t, 0]} e^{\gamma( s-t)} \ell ( \xi(s;t,
x,u^*(\cdot) ) ) \\ > & \max_{ s\in [t, 0]} e^{\gamma( s-t)} \max \bigl( \ell_1 ( \phi_1 (s;t, z_1,v_1^*(\cdot) ) ) , \\ &\ell_2 ( \phi_2 (s;t, z_2,v_2^*(\cdot) ) ) \bigl)\\ = & \max_{ s\in [t, 0]} e^{\gamma( s-t)} \ell ( \bar \xi(s;t,
x,u^*(\cdot) ) ).
\end{align}\] This means \(u^*(\cdot)\) is not optimal, which is a contradiction.*

*Similarly, assume \(\bar{V}_\gamma(x, t) > V_\gamma(x,t)\), then there also exists optimal control signals \(u^*(\cdot) \in \mathbb{U}\), \(v_1^*(\cdot)
\in \mathbb{V}_1\) and \(v_2^*(\cdot) \in \mathbb{V}_2\), s.t. \[\begin{align} \label{lemma1:eqn2} &\max_{ s\in [t, 0]} e^{\gamma( s-t)} \ell ( \xi(s;t,
x,u^*(\cdot) ) ) \notag \\ < & \max_{ s\in [t, 0]} e^{\gamma( s-t)} \max \bigl( \ell_1 ( \phi_1 (s;t, z_1,v_1^*(\cdot) ) ) , \notag\\ &\ell_2 ( \phi_2 (s;t, z_2,v_2^*(\cdot) ) ) \bigl).
\end{align}\qquad{(5)}\] This is obtained following the same procedure as ?? . However, projecting \(\xi(s)\) into the two subsystems spaces, we get \(\bar \phi_1 (s) =
\mathcal{P}_{x,1} (\xi(s))\), \(\bar \phi_2 (s) = \mathcal{P}_{x,2} (\xi(s))\), and we have \[\begin{align} &\max_{ s\in [t, 0]} e^{\gamma( s-t)} \ell ( \xi(s;t, x,u^*(\cdot) ) )\\ =
& \max_{ s\in [t, 0]} e^{\gamma( s-t)} \max \bigl( \ell_1 ( \bar \phi_1 (s)), \ell_2 ( \bar \phi_2 (s) ) \bigl).
\end{align}\] Using ?? , we have \[\begin{align} & \max_{ s\in [t, 0]} e^{\gamma( s-t)} \max \bigl( \ell_1 ( \bar \phi_1 (s)), \ell_2 ( \bar \phi_2 (s) ) \bigl)\\ < & \max_{ s\in [t, 0]} e^{\gamma( s-t)} \max
\bigl( \ell_1 ( \phi_1 (s;t, z_1,v_1^*(\cdot) ) ) , \notag\\ &\ell_2 ( \phi_2 (s;t, z_2,v_2^*(\cdot) ) ) \bigl),
\end{align}\] which shows that either \(v_1^*\) or \(v_2^*\) is not optimal. Combined, we have shown that \(\forall t \leq 0\), \(x\in \mathbb{R}^n\), \(V_\gamma(x,t) = \bar{V}_\gamma(x,t)\). Next we show that \(\bar {V}_\gamma^\infty(x) = V_\gamma^\infty(x)\).*

*Since the CLVF exists for two subsystems, we have \[\begin{align} V_\gamma^\infty(x) & = \lim_{t\rightarrow -\infty} V_\gamma(x,t)\\ & =\lim_{t\rightarrow -\infty} \bar{V}_\gamma(x, t) \\ & = \lim_{t\rightarrow
-\infty} \max \bigl( V_{\gamma,1} (z_1,t),V_{\gamma,2} (z_2,t) \bigl) \\ & = \max \bigl( \lim_{t\rightarrow -\infty} V_{\gamma,1} (z_1,t), \lim_{t\rightarrow -\infty} V_{\gamma,2} (z_2,t) \bigl) \\ & = \max \bigl( V^\infty_{\gamma,1} (z_1),
V^\infty_{\gamma,2} (z_2) \bigl)= V_\gamma(x).
\end{align}\] This also shows that \(\mathcal{D}_\gamma \supseteq \mathcal{P}^{-1}_{x,2}(\mathcal{D}_{\gamma,1}) \bigcap \mathcal{P}^{-1}_{x,2}(\mathcal{D}_{\gamma,2})\).*

*Assume \(\exists x\notin \mathcal{D}_\gamma\), \(x\in \mathcal{P}^{-1}_{x,2}(\mathcal{D}_{\gamma,1}) \bigcap \mathcal{P}^{-1}_{x,2}(\mathcal{D}_{\gamma,2})\). This means the optimal
trajectory \(\xi\) starting from \(x\) satisfies the following two condition simultaneously: 1) \(V_\gamma^\infty(x)\) is infinite and 2) both \(V^\infty_{\gamma,1} (z_1)\) and \(V^\infty_{\gamma,2} (z_2)\) are finite. This is a contradiction. Therefore, we have \(\mathcal{D}_\gamma =
\mathcal{P}^{-1}_{x,2}(\mathcal{D}_{\gamma,1}) \bigcap \mathcal{P}^{-1}_{x,2}(\mathcal{D}_{\gamma,2})\). ◻*

Lemma 1 shows that if there exists no shared control, the reconstructed value function is the CLVF, and their domains are also the same. This means system 1 can be exponentially stabilized to the origin from \(\mathcal{D}_\gamma\).

**Remark 2**. *Numerically, Lemma 1 is easy to implement: after decomposition, it is easy to check if all subsystems have
shared controls. If Lemma 1 is applicable, we simply compute the CLVF for the subsystems and take the max to reconstruct the CLVF for the
original system.*

The essence of this proof is that the optimal control signal for the original system can always be reconstructed from the subsystems’ optimal control signals, i.e. \(u^*(\cdot) = [v_1^*(\cdot);v_2^*(\cdot)]\), and that the projection and back projection between the original system’s trajectory and subsystems’ trajectories are both unique. This is guaranteed because there is no shared control. However, many practical systems cannot be decomposed into fully independent subsystems. From Remark 1, if there are shared controls, there must also be shared states. If two subsystems require conflicting values of the shared control(s), the back projection of subsystems’ trajectories cannot reconstruct the original system’s trajectory.

For the subsystems, denote the ACSS as \(\mathbb{U}_{\text{a},1}(z_1,t)\) and \(\mathbb{U}_{\text{a},2}(z_2,t)\), and their shared control component as \(\mathbb{U}_{\text{a},1}^c (z_1,t)\) and \(\mathbb{U}_{\text{a},2}^c (z_2,t)\). The ACSS for 1 can be reconstructed as \[\begin{align} &\bar {\mathbb{U}}_\text{a}(x,t) = \mathcal{P}^{-1}_{u,1}(\mathbb{U}_{\text{a},1}(\mathcal{P}_{x,1} (x),t) ) \bigcap \mathcal{P}^{-1}_{u,2} ( \mathbb{U}_{\text{a},2}(\mathcal{P}_{x,2} (x),t) )\\ &\bar {\mathbb{U}}_\text{a}^{c}(x,t) = \mathbb{U}_{\text{a},1}^c(\mathcal{P}_{x,1} (x),t) \bigcap \mathbb{U}_{\text{a},2}^c(\mathcal{P}_{x,2} (x),t). \end{align}\] Note that though \(\mathbb{U}_{\text{a},1}^c (z_1,t)\) and \(\mathbb{U}_{\text{a},2}^c (z_2,t)\) are non-empty, \(\mathbb{U}_\text{a}^{c}(x,t)\) might be empty for some \(x\).

**Remark 3**. *One can see that \(\bar {\mathbb{U}}_\text{a}(x,t)\) is non-empty if and only if \(\bar {\mathbb{U}}_\text{a}^c(x,t)\) is non-empty.*

Now, we provide one sufficient condition on exact reconstruction for the case with shared controls.

**Theorem 1**. *Let \(\bar{V}^\infty_\gamma(x) = \max (V^\infty_{\gamma,1} (z_1) ,V^\infty_{\gamma,2} (z_2))\) where \(z_1 = \mathcal{P}_{x,1} (x)\), \(z_2 = \mathcal{P}_{x,2} (x)\). Assume \(\bar {\mathbb{U}}_\text{a}^{c}(x,t)\) is non-empty for all \(t \leq 0\) and \(x\in
\mathcal{S}_{\gamma}\), where \(\mathcal{S}_{\gamma}\) is some level set of \(\bar{V}^\infty_\gamma(x)\). Then, \(\bar{V}^\infty_\gamma(x) =
V_\gamma^\infty(x)\) on \(\mathcal{S}_\gamma\).*

**Proof.* The proof can be obtained similarly to Lemma 1. If \(\bar {\mathbb{U}}_\text{a}^c (x,t)\) is non-empty for all \(t\leq 0\), \(x\in \mathcal{S}_\gamma\), then given the subsystems’ trajectories \(\phi_1 (s)\) and \(\phi_2 (s)\), a trajectory of 1 can again be reconstructed as \[\begin{align} \bar \xi(s) =\mathcal{P}^{-1}_{x,1}(\phi_1 (s)) \bigcap \mathcal{P}^{-1}_{x,2} (\phi_2 (s)). \end{align}\] Similarly, given a trajectory \(\xi(s)\) of 1 , the two subsystems’ trajectories can be obtained as \[\begin{align} \bar \phi_1 (s) = \mathcal{P}_{x,1}(\xi(s)), \quad \bar \phi_2 (s) = \mathcal{P}_{x,1}(\xi(s)). \end{align}\] Therefore, assume \(\bar{V}_\gamma(x, t) < V_\gamma(x,t)\), we can obtain ?? , and assume \(\bar{V}_\gamma(x, t) > V_\gamma(x,t)\), we can obtain ?? . In other words, \(\forall x\in \mathcal{S}_\gamma\), we again construct two contradictions to show that \(\bar{V}_\gamma(x, t) < V_\gamma(x,t)\) and \(\bar{V}_\gamma(x, t) > V_\gamma(x,t)\) cannot happen. The remaining steps are identical to the proof of Lemma 1. ◻*

A direct implication of Theorem 1 is that the reconstructed ACSS (if not empty) is always equal to the ACSS directly computed for 1 . In other words, if \(\bar {\mathbb{U}}_\text{a}^{c}(x,t)\) is non-empty for all \(t \leq 0\) and \(x\in \mathcal{S}_{\gamma}\), then \(\bar {\mathbb{U}}_\text{a}^c(x,t) = \mathbb{U}_\text{a}^c(x,t)\). We provide the following standard results without proof.

**Proposition 1**. *System 1 can be exponentially stabilized to the origin from \({\mathcal{S}}_\gamma\).*

The numerical implementation of Theorem 1 is shown in Alg. 2. The convergence time for the subsystems may differ, therefore, after one TV-CLVF converges, we repeat its final value until the other subsystem’s CLVF converges. Further, since the TV-CLVF is computed backward in time, we reverse the time vector of the TV-CLVF (\(V_\gamma(x,t) = V_\gamma(x, T - t)\) when used to forward propagate the trajectory.

Alg. 2 has a maximum iteration \(\frac{T}{\delta t}\). It first computes the ACS from the subsystem’s value function and reconstructs the ACS for 1 . Then, it checks if the shared control component of the ACS is empty. If not empty, it applies one random control input from ACS, updates the state, and repeats. For any initial state whose ACS is non-empty along the whole trajectory, concatenating the control used, we get one admissible control signal. This state is added to a set \(\mathcal{T}_\gamma\), and \(\mathcal{S}_\gamma\) is the largest level set of \(\bar{V}^\infty_\gamma(x)\) contained in \(\mathcal{T}_\gamma\). Note that we take the maximum between \(V^\infty_{\gamma,1} (z_1,0)\) and \(V^\infty_{\gamma,2} (z_2,0)\). This is because we reverse the time vector at line 3.

Lemma 1 and Theorem 1 both provide exact reconstructions of the CLVF, but on different domains. In this section, we generalize to the case in which the subsystems have shared control, and therefore the sufficient condition in Theorem 1 is not met. Though exact reconstruction of a CLVF in the original state space is impossible, we can still reconstruct a Lipschitz continuous CLF and it ROES using the subsystems’ CLVFs and ACSs. We remove the assumption that \(\ell(x) = \max (\ell_1 (z_1),\ell_2 (z_2))\).

A key difference compared to the exact reconstruction is that we focus on the \(\mathcal{U}_\text{a}(x)\) for \(V_\gamma^\infty(x)\), instead of \(\mathbb{U}_\text{a}(x,t)\) for \(V_\gamma(x,t)\). That is to say, we do not care about how the subsystems TV-CLVF evolve, but only care about its CLVF. Note that for the CLVF, its ACS becomes time-independent: \[\begin{align} \label{eqn:32ACS95clvf} \mathcal{U}_\text{a}(x) = \{ u\in \mathcal{U}: D_xV_\gamma^\infty[f(x)+ g(x) \cdot u] \leq -\gamma V_\gamma^\infty(x) \}. \end{align}\tag{5}\] Denote the ACSs for two subsystems as \(\mathcal{U}_{\text{a},1}(z_1)\), \(\mathcal{U}_{\text{a},2}(z_2)\), and their shared control components as \(\mathcal{U}_{\text{a},1}^c (z_1)\) and \(\mathcal{U}_{\text{a},2}^c (z_2)\). The ACS for 1 can be reconstructed as \[\begin{align} &\bar {\mathcal{U}}_\text{a}(x) = \mathcal{P}^{-1}_{u,1}(\mathcal{U}_{\text{a},1}(\mathcal{P}_{x,1} (x)) ) \bigcap \mathcal{P}^{-1}_{u,2} ( \mathcal{U}_{\text{a},2}(\mathcal{P}_{x,2} (x)) )\\ &\bar {\mathcal{U}}_\text{a}^{c}(x) = \mathcal{U}_{\text{a},1}^c(\mathcal{P}_{x,1} (x)) \bigcap \mathcal{U}_{\text{a},2}^c(\mathcal{P}_{x,2} (x)). \end{align}\]

**Theorem 2**. *Let \(\bar{V}^\infty_\gamma(x) = V^\infty_{\gamma,1} (z_1) + V^\infty_{\gamma,2} (z_2)\) where \(z_1 = \mathcal{P}_{x,1} (x)\), \(z_2 = \mathcal{P}_{x,2} (x)\). Assume \(\bar {\mathcal{U}}_\text{a}^{c}(x)\) is non-empty for all \(x\in \bar {\mathcal{S}}_{\gamma}\), where \(\bar {\mathcal{S}}_{\gamma}\) is some level set of \(\bar{V}^\infty_\gamma\). Then, \(\bar{V}^\infty_\gamma\) is a Lipschitz continuous local CLF for 1 on \(\bar {\mathcal{S}}_\gamma\).*

**Proof.* From 2 , there exist \(v_1 \in \mathcal{U}_{\text{a},1}(z_1)\) and \(v_2 \in \mathcal{U}_{\text{a},2}(z_2)\), such that \[\begin{align} &D_{x} \bar{V}^\infty_\gamma(x) \cdot (f(x)+ g(x) \cdot u)\\ =& D_{x} V^\infty_{\gamma,1} (z_1) \cdot ( f_1(z_1)+g_1(z_1)v_1) + \\ &D_{x} V^\infty_{\gamma,2} (z_2) \cdot ( f_2(z_2)+g_2(z_2)v_2) \\ \leq &-\gamma V^\infty_{\gamma,1} (z_1) -\gamma V^\infty_{\gamma,2} (z_2) \\ =& -\gamma \bar{V}^\infty_\gamma(x), \end{align}\] where \(u= \mathcal{P}^{-1}_{u,1}(v_1) \bigcap \mathcal{P}^{-1}_{u,2} (v_2 )\). In other words, \(\forall x\in \bar {\mathcal{S}}_\gamma\), there exists some \(u\in \mathbb{U}\) s.t. the Lie derivative along 1 is smaller than \(-\gamma \bar{V}^\infty_\gamma(x)\).*

*Further, since both \(V^\infty_{\gamma,1}\) and \(V^\infty_{\gamma,2}\) are positive definite and all level sets are closed, \(\bar{V}^\infty_\gamma(x) =
V^\infty_{\gamma,1}+V^\infty_{\gamma,2}\) must also be positive definite and have closed level sets. We conclude that \(\bar{V}^\infty_\gamma(x)\) is a local CLF. ◻*

Similar to Proposition 1, system 1 can be exponentially stabilized to the origin from \(\bar {\mathcal{S}}_\gamma\).

**Remark 4**. *We use the summation instead of the maximum to reconstruct the value function. This is because using summation does not introduce additional non-differentiable points, whereas for the maximum, the reconstructed value
function is non-differentiable for all \(x\) where \(V^\infty_{\gamma,1} (z_1) = V^\infty_{\gamma,2} (z_2)\).*

For both exact and inexact reconstruction, the controller that guarantees exponential stabilizability can be synthesized by solving the following QP, \[\begin{align} \label{eqn:QP} &\min_{u\in \mathcal{U}} ||u- u_r|| \notag \\ &\text{ s.t. } D_x\bar{V}^\infty_\gamma[f(x)+ g(x) \cdot u] \leq -\gamma \bar{V}^\infty_\gamma(x) . \end{align}\tag{6}\] This QP is guaranteed to be feasible on \(\mathcal{D}_\gamma\) (or \({\mathcal{S}}_\gamma\)) [6]. For systems whose CLVFs can be exactly reconstructed, the constraints in 6 can either be solved directly from the reconstructed CLVF, or from subsystems CLVFs. For systems in which the CLVFs cannot be exactly reconstructed, the control can be determined by solving 6 with the reconstructed CLF.

In this section, we provide numerical examples that validate our theory. This includes examples in 2D, 3D, and 10D. For simplicity, in all examples, the loss function is chosen to be the infinity norm. All simulations are conducted using MATLAB and tool boxes [18], [19].

Consider the following 2D system \[\begin{align} \label{eqn:32nonlinear2D} \dot{x}_1 = x_1^2+u_1, \quad \dot{x}_2 = x_2+u_2, \end{align}\tag{7}\] where \(u_1 \in [-4,4]\) and \(u_2 \in [-1,1]\). The two states evolve purely on their own, therefore each state’s dynamics is one subsystem. We compute the CLVF for both subsystems and also directly for the original system with \(\gamma = 0.1\) and \(0.3\). The results are shown in Fig. 1.

Consider the following 3D system \[\begin{align} \label{eqn:32multlinear3D} \dot{x_1} = x_3 + u_1,\dot{x_2} = x_3 + u_2,\dot{x_3} = u_3, \end{align}\tag{8}\] where \(u_1, u_2\in[-1,1]\) and \(u_3 \in [-0.5,0.5]\) are the control inputs. Take \(z_1 = [x_1,x_3]\), \(z_2 = [x_2,x_3]\), and \(v_1 = [u_1,u_3]\), \(v_2 = [u_2,u_3]\), we can decompose the system into two 2D subsystems. We can apply Alg. 2 to find the domain \(S_\gamma\), such that if initialized inside, the system will stabilize to its equilibrium point. The results are shown in Fig. 3.

Consider the following 3D system \[\begin{align} \label{eqn:32linear3D} \dot{x_1} = x_3 ,\dot{x_2} = x_3,\dot{x_3} = u, \end{align}\tag{9}\] where \(u_1\in[-1,1]\) is the control inputs. Take \(z_1 = [x_1,x_3]\), \(z_2 = [x_2,x_3]\), and \(v_1 = [u_1,u_3]\), \(v_2 = [u_2,u_3]\), we can decompose the system into two 2D subsystems. Although reconstruction to the exact CLVF is not possible, we can still validate the reconstructed lipschitz continuous CLF through a trajectory simulation inside \(S_\gamma\). The results are shown in Fig. 4.

Consider the 10D quadrotor system: \[\begin{align} &\dot{x}_1 = x_2 ,\dot{x}_2 = g\tan{x_3},\dot{x_3} = -d_1 x_3 + x_4, \notag \\ &\dot{x}_4 = -d_0 x_3 + n_0 u_1,\dot{x}_5 = x_6 ,\dot{x}_6 = g\tan{x_7}, \notag \\ & \dot{x}_7 = -d_1x_7 + x_8,\dot{x}_8 = -d_0x_7 + n_0u_2, \notag \\ \; & \dot{x}_9 = x_{10} ,\dot{x}_{10} = u_3, \label{eq:3210D95Quad} \end{align}\tag{10}\] where \((x,y,z)\) denote the position, \((x_2, x_6, x_{10})\) denote the velocity, \((x_3, x_7)\) denote the pitch and roll, \((x_4, x_8)\) denote the pitch and roll rates, and \((u_1, u_2, u_3)\) are the controls. The system parameters are set to be \(d_0 = 10, d_1 = 8, n_0 = 10, k_T = 0.91, g = 9.81\), \(|u_1|, |u_2| \leq \pi/4\), \(u_3\in [-1, 1]\).

This 10D system can be decomposed into three subsystems: \(X\) subsystem with states \([x_1,x_2,x_3,x_4]\), \(Y\) subsystem with states \([y,x_6,x_7,x_8]\), and \(Z\) subsystem with states \([x_9,x_{10}]\). It can be verified that all three subsystems have an equilibrium point at the origin. Further, there’s no shared control or states among subsystems, therefore, Lemma 1 can be used to exactly reconstruct the CLVF using \(\bar \ell(x) = ||x||_\infty\). The result is shown in Fig. 5.

In this paper, we presented an exact and inexact reconstruction method for CLVFs by leveraging ACS, providing a corresponding QP controller, and identifying the domain that can exponentially stabilize the system. We validate our method with a nonlinear 2D system with no shared control for exact reconstruction, as well as a multi-input linear 3D system for the shared control case. Further, we validate the inexact reconstruction method with the single-input linear 3D system by finding the corresponding domain. Lastly, a 10D quadrotor example is provided, showing numerical efficiency.

Through this method, we can scalably compute higher-dimensional CLVFs for a class of nonlinear systems. Future work includes validating high-dimensional CLVFs through neural solvers like DeepReach [11], applying this method to online trajectory planning problems, finding the “smallest control invariant set” defined in the original CLVF work [6], and exploring more on the region where the reconstructed ACS is empty.

[1]

E. D. Sontag, “A ‘universal’ construction of Artstein’s theorem on nonlinear stabilization,” *Systems & control letters*, 1989.

[2]

R. A. Freeman and J. A. Primbs, “Control lyapunov functions: New ideas from an old source,” in *Conf. on Decision and Control*, 1996.

[3]

F. Camilli, L. Grüne, and F. Wirth, “Control Lyapunov functions and Zubov’s method,” *SIAM Journal on
Control and Optimization*, 2008.

[4]

P. Ogren, M. Egerstedt, and X. Hu, “A control lyapunov function approach to multi-agent coordination,” in *Proceedings of the 40th IEEE Conference on Decision and Control
(Cat. No. 01CH37228)*, vol. 2.1em plus 0.5em minus 0.4emIEEE, 2001, pp. 1150–1155.

[5]

Z. Artstein, “Stabilization with relaxed controls,” *Nonlinear Analysis: Theory, Methods & Applications*, 1983.

[6]

Z. Gong, M. Zhao, T. Bewley, and S. Herbert, “Constructing control lyapunov-value functions using hamilton-jacobi reachability analysis,” *IEEE Control Systems
Letters*, vol. 7, pp. 925–930, 2022.

[7]

Z. Gong and S. Herbert, “Robust control lyapunov-value functions for nonlinear disturbed systems,” *arXiv preprint arXiv:2403.03455*, 2024.

[8]

L. C. Evans and P. E. Souganidis, “Differential games and representation formulas for solutions of
Hamilton-Jacobi-Isaacs equations,” *Indiana University Mathematics Journal*, 1984.

[9]

S. Bansal, M. Chen, S. Herbert, and C. J. Tomlin, “Hamilton-jacobi reachability: A brief overview and recent advances,” in *Conference on Decision and Control*.1em
plus 0.5em minus 0.4emIEEE, 2017.

[10]

J. F. Fisac, N. F. Lugovoy, V. Rubies-Royo, S. Ghosh, and C. J. Tomlin, “Bridging hamilton-jacobi safety analysis and reinforcement learning,” in *2019 International
Conference on Robotics and Automation (ICRA)*.1em plus 0.5em minus 0.4emIEEE, 2019, pp. 8550–8556.

[11]

S. Bansal and C. J. Tomlin, “Deepreach: A deep learning approach to high-dimensional reachability,” in *2021 IEEE International Conference on Robotics and Automation
(ICRA)*.1em plus 0.5em minus 0.4emIEEE, 2021, pp. 1817–1824.

[12]

M. Chen, S. Herbert, and C. J. Tomlin, “Exact and efficient hamilton-jacobi guaranteed safety analysis via system decomposition,” in *Int. Conf. on Robot. and
Automat.*1em plus 0.5em minus 0.4emIEEE, 2017.

[13]

M. Chen, S. L. Herbert, M. S. Vashishtha, S. Bansal, and C. J. Tomlin, “Decomposition of reachable sets and tubes for a class of nonlinear systems,” *Trans. on Automatic
Control*, vol. 63, no. 11, 2018.

[14]

S. L. Herbert, S. Bansal, S. Ghosh, and C. J. Tomlin, “Reachability-based safety guarantees using efficient initializations,” in *Conf. on Decision and Control*.1em
plus 0.5em minus 0.4emIEEE, 2019.

[15]

J. Liu, P. Zhao, Z. Gan, M. Johnson-Roberson, and R. Vasudevan, “Leveraging the template and anchor framework for safe, online robotic gait design,” in *2020 ICRA*,
2020.

[16]

P. Holmes, S. Kousik, B. Zhang, D. Raz, C. Barbalata, M. Johnson-Roberson, and R. Vasudevan, “Reachable sets for safe, real-time manipulator trajectory design,” 2020.

[17]

C. He, Z. Gong, M. Chen, and S. Herbert, “Efficient and guaranteed hamilton–jacobi reachability via self-contained subsystem decomposition and admissible control sets,”
*Control Systems Letters*, 2023.

[18]

I. M. Mitchell and J. A. Templeton, “A toolbox of Hamilton-Jacobi solvers for analysis of nondeterministic continuous and hybrid systems,” in
*Hybrid Systems: Computation and Control*, 2005.

[19]

M. Chen, S. Herbert, S. Bansal, and C. Tomlin, “Optimal control helper toolbox,” *URL https://github. com/HJReachability/helperOC*.