\section{Introduction}

Neural networks (NNs) offer promising capabilities for data-driven control of complex systems.
However, formally certifying safety properties of systems that are controlled by NNs remains an open challenge in the field.

To this end, recent work developed reachability analysis techniques for Neural Feedback Loops (NFLs), i.e., dynamical systems with NN control policies~\cite{everett2021reachability,hu2020reach,chen2022one,sidrane2021overt}.
For \textit{forward} reachability analysis, these techniques compute the set of states the system could reach in the future, given an initial state set, trained NN, and dynamics model.
Due to the NNs' high dimensionality and nonlinearities, exact analysis is often intractable.
Thus, current methods instead compute guaranteed over-approximations of the reachable sets, based on relaxations of the nonlinearities in the NN~\cite{everett2021reachability,hu2020reach,chen2022one} or dynamics~\cite{sidrane2021overt,everett2021reachability}.
While forward reachable sets can then be used to check whether the closed-loop system satisfies desired properties (e.g., reaching the goal), forward methods can be overly conservative for obstacle avoidance~\cite{Rober22_CDC}.

Thus, this work focuses on \textit{backward} reachability analysis~\cite{everett2021reachability,Rober22_CDC,Rober23_ACC,Rober22_OJCSYS,vincent2021reachable,bak2022closed}.
The goal of backward reachability analysis is to compute backprojection (BP) sets, i.e., the set of states the will lead to a given target/obstacle set under the given NN control policy and system dynamics, as illustrated in~\cref{fig:cartoon}.
The system can be certified as safe if it starts outside the BP sets, but the non-invertibility of NNs presents a major challenge in calculating these sets.
Recent work~\cite{everett2021reachability,Rober22_CDC,Rober23_ACC,Rober22_OJCSYS} proposes to first compute a backreachable (BR) set, i.e., the set of states that lead to the target set for some control input within known control limits.
These methods then relax the NN controller over the BR set, which is a superset of the BP set, and calculate bounds on the BP set via linear programs.

\begin{figure}[t]
    \centering
    \includegraphics[page=6,width=0.8\linewidth,clip,trim=200 100 20 70]{figures/backward_reachability_diagrams.pdf}
    \caption{Backward reachability analysis calculates the set of states ($\mathcal{P}_{-1}, \ldots, \mathcal{P}_{-6}$ in this illustration) that lead into an obstacle under a given control policy. The system is certified as safe if it starts outside those sets, without requiring any trajectory samples (green). This work develops a new algorithm for bounding these sets for NN control policies.}
    \label{fig:cartoon}
    \vspace{-0.2in}
\end{figure}

However, a key challenge with that formulation is the BP set over-approximations can remain loose.
A major cause of this conservativeness is that the BR set is often large, leading to loose NN relaxations, and thus loose BP set approximations.
Moreover, conservativeness compounds when calculated over multiple timesteps due to the wrapping effect~\cite{neumaier1993wrapping}.
Prior work introduced strategies based on set partitioning~\cite{Rober22_CDC,Rober22_OJCSYS,Rober23_ACC} and mixed integer linear programming (MILP)~\cite{Rober22_OJCSYS} to improve tightness, but these methods are fundamentally hindered by initializing with the BR set.
This paper instead proposes a refinement loop, which, for a particular timestep, iteratively uses the previous BP set estimate to relax the NFL and compute a new BP set, which can lead to much tighter BP set estimates.
At each iteration, this new strategy shrinks the domain over which the NFL is relaxed, which leads to a less conservative relaxation and ultimately tighter bounds on the BP sets.
If desired, this idea could be used in conjunction with set partitioning~\cite{Rober22_CDC,Rober22_OJCSYS,Rober23_ACC}.

Another key limitation of prior work is the use of axis-aligned bounding boxes to represent BP set estimates~\cite{Rober22_CDC,Rober22_OJCSYS,Rober23_ACC}.
These hyperrectangles are computed by solving several linear programs (LPs) over states and controls, but hyperrectangles are often a poor approximation of the true BP sets.
Instead, this paper introduces a new approach that enables directly obtaining a closed-form polytope representation for the BP set estimate.
A key reason why this enables tighter BP set bounds is that the facets of the target set (or future BP sets in the multi-timestep algorithm) can be directly used as the objective matrix for the relaxation algorithm, which also automatically addresses the common issue of choosing facet directions in polytope design.
Furthermore, we show that the NN can be relaxed over these polytope BP sets, which leads to much tighter relaxations and BP set estimates compared to prior work, which inflated polytopes to hyperrectangles.

To summarize, the main contribution of this work is DRIP, an algorithm that provides safety guarantees for NFLs, which includes domain refinement and a polytope formulation to give tighter certificates than state-of-the-art methods.
Numerical experiments on the same data-driven control systems as in prior state-of-the-art works~\cite{Rober22_CDC,Rober22_OJCSYS,Rober23_ACC} demonstrate that DRIP can provide $371\times$ tighter bounds while remaining computationally efficient (${\sim0.3}$s) and can certify obstacle avoidance for a robot with a learned NN policy.

\section{Background}
This section defines the dynamics, relaxations, and sets used for backward reachability analysis.

\subsection{NFL Dynamics}
As in~\cite{Rober22_OJCSYS}, we assume linear time-invariant (LTI) dynamics,
\begin{align}
\begin{split}
    \mathbf{x}_{t+1} & = \mathbf{A} \mathbf{x}_{t} + \mathbf{B} \mathbf{u}_t + \mathbf{c} \triangleq f(\mathbf{x}_{t}, \mathbf{u}_t) \label{eqn:lti_dynamics} \, ,
\end{split}
\end{align}
where $\mathbf{x}_{t} \in \mathds{R}^{n_x}$ is the state at discrete timestep $t$, $\mathbf{u}_t \in \mathds{R}^{n_u}$ is the input, $\mathbf{A} \in \mathds{R}^{n_x \times n_x}$ and $\mathbf{B} \in \mathds{R}^{n_x \times n_u}$ are known system matrices, and $\mathbf{c} \in \mathds{R}^{n_x}$ is a known exogenous input.
We assume $\mathbf{x}_t \in \mathcal{X}$ and $\mathbf{u}_t \in \mathcal{U}$, where $\mathcal{X}$ and $\mathcal{U}$ are convex sets defining the operating region of the state space and control limits, respectively.
The closed-loop dynamics are
\begin{align}
\begin{split}
    \mathbf{x}_{t+1} & = \mathbf{A} \mathbf{x}_{t} + \mathbf{B}\pi(\mathbf{x}_t) + \mathbf{c} \triangleq p(\mathbf{x}_t;\pi) \, ,
    \label{eqn:nfl}
\end{split}
\end{align}
where $\pi(\cdot)$ is a state-feedback control policy, discussed next.

\subsection{NN Control Policy \& Relaxation}\label{sec:background:nn_policy_and_relaxation}

The control policy, $\pi: \mathds{R}^{n_x} \to \mathds{R}^{n_u}$, can be an arbitrary computation graph\footnote{The policies used in this work's experiments are feedforward NNs of the form in~\cite{Rober22_CDC,Rober22_OJCSYS,Rober23_ACC} (i.e., alternating linear transformation and ReLU activation).}, provided that it is composed of primitives that can be linearly relaxed~\cite{xu2020automatic,dathathri2020enabling,jaxverify}.
That is, for each non-linearity $\mathbf{z}_i = \sigma_i(\mathbf{z}_{i-1})$, we need to be able to define $\bm{\alpha}_{l,i}, \bm{\beta}_{l,i}, \bm{\alpha}_{u,i}, \bm{\beta}_{u,i}$ such that
\begin{equation}
    \bm{\alpha}_{l, i} \mathbf{z}_{i-1} + \bm{\beta}_{l,i} \leq \mathbf{z}_i \leq \bm{\alpha}_{u, i} \mathbf{z}_{i-1} + \bm{\beta}_{u,i} \, .
    \label{eq:lin_relaxation_zi}
\end{equation}
For each primitive, values of $\bm{\alpha}_{l,i}, \bm{\beta}_{l,i}, \bm{\alpha}_{u,i}, \bm{\beta}_{u,i}$ depend on the algorithm employed and are derived from intermediate bounds over the activation of the network, which can be obtained by applying the bound computation procedure (described next) to subsets of the network.

CROWN~\cite{zhang2018efficient} and LIRPA algorithms in general~\cite{xu2020automatic} enable backward propagation of bounds (from function output to input, not to be confused with backward in time).
For example, consider evaluating a lower bound on $\mathbf{c}\cdot\mathbf{z}_i + b $, using the relaxations defined in~\cref{eq:lin_relaxation_zi} to replace $\mathbf{z}_i$:
\begin{equation}
  \begin{split}
    \min \mathbf{c} \cdot \mathbf{z}_i + b \ \geq & \min\  [\mathbf{c}]^+ \left(\bm{\alpha}_{u, i} \mathbf{z}_{i-1} + \bm{\beta}_{u,i}\right)\\
                                           & \qquad + [\mathbf{c}]^- \left(\bm{\alpha}_{l, i} \mathbf{z}_{i-1} + \bm{\beta}_{l,i}\right) + b\\
                                      \geq & \min \  \left([\mathbf{c}]^+\bm{\alpha}_{u, i}  + [\mathbf{c}]^- \bm{\alpha}_{l, i}\right) \mathbf{z}_{i-1}\\
                                           & \qquad + \left([\mathbf{c}]^+\bm{\beta}_{u, i}  + [\mathbf{c}]^- \bm{\beta}_{l, i} + b\right) \, ,
  \end{split}
\end{equation}
where \mbox{$[\mathbf{c}]^+ = \max(\mathbf{c}, 0)$} and \mbox{$[\mathbf{c}]^- = \min(\mathbf{c}, 0)$}.
Through backsubstitution, we transformed a lower bound defined as a linear function over $\mathbf{z}_i$ into a lower bound defined as a linear function over $\mathbf{z}_{i-1}$.
By repeated application of this procedure to all operations in the network, given an objective matrix $\mathbf{C}$ defined over the output $\mathbf{y}$ of a computation graph, $\mathbf{y} = g(\mathbf{x})$, we can obtain bounds that are affine in the input $\mathbf{x}$:
\begin{equation}
     \mathbf{M} \mathbf{x} + \mathbf{n} \leq \mathbf{C} \mathbf{y}. \label{eqn:affine_bounds_crown}
\end{equation}
If $\mathbf{x}$ is within some known set $\mathcal{X}$, these bounds can be concretized to an actual value by solving \mbox{$\min_{x \in \mathcal{X}} \mathbf{M} \mathbf{x} + \mathbf{n}$}, and closed-form solutions exist for some forms of $\mathcal{X}$.

Prior CROWN-based backward reachability works~\cite{Rober22_CDC,Rober22_OJCSYS,Rober23_ACC} relaxed the control policy $\pi$, whereas this work directly relaxes the closed-loop dynamics, $p$, drawing inspiration from~\cite{chen2022one}, which focused on forward reachability.

\subsection{Backprojection Sets}
Given target set, $\mathcal{X}_T \subseteq \mathcal{X}$, define 1- and $t$-step BP sets as
\begin{align}
    \mathcal{P}_{-1}(\mathcal{X}_T) &= \{ \mathbf{x}\ \lvert\ p(\mathbf{x};\pi) \in \mathcal{X}_T \}
    \label{eqn:one_step_backprojection_sets} \\
    \mathcal{P}_{t}(\mathcal{X}_T) &= \{ \mathbf{x}\ \lvert\ p^{(-t)}(\mathbf{x}; \pi) \in \mathcal{X}_T \} \, ,
    \label{eqn:backprojection_sets}
\end{align}
where $t<0$, and $p^n$ is the $n^{\mathrm{th}}$ iterate of $p$, i.e. $p^{n+1} \triangleq p \circ p^{n}$.
We will typically drop the argument and simply write $\mathcal{P}_{t}$.

The BP sets contain all states that will lead to the target set within some time horizon.
Due to the NN control policy, computing these sets exactly is difficult.
Instead, this work aims to compute BP over-approximations (BPOAs), since these over-approximations are guaranteed to also contain the true (unknown) BP sets and thus provide a conservative description of all states that will lead to the target set.

\subsection{Backreachable Sets \& BPOAs}
Next, we need a domain over which to apply the relaxation tools from~\cref{sec:background:nn_policy_and_relaxation}.
First introduced in~\cite{everett2021reachability}, one idea is to find the set of previous states $\bar{\mathcal{R}}_t$ possibly leading to $\bar{\mathcal{P}}_{t+1}$ regardless of policy, relax the neural network policy $\pi$ over $\bar{\mathcal{R}}_t$ into affine bounds $\pi^L_t(\mathbf{x}_t)$ and $\pi^U_t(\mathbf{x}_t)$, and then derive bounds on $\mathcal{P}_t$ from this linear relaxation.
Specifically, define the set of all possible previous states, i.e., the BR set, as
\begin{equation}
    \bar{\mathcal{R}}_{t}(\bar{\mathcal{P}}_{t+1}) \triangleq \{ \mathbf{x}\ \lvert\ \exists\mathbf{u} \in \mathcal{U} \mathrm{\ s.t.\ }
     f(\mathbf{x},\mathbf{u}) \in \bar{\mathcal{P}}_{t+1} \} \, ,\
    \label{eqn:backreachable_sets}
\end{equation}
where $\bar{\mathcal{P}}_0 \triangleq \mathcal{X}_T$.
Then, \cite{Rober22_CDC,Rober22_OJCSYS,Rober23_ACC} defined a 1-step BPOA as,
\begin{align}
    \bar{\mathcal{P}}_{-1} \!  \triangleq \!  \{ \mathbf{x}\ \lvert\  \exists\mathbf{u} \in [\pi^L(\mathbf{x}), \pi^U(\mathbf{x})] \mathrm{\ s.t.\ } f(\mathbf{x},\mathbf{u}) \in \mathcal{X}_T \}. \label{eqn:bpoa}
\end{align}
$\bar{\mathcal{P}}_{-1}$ is bounded by the hyper-rectangle, ${[\ubar{\mathbf{x}}_{-1}, \bar{\mathbf{x}}_{-1}]}$, which requires solving the following optimization problems for each state $k \in [n_x]$, with the $i$-th standard basis vector, $\mathbf{e}_{i}$:
\begin{equation}
\begin{split}
& \bar{\mathbf{x}}_{-1; k} = \min_{\mathbf{x},\mathbf{u}\ \text{s.t.}\ f(\mathbf{x}, \mathbf{u}) \in \mathcal{X}_T,\ \mathbf{u} \in [\pi^L(\mathbf{x}), \pi^U(\mathbf{x})]} \mathbf{e}_k^\top \mathbf{x} \\ & \ubar{\mathbf{x}}_{-1; k} = \max_{\mathbf{x},\mathbf{u}\ \text{s.t.}\ f(\mathbf{x}, \mathbf{u}) \in \mathcal{X}_T,\ \mathbf{u} \in [\pi^L(\mathbf{x}), \pi^U(\mathbf{x})]} \mathbf{e}_k^\top \mathbf{x}. \label{eqn:lp_obj}
\end{split}
\end{equation}

\section{Approach}
This section introduces our proposed approach, Domain Refinement Iteration with Polytopes (DRIP), with descriptions of the 3 technical contributions in \cref{sec:approach:refinement_loop,sec:approach:polytope_formulation,sec:approach:polytope_input_bounds} and a summary of the algorithm in \cref{sec:approach:algorithm}.

\subsection{Improved BP Set Representation: Polytope Formulation}\label{sec:approach:polytope_formulation}

\begin{figure}
    \centering
    \begin{subfigure}{\linewidth}
        \centering
        \includegraphics[page=1,width=0.65\linewidth,clip,trim=40 70 200 60]{figures/backward_reachability_diagrams.pdf}
        \caption{Given closed-loop dynamics, $f(\mathbf{x}_{t-1}; \pi)$ and target set, $\bar{\mathcal{P}}_t$, want to compute bounds on true, unknown BP set, $\mathcal{P}_{t-1}$.}
        \label{fig:polytope_representation:setup}
    \end{subfigure}\\
    \begin{subfigure}{\linewidth}
        \centering
        \includegraphics[page=2,width=0.65\linewidth,clip,trim=40 10 200 60]{figures/backward_reachability_diagrams.pdf}
        \caption{Calculate BR set, $\mathcal{R}_{t-1}$, using control limits, $f$, and $\pi$}
        \label{fig:polytope_representation:br}
    \end{subfigure}\\
    \begin{subfigure}{\linewidth}
        \centering
        \includegraphics[page=3,width=0.65\linewidth,clip,trim=40 10 200 60]{figures/backward_reachability_diagrams.pdf}
        \caption{Relax C.L. dynamics within $\mathcal{R}_{t-1}$ to get affine bounds}
        \label{fig:polytope_representation:relax}
    \end{subfigure}\\
    \begin{subfigure}{\linewidth}
        \centering
        \includegraphics[page=4,width=0.65\linewidth,clip,trim=40 10 200 60]{figures/backward_reachability_diagrams.pdf}
        \caption{Return $\bar{\mathcal{P}}_{t-1}$ as interval where region between affine bounds overlaps $\bar{\mathcal{P}}_t$}
        \label{fig:polytope_representation:intersection}
    \end{subfigure}
    \caption{Method for calculating $\bar{\mathcal{P}}_{t-1}$ from $\bar{\mathcal{P}}_{t}$ using \cref{alg:one_step_backprojection}. Given a target set (a), compute BR set (b), relax the closed-loop dynamics (c), to get bounds on the BP set (d).}
    \label{fig:polytope_representation}
\end{figure}

\begin{figure}
    \centering
    \begin{subfigure}{0.45\linewidth}
        \centering
        \begin{tikzpicture}
            \draw[color=blue, thick] (1,0) node[anchor=north]{}
              -- (0,1.5) node[anchor=south]{}
              -- (2,1) node[anchor=north]{}
              -- cycle;
            \draw[color=red, thick, fill=orange, fill opacity=0.2] (0, 0) -- (0, 1.5) -- (2, 1.5) -- (2, 0) -- cycle;
        \end{tikzpicture}
        \caption{\label{fig:hypercube_inp}Existing work}
    \end{subfigure}
    \begin{subfigure}{0.45\linewidth}
        \centering
        \begin{tikzpicture}
            \draw[color=blue, thick, fill=orange, fill opacity=0.2] (1,0) node[anchor=east, opacity=1]{$\mathbf{v}_0$}
              -- (0,1.5) node[anchor=east, opacity=1]{$\mathbf{v}_1$}
              -- (2,1) node[anchor=west, opacity=1]{$\mathbf{v}_2$}
              -- cycle;
              \draw (0.8, 0.8) node[cross,red]{};
              \draw[color=red, thick] (1, 0) -- (0.9, 0.4) node[right=-0.1]{$\sigma_0$} -- (0.8, 0.8);
              \draw[color=red, thick] (0, 1.5) -- (0.4, 1.2) node[right=-0.1]{$\sigma_1$} -- (0.8, 0.8);
              \draw[color=red, thick] (2, 1) --(1.4, 0.9) node[below=-0.1]{$\sigma_2$} -- (0.8, 0.8);
        \end{tikzpicture}
        \caption{\label{fig:simplex_inp}This work}
    \end{subfigure}
    \caption{\label{fig:input_relaxation_comparison} Comparison of how a convex polytope $\mathcal{C}$ (blue) is handled by existing methods (a) and our method (b). The orange region denotes the domain used to relax the NN.}
\end{figure}

\begin{figure*}[t]
    \centering
    \includegraphics[page=5, width=0.8\linewidth,clip, trim=0 150 0 0]{figures/backward_reachability_diagrams.pdf}
    \caption{Architecture to enable relaxation of closed-loop dynamics over a polytope description of the current state set. A linear transformation, $\mathbf{V}$ (stacked polytope vertices, $n_v=3$ here), is added to the front of the closed-loop dynamics. This transformation projects the standard $n_v$-simplex onto the state space, leading to a polytope of states over which the closed-loop dynamics is relaxed. During the relaxation, the function including the vertex projection is used to define relaxation parameters, but the backward mode bound propagation stops before the vertex projection to compute $\mathbf{M} \mathbf{x}_t + \mathbf{n} \leq \mathbf{C}\mathbf{x}_{t+1}$.}
    \label{fig:simplex_polytope}
    \vspace{-0.2in}
\end{figure*}

Instead of relaxing $\pi$ and solving LPs to bound $\mathcal{P}_{-1}$, this work instead relaxes the closed-loop dynamics, $p(\ \cdot\ ; \pi)$, which directly results in a closed-form expression for the halfspace representation (H-rep) of the polytope BPOA.
Specifically, for target set $\mathcal{X}_T = \{ \mathbf{x} \ \lvert \mathbf{A}_T\mathbf{x} \leq \mathbf{b}_T \}$, we apply CROWN~\cref{eqn:affine_bounds_crown} on $p$ over the domain $\mathcal{R}_{-1} = \{ \mathbf{x} \ \lvert \mathbf{A}_\text{relax}\mathbf{x} \leq \mathbf{b}_\text{relax} \}$ with objective $\mathbf{C}=\mathbf{A}_T$.
This gives $\mathbf{M}, \mathbf{n}$ such that
\begin{align}
    \mathbf{M} \mathbf{x}_{t-1} + \mathbf{n} &\leq \mathbf{A}_T \mathbf{x}_{t} \quad \forall \mathbf{x}_{t-1} \in \mathcal{R}_{-1}.
\end{align}
Then, the set of $\mathbf{x}_{t}$ that lead to $\mathcal{X}_T$ is bounded by the constraints $\mathbf{M} \mathbf{x}_{t-1} + \mathbf{n} \leq \mathbf{b}_T$ and $\mathbf{x}_{t-1} \in \mathcal{R}_{-1}$.
Thus, we can write the BPOA as
\begin{align}
    \bar{\mathcal{P}}_{-1} &= \left\{ \mathbf{x}\ \lvert \begin{bmatrix} \mathbf{M} \\ \mathbf{A}_\text{relax} \end{bmatrix} \mathbf{x} \leq \begin{bmatrix} \mathbf{b}_T - \mathbf{n} \\ \mathbf{b}_\text{relax} \end{bmatrix} \right\}. \label{eqn:bpoa_as_polytope}
\end{align}
Whereas \cref{eqn:lp_obj} requires solving LPs, \cref{eqn:bpoa_as_polytope} describes the BPOAs in closed form as polytopes.

\cref{fig:polytope_representation} illustrates this idea for a simple 1-state, 1-control system.
\cref{fig:polytope_representation:setup} sets up the formulation: given a target set $\mathcal{X}_T = \bar{\mathcal{P}}_t$ (green) and closed-loop dynamics $p(\cdot; \pi)$, the goal is to compute bounds on $\mathcal{P}_{t-1}$ (all states at the previous timestep, $\mathbf{x}_{t-1}$, that lead into $\mathcal{X}_T$ at time $t$).
\cref{fig:polytope_representation:br} adds the BR set, $\mathcal{R}_{t-1}$, (red) which is calculated by solving 2 LPs (1 minimization, 1 maximization) to describe the set of $\mathbf{x}_{t-1}$ that lead to $\bar{\mathcal{P}}_t$ for some control within the control limits $\mathcal{U}$.
In~\cref{fig:polytope_representation:relax}, the closed-loop dynamics, $p$, are relaxed over $\mathcal{R}_{t-1}$, which provides bounds (blue) on the affine relationship between subsequent states, $\mathbf{x}_{t-1}$ and $\mathbf{x}_{t}$.
Finally, in~\cref{fig:polytope_representation:intersection}, define $\bar{\mathcal{P}}_{t-1}$ based on the intersection points (magenta) of the relaxed dynamics (blue) and target set (green).
The \textit{upper} limit of $\bar{\mathcal{P}}_{t-1}$ is the point where the \textit{lower} bound on the closed-loop dynamics intersects with the \textit{upper} bound of the target set, and vice versa.
This formulation ensures that $\bar{\mathcal{P}}_{t-1}$ contains only states $\mathbf{x}_{t-1}$ for which \textit{some} realization of the closed-loop dynamics (within the relaxation set) will lead to the target set.

More generally, the intersections in~\cref{fig:polytope_representation:intersection} might lie elsewhere or not exist.
We include $\mathcal{R}_{t-1}$ in the constraints that define $\bar{\mathcal{P}}_{t-1}$, both because $\mathcal{R}_{t-1}$ is the domain over which the relaxation is valid, and to ensure $\bar{\mathcal{P}}_{t-1}$ is closed.

\subsection{Improved Relaxation Domain: Refinement Loop}\label{sec:approach:refinement_loop}

While $\mathcal{R}_{-1}$ provides a domain for relaxing the NN, this domain is often excessively conservative and leads to large BPOAs.
To address this issue, this work leverages the insight that $\mathcal{P}_{-1} \subseteq \bar{\mathcal{P}}_{-1} \subseteq \mathcal{R}_{-1}$.
In particular, after computing $\bar{\mathcal{P}}_{-1}$ using the domain $\mathcal{R}_{-1}$, we can compute a new $\bar{\mathcal{P}}_{-1}' \subseteq \bar{\mathcal{P}}_{-1}$ using the domain $\bar{\mathcal{P}}_{-1}$.
Since $\bar{\mathcal{P}}_{-1} \subseteq \mathcal{R}_{-1}$, the relaxation is often tighter in practice, which can lead to a tighter BPOA.
This refinement can be performed iteratively (e.g., for a specified number of steps or until the reduction in BPOA volume between refinement steps reaches some threshold).

To visualize this idea, the magenta lines in~\cref{fig:polytope_representation:intersection} represent the new relaxation domain over which new affine bounds (blue) would be calculated.
These new affine bounds could be much tighter, because the function's range is much smaller over this new domain.

\setlength{\textfloatsep}{4pt
\begin{algorithm}[t]
 \caption{CROWNSimplex}
 \begin{algorithmic}[1]
 \setcounter{ALC@unique}{0}
 \renewcommand{\algorithmicrequire}{\textbf{Input:}}
 \renewcommand{\algorithmicensure}{\textbf{Output:}}
 \REQUIRE closed-loop dynamics, $f$, polytope vertices $\mathbf{V}$, objective matrix $\mathbf{C}$
 \ENSURE $\mathbf{M}$, $\mathbf{n}$, describing affine bounds from $\mathbf{x}_t$ to $\mathbf{x}_{t+1}$
 \STATE $\mathbf{x}_{t} = \mathbf{V} \cdot \mathbf{s}$
 \STATE $\mathbf{x}_{t+1} = f_{\text{simplex}}(\mathbf{s}) = f(\mathbf{V}\cdot \mathbf{s})$ \label{alg:crown_simplex:add_layer}
 \STATE $\bm{\alpha}, \bm{\beta} \leftarrow \text{defineRelaxationParams}(f_\text{simplex}, [\mathbb{0}_{n_v}, \mathbb{1}_{n_v}], 1)$\label{alg:crown_simplex:fwd_bounds}
 \STATE $\mathbf{M}, \mathbf{n} \leftarrow \text{bwdModeProp}(f_\text{simplex}, \mathbf{x}_{t+1} \rightarrow \mathbf{x}_t, \bm{\alpha}, \bm{\beta}, \mathbf{C})$\label{alg:crown_simplex:bwd_bounds}
 \RETURN $\mathbf{M}, \mathbf{n}$
 \end{algorithmic}\label{alg:crown_simplex}
\end{algorithm}

\subsection{Improved Relaxation over BP Set: Polytope Input Bounds}\label{sec:approach:polytope_input_bounds}

Existing algorithms~\cite{everett2021reachability,Rober22_CDC,Rober22_OJCSYS,Rober23_ACC} simplified the bound computation by assuming the domain over which to relax the network was a hyper-rectangle \mbox{$\mathcal{X} = [\mathbf{l}, \mathbf{u}]^d$}, which enables
\begin{equation}
    \min_{x \in \mathcal{X}}\ \mathbf{M} \mathbf{x} + \mathbf{n} = [\mathbf{M}]^+ \mathbf{l} + [\mathbf{M}]^- \mathbf{u} + \mathbf{n} \, .
\end{equation}
Assuming that the input domain is a convex polytope $\mathcal{C}$, the strategy consists in first solving $2d$ LPs, $l_i = \max_{x \in \mathcal{C}} x_i$ and $u_i = \max_{x \in \mathcal{C}} x_i$, to obtain the hyperrectangle.
However, this can contribute to bound conservativeness (\cref{fig:hypercube_inp}).

Instead, we propose to rely on another form of domain for which concretization of linear bounds is efficient: the simplex.
If we define $\mathcal{X} = \Delta_d = \left\{\mathbf{x} \in \mathcal{R}^d  |\  \mathbf{x} \geq 0, \sum_i^d x_i = 1\right\}$, then
\begin{equation}
    \min_{x \in \mathcal{X}}\ \mathbf{M} \mathbf{x} + \mathbf{n} = \mathbf{M}_{\text{argmin[M]}} + \mathbf{n} \, .
\end{equation}
We can represent a convex polytope $\mathcal{C}$, with a simplex and a transformation based on its vertices $\left\{\mathbf{v}_0, \mathbf{v}_1, \dotsc, \mathbf{v}_n\right\}$:
\begin{equation}
    \mathbf{x} \in \mathcal{C} \iff \ \exists \bm{\sigma} \in \Delta_d\ \text{such that } \mathbf{x} = \bm{\sigma} \mathbf{V} \, ,
\end{equation}
where $\mathbf{V} = \texttt{stack}(\mathbf{v}_0, \mathbf{v}_1, \dotsc, \mathbf{v}_n)$.
\cref{fig:simplex_polytope} shows how this is integrated in the bound propagation: the multiplication by the matrix of vertices is simply prepended to the NN controller as an initial linear layer. As seen in~\cref{fig:simplex_inp}, this is an exact representation of polytope $\mathcal{C}$.

\setlength{\textfloatsep}{2pt
\begin{algorithm}[t]
 \caption{Domain Refinement with Polytopes (DRIP)}
 \begin{algorithmic}[1]
 \setcounter{ALC@unique}{0}
 \renewcommand{\algorithmicrequire}{\textbf{Input:}}
 \renewcommand{\algorithmicensure}{\textbf{Output:}}
 \REQUIRE target set $\mathcal{X}_T$, policy $\pi$, dynamics $p$, iterations $n_\text{iters}$
 \ENSURE BP set approximation $\bar{\mathcal{P}}_{-1}(\mathcal{X}_T)$
    \STATE $\mathbf{A}_{\mathcal{X}_T}, \mathbf{b}_{\mathcal{X}_T} \leftarrow \mathcal{X}_T = \{\mathbf{x}\ |\ \mathbf{A}_{\mathcal{X}_T} \mathbf{x} \leq \mathbf{b}_{\mathcal{X}_T}\}$
    \STATE $\bar{\mathcal{P}}_{-1} \leftarrow\! \mathcal{R}_{-1} = \mathrm{backreach}(\mathcal{X}_T, \mathcal{U}, p(\cdot; \pi))$ \label{alg:one_step_backprojection:backreach}
    \FOR{$i$ in $\{1, 2, \ldots, n_\text{iters}\}$} \label{alg:one_step_backprojection:loop}
        \STATE $\mathbf{V} \leftarrow \mathrm{findVertices}(\bar{\mathcal{P}}_{-1})$ \label{alg:one_step_backprojection:vrep}
        \STATE $\mathbf{M}, \mathbf{n} \leftarrow \mathrm{CROWNSimplex}(p(\cdot; \pi), \mathbf{V}, \mathbf{A}_{\mathcal{X}_T})$ \label{alg:one_step_backprojection:crown}
        \STATE $\bar{\mathcal{P}}_{-1} \leftarrow \left\{\mathbf{x}\ |\ \begin{bmatrix} \mathbf{M} \\ \mathbf{A}_{\bar{\mathcal{P}}_{-1}} \end{bmatrix} \mathbf{x} \leq \begin{bmatrix} \mathbf{b}_T - \mathbf{n} \\ \mathbf{b}_{\bar{\mathcal{P}}_{-1}} \end{bmatrix}\right\}$\label{alg:one_step_backprojection:new_hrep2}
    \ENDFOR
 \RETURN $\bar{\mathcal{P}}_{-1}$
 \end{algorithmic}\label{alg:one_step_backprojection}
\end{algorithm}
\begin{algorithm}[t]
 \caption{Multi-Step DRIP}
 \begin{algorithmic}[1]
 \setcounter{ALC@unique}{0}
 \renewcommand{\algorithmicrequire}{\textbf{Input:}}
 \renewcommand{\algorithmicensure}{\textbf{Output:}}
 \REQUIRE target set $\mathcal{X}_T$, policy $\pi$, dynamics $p$, iterations $n_\text{iters}$, time horizon $\tau$
 \ENSURE BP set approximations $\bar{\mathcal{P}}_{-\tau:0}(\mathcal{X}_T)$
    \STATE $\bar{\mathcal{P}}_{0}(\mathcal{X}_T) \leftarrow \mathcal{X}_T$ \label{alg:backprojection:initialize}
    \FOR{$t$ in $\{-1, -2, \ldots, -\tau\}$} \label{alg:backprojection:timestep_for_loop}
        \STATE $\bar{\mathcal{P}}_{t}(\mathcal{X}_T) \leftarrow \text{DRIP}(\bar{\mathcal{P}}_{t+1}(\mathcal{X}_T), \pi, p, n_\text{iters})$ \label{alg:backprojection:one_step_backproj_iter}
    \ENDFOR
 \RETURN $\bar{\mathcal{P}}_{-\tau:0}(\mathcal{X}_T)$ \label{alg:backprojection:return}
 \end{algorithmic}\label{alg:backprojection}
\end{algorithm}

\subsection{Algorithm Overview}\label{sec:approach:algorithm}

To summarize an implementation of the proposed approach, we first describe the method for relaxing $p$ (\cref{alg:crown_simplex}), then show how this is used to calculate the BPOA for a single timestep (\cref{alg:one_step_backprojection}), then describe how to compute BPOAs over multiple timesteps (\cref{alg:backprojection}).

\begin{table*}
    \centering
    \captionof{figure}{Comparison of BP set bounds. Given a target set (red), proposed methods lead to much tighter bounds (dashed lines) on true BP sets (solid lines), by combining multiple iterations (bottom row), polytope target sets (middle column), and simplex bounds during relaxation (right column).}
    \label{fig:double_integrator_sets}
    \setlength{\tabcolsep}{0pt}
    \begin{tabular}{cM{58mm}M{58mm}M{58mm}}

        & LP solver & Polytope target sets (Proposed) & Polytope target sets + simplex bounds (Proposed)  \\

        \rotatebox[origin=c]{90}{1 Iteration} &
        \begin{subfigure}{\linewidth}\centering\includegraphics[width=\linewidth]{figures/reachable_sets/di_JaxBackwardCROWNLPIterativeJitted_None_1.png}\vspace{-0.2in}\caption{Prior Work: BReachLP~\cite{Rober22_CDC}}\label{fig:double_integrator_sets:lp_1}\end{subfigure}
        &
        \begin{subfigure}{\linewidth}\centering\includegraphics[width=\linewidth]{figures/reachable_sets/di_JaxBackwardCROWNRectangleIterativeJitted_None_1.png}\vspace{-0.2in}\caption{Proposed: DRIP-}\label{fig:double_integrator_sets:rect_1}\end{subfigure}
        &
        \begin{subfigure}{\linewidth}\centering\includegraphics[width=\linewidth]{figures/reachable_sets/di_JaxBackwardCROWNPolytopeIterativeJitted_None_1.png}\vspace{-0.2in}\caption{Proposed: DRIP}\label{fig:double_integrator_sets:poly_1}\end{subfigure}
        \\
        \rotatebox[origin=c]{90}{5 Iterations (Proposed)} &
        \begin{subfigure}{\linewidth}\centering\vspace{-0.2in}\includegraphics[width=\linewidth]{figures/reachable_sets/di_JaxBackwardCROWNLPIterativeJitted_None_5.png}\vspace{-0.2in}\caption{Proposed: BReachLP-Iterate}\label{fig:double_integrator_sets:lp_5}\end{subfigure}
        &
        \begin{subfigure}{\linewidth}\centering\vspace{-0.2in}\includegraphics[width=\linewidth]{figures/reachable_sets/di_JaxBackwardCROWNRectangleIterativeJitted_None_5.png}\vspace{-0.2in}\caption{Proposed: DRIP-}\label{fig:double_integrator_sets:rect_5}\end{subfigure}
        &
        \begin{subfigure}{\linewidth}\centering\vspace{-0.2in}
        \begin{tikzpicture}
        \node[inner sep=0pt] (plot) at (0,0)
            {\includegraphics[width=\linewidth]{figures/reachable_sets/di_JaxBackwardCROWNPolytopeIterativeJitted_None_5.png}};
        \node[inner sep=0pt,label={[yshift=-0.2cm]Timestep}] (legend) at (0.3,-0.7)
            {\includegraphics[width=0.8\linewidth]{figures/reachable_sets/di_reachable_set_legend.png}};
        \end{tikzpicture}
        \vspace{-0.3in}\caption{Proposed: DRIP}\label{fig:double_integrator_sets:poly_5}\end{subfigure}
        \\
        \vspace{-0.2in}
    \end{tabular}
    \vspace{-0.3in}
\end{table*}

Given an objective matrix $\mathbf{C}\in \mathds{R}^{h \times n_x}$, \cref{alg:crown_simplex} aims to compute $\mathbf{M}, \mathbf{n}$ such that $\mathbf{C}\mathbf{x}_t \leq \mathbf{M} \mathbf{x}_{t-1} + \mathbf{n}$ $\forall \mathbf{x}_{t-1} \in \texttt{conv\_hull}(\mathbf{V})$.
Add a linear transformation to the front of the closed-loop dynamics to obtain a new function, \mbox{$f_\text{simplex}: \Delta_{n_v} \to \mathcal{X}$} (\cref{alg:crown_simplex:add_layer}).
Then, define the parameters of each relaxation $\bm{\alpha}_{l,i}, \bm{\beta}_{l,i}, \bm{\alpha}_{u,i}, \bm{\beta}_{u,i}$, from \cref{eq:lin_relaxation_zi}, specifying the standard $n_v$-simplex as the bounds on the function's input (\cref{alg:crown_simplex:fwd_bounds}).
Finally, compute backward bounds from the objective, $\mathbf{C}$, to the current state, $\mathbf{x}_{t}$, stopping before reaching the added linear transformation (\cref{alg:crown_simplex:bwd_bounds}).

\cref{alg:one_step_backprojection} describes the proposed method for computing a BPOA, $\mathcal{P}_{-1}$, which uses \cref{alg:crown_simplex}, $\mathcal{X}_T$ (represented as a polytope), and $p$.
First, extract $\mathbf{A}_T, \mathbf{b}_T$ as the H-rep of the target set polytope.
Then, the BR set $\mathcal{R}_{-1}$ is computed by solving $2n_x$ LPs (\cref{alg:one_step_backprojection:backreach}), and the BPOA, $\bar{\mathcal{P}}_{-1}$, is initialized as $\mathcal{R}_{-1}$.
To refine that estimate, loop for $n_\text{iters}$ (\cref{alg:one_step_backprojection:loop}).
At each iteration, compute the vertex representation (V-rep) of the current BPOA (\cref{alg:one_step_backprojection:vrep}) and run \cref{alg:crown_simplex} with $\mathbf{C}=\mathbf{A}_T$ (\cref{alg:one_step_backprojection:crown}).
Then, set the parameters of the H-rep of the refined BPOA using the CROWN relaxation and target set along with the prior iteration's BPOA.
Update $\bar{\mathcal{P}}_{-1}$ and the loop continues.
After $n_\text{iters}$, $\bar{\mathcal{P}}_{-1}$ is returned.

To calculate BPOAs for a time horizon $\tau$ into the past, \cref{alg:backprojection} simply makes $\tau$ calls to \cref{alg:one_step_backprojection}, each time using the next step's BPOA as the target set.

\section{Results}

This section demonstrates DRIP on two simulated control systems (double integrator and mobile robot), implemented with the Jax framework~\cite{jax2018github} and \texttt{jax\_verify}~\cite{jaxverify} library.

\subsection{Ablation Study}

\begin{figure}[t]
    \centering
    \begin{tikzpicture}
        \node[inner sep=0pt] (plot) at (0,0)
            {\includegraphics[width=0.6\linewidth]{figures/runtime_vs_error_2022_11_16__11_43_10_table.png}};
        \node[inner sep=0pt] (legend) at (1.0,-0.3)
            {\includegraphics[width=0.3\linewidth]{figures/runtime_vs_error_2022_11_16__11_43_10_table_legend.png}};
    \end{tikzpicture}
    \caption{Runtime vs. Error. Compared to prior work (BReachLP), the refinement loop lowers the approximation error with more computation time (BReachLP-Iterate). Incorporating  target set facets reduces the error and runtime (DRIP-), and incorporating simplex bounds further improves the error with similar runtime (DRIP). Overall, there is a $371\times$ improvement in error for similar runtime (0.3s).}
    \label{fig:runtime_vs_error}
\end{figure}

Using the double integrator system and learned policy from~\cite{everett2021reachability}, \cref{fig:double_integrator_sets} demonstrates the substantial improvement in reachable set tightness enabled by our proposed method.
For a target set (red) $\mathcal{X}_T = [4.5, 5.0] \times [-0.25, 0.25]$, the true BP sets are shown in solid colors for each timestep, and the BPOAs are shown in the same colors with dashed lines.
For example, \cref{fig:double_integrator_sets:lp_1} implements the prior work~\cite{Rober22_CDC}, in which BPOAs become very conservative after a few timesteps.
By adding the refinement loop (\cref{sec:approach:refinement_loop}), \cref{fig:double_integrator_sets:lp_5} shows substantial improvement with $n_\text{iters} = 5$.

However, increasing $n_\text{iters} > 5$ does not improve the results with the prior LP and hyper-rectangular BPOA formulation.
Instead, the middle column shows the impact of the proposed formulation from \cref{sec:approach:polytope_formulation}, in which the closed-loop dynamics are relaxed to directly provide BPOAs as polytopes.
A key reason behind this improvement is the use of the target set's facets as the objective matrix during the backward pass of the CROWN algorithm.
When increasing $n_\text{iters}$ to 5, we see nearly perfect bounds on the first three BP sets, with tighter yet still somewhat loose bounds on the final two timesteps.

The impact of the polytope domain used in the relaxation (\cref{sec:approach:polytope_input_bounds}) is shown in the rightmost column.
Because the first iteration always uses the hyper-rectangular $\mathcal{R}_{-t}$, we expect to see no difference between \cref{fig:double_integrator_sets:rect_1} and \cref{fig:double_integrator_sets:poly_1}.
However, \cref{fig:double_integrator_sets:poly_5} shows much tighter BPOAs in the last two timesteps when compared to \cref{fig:double_integrator_sets:rect_5}.

Overall, the massive improvement in BPOA tightness between the prior work~\cite{Rober22_CDC} (\cref{fig:double_integrator_sets:lp_1}) and the proposed approach (\cref{fig:double_integrator_sets:poly_5}) would enable a practitioner to certify safety when starting from a larger portion of the state space.

\begin{figure}[t]
    \centering
    \begin{subfigure}{0.45\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/ground_robot_bp_zoomed_in.png}
        \caption{1-step BPOAs (blue) for 15 iterations with obstacle (red)}
        \label{fig:ground_robot:bp}
    \end{subfigure}\hfill
    \begin{subfigure}{0.45\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/ground_robot_bp.png}
        \caption{400 rollouts with target set (red) and BPOAs (blue)}
        \label{fig:ground_robot:trajs}
    \end{subfigure}
    \caption{Ground robot policy certification. With $n_\text{iters}=15$, $\bar{\mathcal{P}}_{-1} \subseteq \mathcal{X}_T$, which certifies that the system will never collide with the obstacle (starting anywhere outside the obstacle).}
    \label{fig:ground_robot}
\end{figure}

\subsection{Runtime vs. Error Tradeoff}

\cref{fig:runtime_vs_error} provides quantitative analysis of the tradeoff between tightness and computational runtime for the various methods.
Final step error is the ratio $\frac{A_\text{BPOA} - A_\text{BP}}{A_\text{BP}}$, where $A_{\mathcal{A}}$ denotes the area of set $\mathcal{A}$, and the runtime is plotted as the mean and shaded standard deviation over 20 trials (with the worst one discarded).
For the 3 variants of the proposed algorithms, by increasing $n_\text{iters}$, we observe an improvement in error at the cost of additional runtime.

The red triangle corresponds to \cref{fig:double_integrator_sets:lp_1}~\cite{Rober22_CDC}.
The pink curve corresponds the leftmost column of~\cref{fig:double_integrator_sets}, which adds a refinement loop on the relaxation domain (\cref{sec:approach:refinement_loop}).
The blue curve corresponds to the middle column of~\cref{fig:double_integrator_sets}, which uses polytope BPOAs (\cref{sec:approach:polytope_formulation}).
The green curve corresponds to the rightmost column of~\cref{fig:double_integrator_sets}, which uses polytope relaxation domains (\cref{sec:approach:polytope_input_bounds}).

A key cause of the computation time reduction between red/pink and green/blue is that the new formulation replaces many of the LPs with closed-form polytope descriptions.
Overall, we observe a $371\times$ improvement in the error.
Moreover, this experiment did not require any set partitioning and was still able to provide very tight BPOAs.

\subsection{Ground Robot Policy Certification}

Next, we use DRIP to certify that a robot will not collide with an obstacle for any initial condition in the state space (outside of the obstacle itself).
We use the 2D ground robot dynamics, policy, and obstacle from~\cite{Rober22_CDC,Rober22_OJCSYS}.
The target set is partitioned into 4 cells (2 per dimension), which we note is $4\times$ fewer cells than in~\cite{Rober22_CDC,Rober22_OJCSYS}.
For each cell, the 1-step BPOA is calculated for $n_\text{iters}$, and the returned BPOA is the convex hull of the union of vertices for each cell's BPOA.

\cref{fig:ground_robot:bp} shows the target set (red) and BPOA for $n_\text{iters}=\{0, 1, \ldots, 15\}$ (blue, darker corresponds to larger $n_\text{iters}$).
When $n_\text{iters}=15$, $\bar{\mathcal{P}}_{-1} \subseteq \mathcal{X}_T$.
By Corollary A.2 of~\cite{Rober22_OJCSYS}, this certifies that the system can only reach the obstacle if it starts from within the obstacle, and there is no need to compute BPOAs for further timesteps.
\cref{fig:ground_robot:trajs} shows the BPOAs with rollouts of the closed-loop dynamics starting from 400 random states (green $\rightarrow$ yellow as $t \uparrow$), although sampling trajectories is not necessary given the certificate.

\section{Conclusion}

This paper proposed a new backward reachability algorithm, DRIP, for formal safety analysis of data-driven control systems.
DRIP advances the state-of-the-art by introducing ideas to shrink the domain over which the closed-loop dynamics are relaxed and leverage polytope representations of sets at both the input and output of the relaxation.
These innovations are shown to provide 2 orders of magnitude improvement in bound tightness over the prior state-of-the-art with similar runtime.
Future work will investigate tighter simplex bounds~\cite{behl2021overcoming} and convergence properties.

\bibliographystyle{IEEEtran}


\section{Introduction}

Intense laser-induced strong-field ionization of the OCS molecule has earned a significant attention in the past decade, where the focus of several experimental works has been on the investigation of the total ionization yields (TIYs) as a function of the orientation angle, $\beta$, of the molecule with respect to the linear polarization of the intense laser pulse ~\cite{HansenJPhysB2011,PhysRevA.89.013405,Johansen_2016,PhysRevA.98.043425}.
The minimum TIY was obtained at $\beta=0^\circ$, that is when the laser polarization is parallel to the molecular axis. At  $\beta=0^\circ$, suppression of ionization is understood since the laser polarization runs through the nodal plane of the highest occupied orbital (HOMO). The maximum TIY from the HOMO of OCS was reported at $\beta=90^\circ$, that is when the laser polarization is perpendicular to the molecular axis. This measured orientation dependence of the TIYs presented a challenge to ionization models. Neither the molecular Ammosov-Delone-Karinov model~\cite{PhysRevA.66.033402},  the strong-field approximation~\cite{sfa1,sfa2,sfa3}, the weak-field asymptotic theory~\cite{Tolstikhin2011,PhysRevA.87.013406}, the Stark-corrected molecular tunneling theory~\cite{PhysRevA.82.053404,Holmegaard2010}, or
the Stark-corrected molecular tunneling theory~\cite{PhysRevA.82.053404,Holmegaard2010} could reproduce these experimental results for the OCS molecule.
The behaviour of the measured $\beta$-resolved TIY could, however, be predicted by time-dependent density-functional theory accompanying the experiment in reference \cite{PhysRevA.98.043425}. Recently, we showed that the TIY for OCS as a function of $\beta$ can also be captured by an extension of the single-active-electron (SAE) time-dependent Schr\"odinger equation (TDSE) approach~\cite{PhysRevA.102.063111}. The methodology consists of supplementing the SAE TDSE with an account for the multielectron polarization (MEP) effect as described by the induced-dipole potential, reflecting the laser-induced polarization of the electrons in the cation. With this approach, we have also explained strong-field ionization results in CO~\cite{PhysRevA.101.013433}, as well as in O$_2$, CO$_2$ and  CS$_2$~\cite{PhysRevA.106.013117}; see also the related works by another groups based on very similar approaches~\cite{PhysRevA.95.023407,0953-4075-51-10-105601,PhysRevA.105.023106}.

So far an investigation of the ability to predict observables that are highly influenced by excited states in strong-field ionization has not been performed in OCS with the effective SAE TDSE approach including MEP effects. Insights into such a possibility is nevertheless highly desirable,  as this could allow one to apply the methodology with confidence for a wider range of observables. Now, experiments published in reference~\cite{Yu_2017} lend us such a possibility for studying excited-state dynamics. In that work, ionization of OCS was considered at different intensities for driving wavelengths of 800 nm and 400 nm. Signature peaks in the intensity-dependent above-threshold ionization (ATI) spectra that did not shift as expected from a direct-ionization model were identified and associated with ionization from excited  states in OCS, in a resonance-enhanced multiphoton ionization (REMPI) process. In this work, we will use the characteristics of the  experimental findings to explore the performance of the SAE TDSE approach including MEP for dynamics that involve excited states. This investigation will, as in the experiment, observe the behaviour of the ATI spectra with varying peak laser intensity. This intensity scan allows us to see effects associated with excited states shifting in and out of resonance in an intensity-dependent manner.
The comparison between our results and the experiment indicates that our SAE potential for OCS including MEP captures the essential physics needed for a description of strong-field ionization involving REMPI mechanisms.

The theoretical and computational models are presented in section \ref{compdet}, followed by results and discussion in section \ref{res} and conclusions in section \ref{conc}. Atomic units are used throughout unless otherwise stated.

\section{Theoretical and Computational Models}
\label{compdet}
In this section, we give a brief overview of the procedures for generating the SAE potential and the wavefunction for the HOMO of OCS, full details are available elsewhere~\cite{PhysRevA.102.063111}. We also provide a brief overview of our TDSE method including the MEP correction; see references \cite{PhysRevA.101.013433,PhysRevA.102.063111,PhysRevA.106.013117} for further description.

The SAE potential for OCS was obtained from quantum chemistry calculations using density-functional theory (DFT), following the procedure described in references \cite{PhysRevA.81.033416,PhysRevA.102.063111}. This potential was expanded in partial waves as $V(\vec{r})=\sum_{l,m=0}^{l_{max}}V_{l0}(r)Y_{l0}(\theta,\phi)$ where $l_{max}=20$. All partial waves in the expansion of $V(\vec{r})$  have $m=0$ since the molecule is linear and the potential is invariant under rotations around the molecular axis. Based on our SAE potential, the HOMO of OCS is a $\pi$ orbital with energy $-0.35$~a.u., which is slightly above a reference orbital energy of $-0.42$ a.u. as obtained from quantum chemistry calculations~\cite{gamess}. The wavefunction corresponding to the HOMO of OCS was obtained using the split-operator spectral method ~\cite{PhysRevA.38.6000}, see reference \cite{PhysRevA.102.063111} for more details. The obtained wavefunction, after being checked for stability on the grid, is used as the initial state for the solution of the TDSE.

The MEP term is approximated as $- \vec{\mu}_{ind}\cdot \vec{r} /r^3 \approx -\alpha_{||} \vec{E}(t)\cdot\vec{r}/r^3$, where $\alpha_{||}$ is the static polarizability of the cation parallel to the laser polarization axis. Reference ~\cite{PhysRevA.106.013117} suggested that there exist a non-vanishing contribution to the induced dipole from the polarizability induced perpendicular to the laser polarization $\alpha_{\perp}$. For the OCS molecule, however, at the orientation angle considered here ($\beta=0^\circ$, corresponding to alignment of the molecular axis with the linear polarization axis), the contribution from $\alpha_{\perp} = (\alpha_{zz}-\alpha_{xx})\sin(\beta)\cos(\beta)$ vanishes. At other orientation angles, our calculations of the induced dipole moment show that  the contribution from $\alpha_{\perp}$ is small and can be neglected. The polarizability $\alpha_{||}$ was computed based on the polarizability components for the OCS$^+$ cation which were reported in our previous work~\cite{PhysRevA.102.063111}. At short range,  a cutoff radius is selected at a radial distance $r_{c}=\alpha_{||}^{1/3}$ such that the MEP cancels the external field at $r\le r_c$~\cite{0953-4075-51-10-105601,PhysRevLett.95.073001,doi:10.1080/09500340601043413}. The reasoning behind this cutoff is that electrons remaining in the cation polarize and thereby set up a field that counteracts the externally applied field. The expression of $r_c$ in terms of $\alpha_{||}$ is determined by setting the sum of the time-dependent interaction potentials equal to zero, $\vec{E}(t)\cdot\vec{r}(1-\alpha_{||}/r_c^3) =0$.

The TDSE is solved for the HOMO of OCS in an effective one-electron potential describing the interaction with the nucleus, the remaining electrons, and the external field, including the induced-dipole MEP term. The TDSE is propagated in the length gauge  with a combined split-operator~\cite{PhysRevA.38.6000}, Crank-Nicolson method; see reference~\cite{Kjeldsen2007a}. The electric dipole field of the laser pulse, $\vec E(t)$, linearly-polarized along the laboratory-frame $z$-axis, is defined in terms of the vector potential $\vec A(t)$ by
\begin{equation}
    \vec{E}(t) = -\partial_t \vec{A}(t) = -\partial_t \left(\frac{E_0}{\omega}\sin^2(\pi t/\tau)\cos(\omega t+\varphi) \right)\hat{z},
    \label{E_field}
\end{equation}
where $E_0$ is the field amplitude, $\omega$ the angular frequency, and $\varphi$ the carrier-envelope phase (CEP) for a laser pulse with duration $\tau$. We consider pulses containing 5 and 8 cycles at angular frequencies of $\omega=0.057$ and 0.114~a.u., corresponding to wavelengths of 800 and 400~nm, respectively. The CEP is kept fixed  ($\varphi=-\pi/2$). Convergence of the results is ensured by using a time-step of $\delta t =0.001$~a.u. in the time propagation and a radial grid containing up to 8192 points and extending to 320 a.u. The size of the angular basis set is limited by setting $l_{max}=60$ in the partial wave expansion of the wavefunction.  The ATI spectra were produced by projecting the wavefunction at the end of the laser pulse on scattering states of the Coulomb potential in the asymptotic region~\cite{Madsen2007}.

\section{Results and discussion}
\label{res}

\begin{figure*}
\includegraphics[width=1.0\textwidth]{OCS_2DATI_800nm_Exp_theo.PNG
\caption{\label{fig1} (a) Experimental and (b) Theoretical 2D ATI spectra of OCS at intensities in the range from 28.4 to 94.3 TW/cm$^2$ for 800-nm laser fields. Experimental results are reproduced from reference~\cite{Yu_2017}. The theoretical results were obtained for OCS aligned along the linear polarization direction, i.e., at orientation angle $\beta=0^\circ$. The solid, black lines in (a) and (b) show the kinetic energy positions predicted by equation (\ref{KE}), the integer values $n$ entering that formula are given in (a) and (b). The dashed, black, horizontal lines indicate the position of intensity-independent ATI features (peak1, peak2, and peak3 in (a), and peak1 and peak2 in (b)) associated with resonant ionization.}
\end{figure*}

\subsection{Ionization by 800-nm laser}
In figure~\ref{fig1}~(a) we show intensity-dependent ATI spectra for randomly-oriented OCS from recent experiments~\cite{Yu_2017}. To guide the eye, the spectra are complemented with black, full lines that show how the energy of the ATI peaks would change with laser intensity, assuming the electrons were ionized directly from the HOMO without participation of resonances. The kinetic energy of the photoelectron for each line is calculated using the equation
\begin{equation}
  E_{kin}= n\omega - I_p - U_p,
  \label{KE}
\end{equation}where $n$ is the number of absorbed photons, $I_p$ the ionization potential, and the pondermotive shift $U_p = E_0^2/(4\omega^2)$. The first ATI peak, observed at laser intensities in the range 30$-$45~TW/cm$^2$ corresponds to the OCS molecule absorbing $n=9$ photons (as indicated by the label $n=9$). For the OCS molecule in the intensity range 30$-$60~TW/cm$^2$, the $n=9, 10, 11$, and 12 ATI peaks  shift to lower energy upon increasing the laser intensity, suggesting that at these intensities the HOMO electron is directly ionized into the continuum, without contributions from excited states. At laser intensities exceeding 60~TW/cm$^2$, however, new features appear in the spectra which have been associated with ionization via excited states. These are labeled peak1, peak2, and peak3 in figure~\ref{fig1}~(a) and centered at fixed kinetic energies of 0.7, 1.2 and 2.2 eV. The feature peak1 is clearly visible in the intensity range 60$-$80~TW/cm$^2$. The feature peak3 is a progression of peak1 caused by absorption of an additional photon. These features (peak1 and peak3) have been attributed in reference~\cite{Yu_2017} to (6+4)- and (6+5)-photon resonant ionization through the $^1\Delta$ excited state of OCS: at the considered intensities the HOMO will be excited into the $^1\Delta$ excited state upon absorption of 6 photons and ionized into the continuum upon absorption of additional 4 photons. Regarding the feature labeled peak2, it is fixed at a kinetic energy of 1.2~eV and, compared with peak1 and peak3, appears in the ATI spectra at higher intensities in the range 80$-$95~TW/cm$^2$. This feature (peak2) has been associated with (8+3)-photon ionization through the $^1\Pi$ excited state of OCS~\cite{Yu_2017}. The energy positions of these resonance-mediated features (peak1, peak2, and peak3) are intensity-independent. In reference~\cite{Yu_2017} this characteristic is explained by assuming that the excited resonance states are AC-Stark shifted by the ponderomotive potential, and therefore their energy follow the ponderomotive-potential-shifted ionization threshold, $I_p+U_p$. We note that while the ionization threshold and highly-excited states converging to this threshold are AC-Stark shifted by $U_p$, it is not to be expected that excited states located  as far as $\sim 3, 4$ or 5 photons below threshold follow the simple $U_p$ Stark shift.  Such reasoning could therefore lead to a questioning of the accuracy of the assignment of the relatively low-lying excited states to the REMPI pathways discussed above.

The intensity-dependent ATI spectra of OCS, obtained from our TDSE calculations, are shown in figure~\ref{fig1}~(b). The numerical cost of the extensive scan over laser peak intensities prevents a consideration of laser focal volume and orientation effects. Accordingly, we only show results for OCS with the molecular axis aligned with the linear polarization direction, i.e.,  at  fixed orientation angle $\beta=0^\circ$, and fixed peak intensities. In spite of these restrictions, we find qualitative similarities between our theoretical results (cf. figure~\ref{fig1}~(b)) and the experimental data: at laser intensities in the range 30-60~TW/cm$^2$, the ATI peak positions and their shift with laser intensity suggest that the HOMO electron is directly ionized without participation of high-lying resonances, in agreement with the experimental findings of reference~\cite{Yu_2017}. Moreover,  features associated with resonant ionization via $U_p$ Stark shifted high-lying states, in particular peak1 and peak2, are clearly seen in our theoretical calculations in the same intensity range, 60$-$95~TW/cm$^2$, as reported experimentally.
The feature peak3 is not well resolved in our TDSE calculations, therefore it will not be discussed here.

Now, a comparison between our theoretical results and the experiment reveal some differences. First of all our SAE model for OCS underestimates the energy of the HOMO: our prediction of the HOMO energy is -0.35~a.u. compared with a literature value of -0.42~a.u., see section~\ref{compdet}. This means that we need less photons to ionize than experimentally, and the observed shift in energy explains the difference in the photon-absorption numbers assigned to the ATI peaks in figure 1 following from the direct ionization model of equation (\ref{KE}). Our results show that it takes a minimum of 8 photons to direct ionization from the HOMO orbital. So, the theoretical direct ATI channels in figure 1(b) are labeled by $n=8, 9, 10, 11, 12$.
Moreover, turning to the resonance features (peak1 and peak2) in figure~\ref{fig1}, the energies of these peaks from TDSE calculations seems to be shifted to lower values in comparison with the experiment: peak1 (peak2) is observed experimentally at electron energy of $\sim 0.7$~eV ($\sim 1.2$ eV), whereas TDSE calculations predict this peak at $\sim 0.5$~eV ($\sim 0.9$ eV).  A possible explanation for this discrepancy is errors in the energies of excited states of OCS predicted by the SAE model.
This is both  plausible and consistent with our preceding discussion of the SAE performance for OCS, namely that our estimate of the ionization potential for the HOMO orbital being somewhat lower than the literature value. We will come back to the experimental and theoretical excited-state energy positions in section 3.3 below.

\subsection{Ionization by 400-nm laser}

\begin{figure*}
\includegraphics[width=1.0\textwidth]{OCS_2DATI_400nm_Exp_theo.PNG
\caption{\label{fig3} (a) Experimental and (b) Theoretical 2D ATI spectra of the oriented OCS molecule at orientation angle $\beta=0^\circ$ and intensities in the range from 3.5 to 35 TW/cm$^2$ for 400-nm laser fields. The experimental results reproduced from reference~\cite{Yu_2017}. The solid, black lines in (a) and (b) indicate the positions of the ATI peak energies predicted by equation (\ref{KE})
based on direct ionization model without contribution from excited states; the integer values $n$ entering equation~(\ref{KE}) are given in (a) and (b). The dashed lines indicate the positions of intensity-independent ATI features (peak4) associated with resonant ionization via  excited states. }
\end{figure*}

In figure \ref{fig3}(a), we show the experimental 2D ATI spectra of OCS at 400-nm wavelength as obtained from reference \cite{Yu_2017}. From the figure, one can see 4 ATI peaks resulting from absorption of, respectively, $n=4, 5, 6$, and 7 photons. In reference~\cite{Yu_2017}, analysis of contributions from excited states was only considered for the lowest ($n=4$) ATI peak, it was not considered for the higher order ATI peaks ($n=5, 6$, and 7) simply because the ionization yield of these ATI peaks is diminished by going to lower laser intensities. Here, we will follow the same approach and will only focus on the imprint of resonance ionization in the first ATI peak. In figure~\ref{fig3}~(a), the electron kinetic energies calculated based on the direct ionization model of eqaution (\ref{KE}) are shown as solid, black lines. According to the direct ionization model, the $n=4$ ATI peak is expected to shift slightly to lower kinetic energy with increasing intensity. However, at the considered intensity range (3.5$-$34~TW/cm$^2$), the ATI peak seems to be fixed at an intensity-independent energy of about 0.9 eV. This ATI feature
was labeled peak4 in figure~\ref{fig3}~(a) and its energy position is represented by the dashed, black  line in the figure. This feature was been attributed in reference~\cite{Yu_2017} to (3+1)-photon resonant ionization via the 4p$\sigma$ $\Pi$ Rydberg state.

In figure~\ref{fig3}(b), we show the 2D ATI spectra for OCS as obtained from the present TDSE calculations at a fixed orientation angle of $\beta=0^\circ$. The OCS molecule was probed by 400-nm laser light at intensities in the range 3.5 – 35 TW/cm$^2$. With increasing intensity, our calculations do not show the higher-order ATI peaks ($n=5$, $n=6$, $n=7$) as clearly as  observed experimentally. This is most likely because our calculations were conducted at fixed molecular orientation.
The black solid lines indicate the intensity-dependent positions of the ATI peaks in the case of direct ionization, without contribution from excited states, i.e., the positions predicted from equation (\ref{KE}). Our TDSE calculations of the first ATI peak ($n=4$ peak) show that this peak corresponding to peak4 in figure 2(a) is fixed at an energy of 2.1 eV in particular at laser intensities up to 20~TW/cm$^2$.
We note that in the experimental measurements, the resonance peak (peak4) was observed at a kinetic energy of 0.9 eV. As discussed in the previous section, a possible explanation for this discrepancy is errors in the energies of excited states of OCS as predicted by the SAE model. This aspect will be discussed further in the following.

\begin{figure*}
\includegraphics[width=0.5\textwidth]{ocs_energy_spec_400vs800.png
\caption{\label{fig4} Population in $\Pi$-symmetry excited states of the OCS upon probing the molecule by five-cycle pulse with laser peak intensities of 9~TW/cm$^2$ (at 400 nm; dashed line) and 88~TW/cm$^2$ (at 800 nm wavelength; solid line).}
\end{figure*}

\subsection{Population in excited states and intensity-independent peak features in the ATI spectra }

As can be seen from the experimental results at 800-nm wavelength in figure~\ref{fig1}~(a), the resonance features peak1 and peak2,
are observed at different intensity ranges. In particular, peak1
is observed in the intensity range 60$-$80~TW/cm$^2$ and peak2
in the range 80$-$95~TW/cm$^2$. Regarding the experimental results at 400-nm wavelength in figure~\ref{fig3}~(a), the resonance feature (peak4) is observed at intensities in the range 5$-$15~TW/cm$^2$, see dashed line in figure~\ref{fig3}~(a).

In the previous section, we established a qualitative agreement between the theoretical intensity-dependent ATI spectra obtained by our TDSE calculations and experimental results. In particular, the ATI features associated with resonance ionization were clearly identified in our calculations.
In the following, we conduct analysis of excited state population based on our TDSE calculations at laser peak intensities where peak2 (at 800~nm) and peak4 (at 400~nm) are observed. To help identify which excited states contribute to ionization at the $\beta=0^\circ$ orientation,  we show in figure~\ref{fig4} population of the excited states of OCS after probing the molecule by five-cycle (800~nm) and eight-cycle (400~nm) pulses with a peak laser intensity of 9~TW/cm$^2$ at 400-nm and 88~TW/cm$^2$ at 800-nm wavelength. At $\beta=0^\circ$, the laser only couples excited states with the same $\pi$ symmetry as the HOMO of OCS. Notice that we use the lower-case letter $\pi$ to describe the symmetry of the active orbital and we will use the wording 'orbital' and not 'state', since the latter is reserved for the many-electron complex. As can be seen from figure \ref{fig4}, when our  OCS  model is probed by the 800-nm pulse, three excited orbitals of $\pi$ symmetry are populated, with the most prominent at an energy of -0.022~a.u. Notice that ionization from this excited orbital by absorption of a single photon results in the first ATI peak at energy of $E_{kin}=-0.022+ \omega~\approx 1.0$~eV. The expected $E_{kin}$ value is in agreement with our theoretical calculations of the peak2 position in figure~\ref{fig1}(b). Notice that the excited orbitals we identify at 800-nm will ionize by a single-photon absorption. However, in reference~\cite{Yu_2017}, the resonance features was assigned to an excited state that would require absorption of additional 3 photons to ionize. This means that this excited-state resonance is located relative far away from the ionization threshold, and it is not obvious that it will satisfy the $U_p$-energy shifting assumed in the data analysis.

In the case when our OCS model is probed by 400-nm pulses with a peak intensity of 9~TW/cm$^2$, the excited orbital at energy of -0.029~a.u. shifts into 3-photon resonance with the HOMO. Ionization by single-photon absorption from this excited orbital results in the first ATI peak at energy of $E_{kin}\approx 2.2$~eV. The expected $E_{kin}$ value from this analysis is in agreement with our theoretical predictions of the $n=4$ position in figure~\ref{fig3}(b). Notice that the results of these analyses are also consistent with the assumption that the AC-Stark shift of the resonantly highly excited orbitals can be well-described by $U_p$. Accordingly, the kinetic energies of electrons ionized from these orbitals remain constant upon increasing the laser intensity as was observed for OCS experimentally~\cite{Yu_2017} and as is now confirmed by the present TDSE results in figures~\ref{fig1} and \ref{fig3}.

\section{Summary and Conclusions}
\label{conc}
Recent experimental work~\cite{Yu_2017} addressed the contribution from excited states to strong-field ionization of randomly oriented OCS at wavelengths of 800 and 400 nm. They associated ATI peaks with ionization from certain excited states. We conducted theoretical analysis of the contribution of excited orbitals to strong-field ionization from the HOMO of the OCS molecule based on TDSE calculations including MEP effects at fixed orientation angle $\beta=0^\circ$. Our theoretical results at both laser wavelengths are in qualitative agreement with the experimental findings. We notice that this agreement is obtained although the present calculations were conducted at fixed molecular orientation of $\beta=0^\circ$ (only $\pi$-symmetry excited orbitals are considered at this orientation) and without accounting for focal-volume averaging. The ability of the present approach to capture essential aspects of excited-state physics under a strong external field, will pave the way for further theoretical investigations of nonlinear processes involving excited molecular states. To mention a few recent applications, this could, for example, include investigations of the role of highly-excited molecular states in analysis of multiphoton-induced free-induction decay in St\"uckelberg spectroscopy ~\cite{PhysRevResearch.4.023135} or in high-order harmonic generation by counterrotating circularly polarized fields~\cite{PhysRevLett.129.173202}.

\section{Acknowledgement}
The time propagation was performed on the cluster of the Centre of Scientific Computing Aarhus (CSCAA). All wavepacket analysis were performed using the Phoenix High Performance Computing facility at the American University of the Middle East (AUM), Kuwait.

\section{Bibliography}
\providecommand{\noopsort}[1]{}\providecommand{\singleletter}[1]{#1}%