\section{Introduction}

The recently discovered uranium-based unconventional superconductor  UTe$_{\rm 2}$~has attracted a lot of attention as a promising candidate for spin-triplet pairing and topological superconductivity \cite{ran2019}. Yet, despite intense experimental efforts (see Ref.~\cite{aoki2022}), the nature of the superconducting order parameter (OP) of UTe$_{\rm 2}$~remains elusive.
Among its numerous peculiarities, there is growing evidence that UTe$_{\rm 2}$~hosts different superconducting phases as a function of applied hydrostatic pressure  \cite{braithwaite2019,thomas2020,aoki2020,ran2020,ran2021} and magnetic field \cite{ran2019_2,knebel2019,aoki2020,ran2020,ran2021}. However, the number of superconducting transitions at ambient pressure and zero magnetic field seems to depend on sample details and is thus still debated \cite{aoki2022}.

At zero field and ambient pressure, several studies suggested that UTe$_{\rm 2}$~is a chiral superconductor based on the observation of chiral surface states \cite{jiao2020,bae2021}, a gap structure with point nodes \cite{metz2019}, and a broken time reversal symmetry (TRS) in the superconducting state \cite{hayes2021,wei2022}.
 Due to the presence of two distinct thermodynamic superconducting anomalies in some samples \cite{hayes2021}, UTe$_{\rm 2}$~is proposed to possess two superconducting OPs. Additionally, the trainability of the polar Kerr signal with magnetic field along the crystallographic $c$ axis suggests that the product of these two superconducting OPs transforms as the $B_{1g}$ irreducible representation (irreps) of the orthorhombic point group D$_{2h}$ \cite{hayes2021,wei2022}. Importantly, because this point group has no multi-dimensional irreps, a description in terms of a two-component gap necessarily requires near-degenerate superconducting instabilities.

While some samples display two features in the specific heat across the superconducting transition, UTe$_{\rm 2}$~crystals that display an optimal superconducting transition temperature $T_{\rm c} = 2\, \rm K$ and large residual resistivity ratios host a single thermodynamic superconducting transition, as manifested by a single jump in the specific heat \cite{cairns2020,rosa2021}.
On the one hand, this could be an indication that the two OPs condense at very close temperatures that cannot be resolved in the specific heat. On the other hand, there has been no report of broken TRS in the superconducting phase of these samples showing a single transition. As a result, the nature of the superconducting OP is still unclear.

 In samples with a single superconducting transition (i.e. one jump in the specific heat), pressure splits it into two thermodynamic transitions that have opposite pressure dependence above $0.3\, \rm GPa$ \cite{braithwaite2019}. For samples with two peaks in the specific heat at ambient pressure, four peaks are observed above $0.3\, \rm GPa$ \cite{thomas2020}.
Moreover, upon further increasing the pressure, an antiferromagnetic phase emerges \cite{braithwaite2019,thomas2020}, which contrasts with the ferromagnetic fluctuations expected at ambient pressure \cite{sundar2019,ran2019}. The connection between this pressure-induced change in magnetic fluctuations and the pressure-induced splitting of the superconducting transitions remains unclear \cite{ishizuka2021}.

 \begin{figure*}[t]
\centering
\includegraphics[width=1\textwidth]{Fig1.pdf}
\caption{Temperature dependence of the heat capacity $C/T$ (top panels) and electrical resistivity (bottom panels) along the indicated crystal directions ($\rho_{100}$, $\rho_{110}$ and $\rho_{001}$) at the indicated uniaxial stress values and orientations [(a), (d): $\sigma_{100}$; (b), (e): $\sigma_{110}$; (c), (f): $\sigma_{001}$]. ac calorimetry data were recorded with an excitation frequency $f\approx 20\, \rm Hz$. Colored dashed lines mark, for each stress value, the average temperature of the sharp rise of the jump in $C/T$ at $T_{{\rm c},ac}$ on the top panels and the temperature $T_{{\rm c},\rho}$ below which the resistivity of the sample is zero on the bottom panels.  Colored arrows show the trends of $T_{{\rm c},ac}(\sigma)$ and $T_{{\rm c},\rho}(\sigma)$ upon increasing the compressive stress. Black arrows on the top panels show the onset [(a), (b)] and end [(c)] of the heat capacity jump of the lowest stress curve (black lines). The sketches on the top panels show the direction of the applied stress. Resistivity values $\rho_{110}$ and $\rho_{001}$ are scaled ($\times 2$ and $\times 1/12.5$ respectively) for clarity.  The difference between the resistivity curves at the two lowest $\sigma_{110}$ values is due to a crack in the sample, not to reversible piezoresistivity.}
\label{Fig1}
\end{figure*}

Uniaxial stress has proved to be a powerful tool to study multi-component superconductors. This technique has been used extensively to study the phase diagram and the OP of $\rm Sr_2RuO_4$  \cite{li2020,grinenko2021}. While this material was initially thought to display spin-triplet pairing \cite{ishida1998}, recent NMR data performed under strain demonstrated it to be actually a singlet superconductor \cite{pustogow2019}.

Here we investigate the nature of the superconducting OP of UTe$_{\rm 2}$~by measuring the low temperature ac heat capacity and the electrical resistivity of single crystals under uniaxial stress $\sigma_{100}$, $\sigma_{110}$ and $\sigma_{001}$  respectively applied along the [100], [110] and [001] crystallographic directions (details in the Methods section of Supplemental Materials). Our unstressed samples display a single transition according to the specific heat.
 We report that the $T_{\rm c}$ value extracted from calorimetry decreases with compressive $\sigma_{100}$ and $\sigma_{110}$ but increases with compressive $\sigma_{001}$.
   Most importantly, we show that a symmetry-breaking in-plane shear stress $\sigma_{xy}$ does not induce any observable splitting of the superconducting transition. We discuss different scenarios to explain these results, including one in which the superconducting OP of UTe$_{\rm 2}$~is single-component and another one in which the two components belong to different symmetry channels from those proposed previously \cite{hayes2021,wei2022}.
  Aside from its effect on $T_{\rm c}$, $c$-axis stress induces a significant piezoresistivity, presumably caused by
 the reduction of the energy scale corresponding to the feature observed at $T^\star\approx 15\, \rm K$, which has been attributed to the of onset short-range magnetic correlations or anisotropic Kondo coherence \cite{willa2021,kang2021,eo2021}.

\begin{figure*}[t]
\centering
\includegraphics[width=1\textwidth]{Fig2.pdf}
\caption{\textbf{(a)} Uniaxial stress dependence of the superconducting transition temperatures  $T_{{\rm c},\rho}$ (circles) and $T_{{\rm c},ac}$ (crosses) obtained for applied  $\sigma_{100}$ (red),  $\sigma_{110}$ (blue) and $\sigma_{001}$ (green) directions. Full symbols correspond to data from Fig.~\ref{Fig1} and open symbols from panel (d) at higher $\sigma_{001}$ values. Colored dashed lines are linear fits to $T_{{\rm c},ac}(\sigma)$, whose slopes $\frac{{\rm d}T_{{\rm c},ac}}{{\rm d}\sigma}$ are displayed on the lower left corner of the figure.
\textbf{(b) \& (c)} Temperature dependences of the heat capacity $C/T$ under compressive stress $\sigma_{100}=0.050\,\rm GPa$ [panel (b)] and $\sigma_{001}=0.200\,\rm GPa$ [panel (c)], at the indicated oscillating current heater frequencies. 5-8 times higher powers were used to improve resolution and raw data (light colors) are smoothed (dark colors) to better show the different features. For clarity, amplitudes have been renormalized and offset to show a comparable jump amplitude. Dashed lines indicate the temperature of the main heat capacity jump, that shows no frequency dependence. Colored arrows show the evolution of the onset temperature of the foot above $T_{\rm c}$~[panel (b)] and shoulder below $T_{\rm c}$~[panel (c)].
 \textbf{(d) \& (e)} Temperature dependence of the $c$-axis resistivity $\rho_{001}$ at different $c$-axis stress levels $\sigma_{001}$.  (d) Resistivity below $5\, \rm K$ and plotted versus $T^2$ together with $\rho_{001} = \rho_0 +AT^2$ fits to the data between $2.6\, \rm K$ and $4.2\, \rm K$ (dashed lines). The inset shows the $\sigma_{001}$ dependance of the $\rho_0$ and $A$ parameters extracted from the fits. (e) Stress evolution of the position of the peak at $T^\star \approx 15\, \rm K$ (colored arrows). The inset shows the data at $\sigma_{001}=-0.030\, \rm GPa$ below room temperature. }
\label{Fig2}
\end{figure*}

The ac calorimetry data displayed in Figs.~\ref{Fig1}~(a)~-~(c) show the temperature dependence of the heat capacity plotted as $C/T$ of three samples at the indicated $\sigma_{100}$, $\sigma_{110}$ and $\sigma_{001}$ uniaxial stress values. Note that in our convention positive $\sigma$ means compressive strain, whereas negative means tensile. The curves at the lowest stress show a single and sharp transition at the thermodynamic superconducting critical temperature $T_{{\rm c},ac}$, defined as the average temperature of the sharp rise in $C/T$ that occurs when most of the sample becomes superconducting (dashed vertical lines). A single superconducting transition is in agreement with the characterization data of the unstressed samples (S.M.~Fig.~2) and with the results of Refs.~\cite{cairns2020,rosa2021}.

Figures~\ref{Fig1}~(d)~-~(f) show the temperature evolution of the electrical resistivity $\rho_{100}$, $\rho_{110}$ and $\rho_{001}$ with current along the applied stress direction. These measurements were carried out simultaneously with ac calorimetry. For all  stress and current directions, the resistivity at the lowest stress value displays a sharp transition to the superconducting state at a resistive critical temperature $T_{{\rm c},\rho}$ (below which $\rho=0$) that is in good agreement with the one extracted from heat capacity. The small difference between them (less than $0.1\, \rm K$) is easily explained by different definitions of $T_{{\rm c},ac}$ and $T_{{\rm c},\rho}$.

We observe two main effects upon application of compressive uniaxial stress. First,  $T_{{\rm c},ac}$ changes monotonically upon applying stress along all directions. For $\sigma_{100}$ and $\sigma_{110}$, $T_{{\rm c},ac}$ decreases with increasing stress, while $T_{{\rm c},ac}$ increases with increasing $\sigma_{001}$. The evolution of  $T_{{\rm c},ac}$ and $T_{{\rm c},\rho}$ with $\sigma_{100}$, $\sigma_{110}$ and $\sigma_{001}$ is summarized in the stress versus temperature phase diagram shown in Fig.~\ref{Fig2}~(a). As expected for bulk superconductivity, the stress evolution of $T_{{\rm c},\rho}$ tracks that of $T_{{\rm c},ac}$ for most stress values; the difference between $T_{{\rm c},ac}$ and $T_{{\rm c},\rho}$ for $\sigma_{100}>0.250\, \rm GPa$ will be discussed later.

Using Erhenfest's relation, we find that the opposite trend of $T_{{\rm c},ac}$ with $\sigma_{100}$ (and $\sigma_{110}$) as compared to $\sigma_{001}$ is in agreement with a previous thermal expansion study \cite{thomas2021} that reported a negative jump at $T_{\rm c}$ in the linear thermal expansion coefficient along $[100]$ (and $[010]$), but a positive jump for the coefficient along the $[001]$ direction.
These results suggest that the superconducting state observed at ambient pressure is favored by a smaller $c$-axis length  and larger $a$-axis length.
Because a U-U dimer (shortest U-U distance) in the crystal structure of UTe$_{\rm 2}$~is along the $c$ axis, and uranium chains run along the $a$ axis, our results support previous theoretical arguments that the $c$-axis dimer is key to the formation of the superconducting state in UTe$_{\rm 2}$~\cite{shishidou2021,miao2020}.

The second effect of applied uniaxial stress is the appearance of a foot above $T_{{\rm c},ac}$ for $\sigma_{100}$ and $\sigma_{110}$ and a shoulder below  $T_{{\rm c},ac}$~for $\sigma_{001}$ [see Figs.~\ref{Fig1}~(a)~-~(b) and Fig.~\ref{Fig1}~(c), respectively]. As the number of  superconducting transitions in UTe$_{\rm 2}$~remains a central question \cite{rosa2021,hayes2021}, the origin of these features must be understood. It is clear from Figs.~\ref{Fig1}~(a)~-~(b) that the foot above $T_{{\rm c},ac}$ for all $\sigma_{100}$ and $\sigma_{110}$  values coincides with the onset temperature of the superconducting anomaly of the lowest stress curve. For $\sigma_{001}$, the shoulder for all stress values extends above the temperature of the $C/T$ jump at $T_{{\rm c},ac}$ of the lowest stress curve, as shown in Fig.~\ref{Fig1}~(c). The most likely explanation for these features is thus an inhomogeneous stress in the samples, especially in regions of the samples that extend underneath the mounting plates of the stress cell.

To test
 the possibility of inhomogeneous stress, we measured heat capacity at different excitation frequencies $f$. As  $f$ increases, ac calorimetry probes an increasingly smaller, more homogeneously stressed volume of the sample \cite{li2020}.
   The results for $\sigma_{100}=0.050\, \rm GPa$ and $\sigma_{001}=0.200\, \rm GPa$, which are the stress directions that cause the largest effects, are respectively shown in Figs.~\ref{Fig2}~(b) and (c). At the lowest excitation frequency ($f=10\,\rm Hz$), the foot and shoulder features are clearly visible.
    As frequency increases, the position of the main heat capacity jump at $T_{{\rm c},ac}$ remains essentially unchanged. In comparison, the shoulder and foot onset temperatures seem to progressively merge with $T_{{\rm c},ac}$, as expected when probing a more homogeneously stressed part of the sample.
 This shows that stress inhomogeneity is responsible for these two artifacts
  and that there is a \textit{single} heat capacity jump related to bulk superconductivity at the indicated stress levels.

Our main finding
is the absence of splitting of the superconducting transition for applied stress along [100], [110] and [001], as seen in Figs.~\ref{Fig1}~(a)~-~(c).
This is especially meaningful for $\sigma_{110}$ which, unlike $\sigma_{100}$ and $\sigma_{001}$, breaks the orthorhombic symmetry of the crystal to monoclinic by lowering the point group symmetry from D$_{2h}$ to C$_{2h}$. By expressing the stress tensor of $\sigma_{110}$ into the basis of the main crystallographic axes, one finds $\sigma_{110}=0.68\sigma_{100}+0.32\sigma_{010}+0.47\sigma_{xy}$, where $\sigma_{xy}$ is the in-plane shear stress and $\sigma_{010}$ stress along the [010] direction.
Because the resulting shear strain $\epsilon_{xy}$ transforms as the $B_{1g}$ irrep of D$_{2h}$,  it is expected to split the transition temperatures of the two proposed nearly-degenerate OPs \cite{hayes2021}, since their product also transforms as $B_{1g}$. The difference between the two nearly-degenerate superconducting transition temperatures $\Delta T_{\rm c}$ is expected to follow (to leading order in the applied strain):
\begin{equation*}
\Delta T_{\rm c}= \sqrt{\Delta T_{{\rm c}( \epsilon_{xy}=0)}^{2}+\lambda^2 \epsilon_{xy}^2},
\end{equation*}
where $\Delta T_{{\rm c}( \epsilon_{xy}=0)}$ is the unstressed splitting of the superconducting transition temperatures  ($\Delta T_{{\rm c}( \epsilon_{xy}=0)}\approx0$  in our case, since a single transition is observed in the unstressed samples) and $\lambda$ is a coupling constant (see S. M. section E for details of the calculations).

The fact that we do not observe any noticeable splitting of the superconducting transition upon application of $\sigma_{110}$ [see Fig.~\ref{Fig1}~(b)] may
be explained by different scenarios.  One possibility is that the OP is different from the one proposed previously, either because there is only a single superconducting OP or because the nearly-degenerate OPs belong to symmetry channels that would not allow for a coupling that is linear in $\epsilon_{xy}$. Another possible explanation for the absence of detectable splitting of the superconducting transition induced by $\sigma_{110}$ would be a $\lambda$ value too small to cause any appreciable splitting for $\sigma_{110}<0.355\, \rm GPa$.
 Using the elastic tensor obtained from density functional theory (DFT) detailed in S.M. section D,  we find that at that stress level, the induced $\epsilon _{xy}\approx -0.6\%$ would not lead to a $\Delta T_{\rm c}$ greater than $0.1 \, \rm K$ for $\lambda<16\, \rm K$). In this regard, we note that, as detailed in Ref.~\cite{rosa2021}, the presence of two superconducting anomalies in some unstrained samples could be related to a lower sample quality.

If  UTe$_{\rm 2}$~hosts a single superconducting OP, regardless of its symmetry, one expects $T_{\rm c}$ to evolve quadratically with shear strain $T_{\rm c}= T_{\rm c}^{(\epsilon_{xy}=0)}+\lambda \epsilon_{xy}^2$. From Fig.~\ref{Fig2}~(a), we see that $T_{{\rm c},ac}$ evolves linearly with $\sigma_{110}$, like $\sigma_{100}$ and $\sigma_{001}$, but with a different slope. This implies that the $\sigma_{110}$ response is dominated by the symmetry-preserving stress along the main axes, $\sigma_{100}$ and $\sigma_{010}$, rather than the symmetry-breaking stress $\sigma_{xy}$. Symmetry considerations imply that $T_{\rm c}$ should show a linear dependence on $\sigma_{100}$ and $\sigma_{010}$.

  From Fig.~\ref{Fig2}~(a), we determine $\frac{{\rm d} T_{{\rm c},ac} }{ {\rm d} \sigma_{100}}\approx-0.87\,\rm K/GPa$ and $\frac{{\rm d} T_{{\rm c},ac} }{ {\rm d} \sigma_{001}}\approx+0.56\,\rm K/GPa$.

Their sum gives $\sim-0.31\,\rm K/GPa$, which
 upon comparison with both prior  hydrostatic pressure studies \cite{braithwaite2019,thomas2020} ($\frac{{\rm d} T_{\rm c} }{ {\rm d} P}\approx-0.5\,\rm K/GPa$) and thermal expansion and specific heat ($\frac{{\rm d} T_{\rm c} }{ {\rm d} P}\approx-0.49\,\rm K/GPa$) using Erhenfest's relation \cite{thomas2021} suggests that the evolution of $T_{\rm c}$ under applied $\sigma_{010}$ is smaller than that of the two other axes.
 Using the DFT elastic tensor shown in S.M. section D, we find that the evolution of $T_{\rm c}$ with $\sigma_{100}$ and $\sigma_{001}$ cannot be explained in terms of strain along a single direction through Poisson expansion. This suggests that there is no dominant strain direction controlling $T_{\rm c}$ in UTe$_{\rm 2}$.
It would be interesting to apply higher $c$-axis stress, which seems to be tuning the system in a different direction than hydrostatic pressure \cite{braithwaite2019,thomas2020}.

Under hydrostatic pressure, the superconducting transition splits into two above $0.3\, \rm GPa$, leading to a slight slope change of the lower $T_{\rm c} (P)$ and an initial enhancement of the higher $T_{\rm c} (P)$. The latter undergoes a drastic suppression for $P>1.2\, \rm GPa $ \cite{braithwaite2019,thomas2020} due to the emergence of a magnetic ground state. Here, the absence of splitting of the superconducting transition and the linear evolution of  $T_{{\rm c},ac}$ with uniaxial stress approaching $0.3\, \rm GPa$ ($\sigma_{001}$) or exceeding this value ($\sigma_{100}$ and $\sigma_{110}$) suggests that higher stress levels would be required to drive the system to a regime with a different ground state.
However, for $\sigma_{100}>0.250\, \rm GPa$, $T_{{\rm c},\rho}$ starts to increase with increasing stress, in contrast to the behavior of $T_{{\rm c},ac}$ [see Fig.~\ref{Fig2}~(a)]. In addition, the resistive superconducting transition shows substantial broadening upon application of $\sigma_{100}$ whereas its width barely increases with $\sigma_{110}$ and $\sigma_{001}$, as shown in Figs.~\ref{Fig1}~(d)~-~(f). This behavior was verified in another sample with applied $\sigma_{100}$ (see S.M.~Fig.~3).

  For hydrostatic pressures just above $0.3\, \rm GPa$, the emerging superconducting transition that splits from the main $T_{\rm c}(P)$ curve displays a positive ${\rm d} T_{\rm c} / {\rm d} P$ as well as both an initially small signature in heat capacity and a significant broadening in resistivity \cite{braithwaite2019,thomas2020}. A similar scenario could then take place under applied  $\sigma_{100}$, based on the analogous resistive behavior. In this case, a stress value of $\sigma_{100}\approx 0.3 \, \rm GPa$ could be just enough to drive the system towards the regime in which a split superconducting transition emerges. This would mean that shorter U-U distance in the chains along the $a$ axis could be a key ingredient for the enhancement of the second superconducting phase observed at high pressures. Alternatively, this effect may also be caused by the presence of filamentary or surface superconductivity. The application of higher stress would be useful
 to distinguish between the two scenarios, since only in the first one the higher-temperature transition observed in resistivity would be manifested in the heat capacity.

Finally, we turn to the pronounced piezoresistivity  observed above $T_{{\rm c},\rho}$ for stress and current along the $c$ axis [see Fig.~\ref{Fig1}~(f)]. This effect was not observed for stress applied along the other crystal directions [Figs.~\ref{Fig1}~(d)~and~(e)].
By fitting the normal-state $\rho_{001}$ over the extended temperature range shown in Fig.~\ref{Fig2}~(d) to $\rho_{001} = \rho_0 +AT^2$, we observe an
enhancement, upon increasing compressive $\sigma_{001}$, of both the coefficient associated with electron-electron scattering $A$ and the residual resistivity $\rho_0$.
 The observed enhancement of $A$ over a stress range of about $0.5\, \rm GPa$ is non-trivial ($30 \%$)
  but smaller that the approximately factor of two increase in $A$ with applied hydrostatic pressure of $0.56\, \rm GPa$ \cite{braithwaite2019} and much smaller than the $\sim 1000\%$ increase in $A$ at a metamagnetic transition near $ 32\, \rm T$ for magnetic fields applied close to the $b$ axis at atmospheric pressure \cite{knafo2021}.
   These comparisons suggest that at $\sigma_{001}\approx0.5\, \rm GPa$, the system is still away from any possible quantum critical point.
Additionally, we find that the peak in $\rho_{001}$ around $T^ \star\approx15\, \rm K$ that was previously reported in Ref.~\cite{eo2021}, shifts towards lower temperatures as $\sigma_{001}$ increases,
  as displayed in Fig.~\ref{Fig2}~(e). This peak has been attributed to either a Kondo energy scale \cite{kang2021} or
  to the onset of short-range magnetic correlations \cite{willa2021}. The large low temperature $c$-axis piezoresistivity could therefore be caused by the decrease of the energy scale associated with the peak in $\rho_{001}$ at $T^ \star$ with $\sigma_{001}$.

In conclusion, our combined low temperature ac calorimetry and electrical transport measurements under uniaxial stress show that the single thermodynamic superconducting transition at $T_{{\rm c},ac}$ in UTe$_{\rm 2}$~has opposite evolution under applied compressive stress $ \sigma_{100},\sigma_{110}$ ($\frac{{\rm d}T_{{\rm c},ac}}{{\rm d}\sigma_{100}},\frac{{\rm d}T_{{\rm c},ac}}{{\rm d}\sigma_{110}}<0$) and  $ \sigma_{001}$ ($\frac{{\rm d}T_{{\rm c},ac}}{{\rm d}\sigma_{001}}>0$). This result is consistent with conclusions from previous thermal expansion studies and suggests that superconductivity in UTe$_{\rm 2}$~is favored by a smaller $c$-axis length and larger $a$-axis length.
Additionally, we observe a relatively large piezeoresistivity above $T_{\rm c}$ for current and stress applied along the $c$-axis, suggesting that aside from enhancing superconductivity, $\sigma_{001}$ is responsible for a decrease in a different and yet unidentified normal-state energy scale \cite{willa2021,kang2021,eo2021}. Applying higher $ \sigma_{001}$ values could be interesting to determine if the system can be tuned towards a different ground state.
 For $\sigma_{100}>0.25 \, \rm GPa$, we observe a difference between the  $T_{{\rm c},ac}$ and $T_{{\rm c},\rho}$ values, similar to what is observed under hydrostatic pressure exceeding $0.3 \, \rm GPa$, where the signatures of a second superconducting phase become unambiguous \cite{braithwaite2019,thomas2020}.
Our results suggest that a stress level $\sigma_{100}\approx 0.25 \, \rm GPa$ could be just enough to drive the system towards this regime.
Finally, through the application of a symmetry breaking shear stress, $\sigma_{xy}$, we fail to observe clear signatures for two split superconducting transitions, as one would have expected for the case of two nearly-degenerate superconducting OPs whose product transforms as $B_{1g}$. This implies either that the coupling between shear strain and the superconducting OPs is very small or that the superconducting OP of UTe$_{\rm 2}$~is different from the one proposed previously. In the latter scenario, TRS breaking might be explained by the condensation of a sub-leading superconducting instability near dislocations and other lattice defects, similarly to what has been recently proposed to explain TRS breaking in Sr$_2$RuO$_4$ \cite{willa2021_2}. To disentangle these different scenarios, it would be interesting to perform Kerr effect measurements as in Ref.~\cite{hayes2021} on crystals showing a single superconducting transition.

The development of ac calorimetry under uniaxial stress was supported by the Laboratory Directed Research and Development program. The remaining experimental work and crystal synthesis at Los Alamos were supported by the U.S. Department of Energy (DOE), Office of Basic Energy Sciences, Division of Materials Science and Engineering project "Quantum Fluctuations in Narrow Band Systems". JXZ (density functional theory calculations) was supported by Quantum Science Center, a U.S. DOE Office of Science National Quantum Information Science Research Center, and in part by Center for Integrated Nanotechnologies, a U.S. DOE Office of Basic Energies Science user facility, in partnership with the LANL Institutional Computing Program for computational resources. CS and AH (crystal synthesis) acknowledge support from UK EPSRC grant EP/P013686/1. RMF (phenomenological modeling) was supported by the U.S. DOE, Office of Science, Basic Energy Sciences, Materials Science and Engineering Division, under award no. DE-SC0020045.

\bibliographystyle{apsrev4-2}

\section*{Supplementary Material}

\subsection{Methods}

UTe$_{\rm 2}$~single crystals were grown using the chemical vapor transport method described in references~\cite{rosa2021,cairns2020}. Samples used in this study display a single superconducting transition in heat capacity  with $T_{\rm c}\approx 1.8\, \rm K$ (see S.M. Fig.~\ref{UTe2_Cp_comparison}). Bar shaped samples with typical dimensions $\rm 2\, mm \times 0.5 \, mm \times 0.15\, mm$ were prepared, with longest dimensions along the [100], [110] and [001] crystallographic directions using a Laue diffractometer (two different [100] samples and one for each of the other directions were studied,  as shown on S.M. Fig.~\ref{SM_Fig3}). Samples properties are shown in S.M. Table.~\ref{table1}.

 \begin{table}[H]
\caption{Properties of the studied samples: name, unstressed critical temperature $T_{\rm c}$ from heat capacity measurement (from S.M. ~Fig.~\ref{UTe2_Cp_comparison}), sample dimensions, inside and outside distance between the electrical contacts, gap between the stress cell clamps and slope of $T_{{\rm c},ac}(\sigma)$ for the associated uniaxial stress direction.}
\begin{ruledtabular}
\begin{tabular}{cccccc}
Sample name & Unstressed $T_{\rm c}\,$(K) & $l\times w\times t\, (\rm \mu m^3)$ &Contact distance $(\rm \mu m)$ & Gap $(\rm \mu m)$ & $\frac{{\rm d}T_{{\rm c},ac}}{{\rm d}\sigma}\, (\rm K/GPa)$ \\
\hline
100 & $1.85$ & $2307\times462\times170$ &$248-547$& 1168  & -0.87 \\
100 B & 1.66 & $1800\times380\times149$ &$235-615$ &1181 & -0.88 \\
\hline
110 & 1.87 & $3000\times624\times150$ &$264-588$& 2013 & -0.58 \\
\hline
001 & 1.77& $1400\times303\times195$ &$172-404$& 850 & +0.56 \\
\end{tabular}
\end{ruledtabular}
\label{table1}
\end{table}

Uniaxial stress was applied along the length of the samples using a commercial stress cell (Razorbill FC100), resulting in stresses $\sigma_{100}=\sigma_{xx}$, $\sigma_{110}\approx 0.68 \sigma_{xx}+ 0.32 \sigma_{yy}+0.47 \sigma_{xy}$ and $\sigma_{001}=\sigma_{zz}$.  $\sigma_{110}$  was decomposed in the basis of the main axes by rotating the $\sigma_{110}$ stress tensor by an angle $\theta=\arctan{a/b} \approx 34.16^\circ$ in the $ab$ plane. Stress levels reaching $\sigma\approx 0.3-0.4\, \rm GPa$ could be obtained without damaging the samples. Samples were glued between the clamps of the stress cell using Stycast 2850FT epoxy with catalyst 9. Electrical contact to the sample was made using platinum wires and silver paint deposited on gold pads sputtered directly onto the samples. Electrical resistance was measured using a Lakeshore 372 resistance bridge. Temperature oscillations for the AC heat capacity measurements were generated using a constantan wire heater ($R\approx \rm 15 \, \Omega$) using powers of order $P_{AC}\approx \rm 1-10 \, \mu W$ at frequencies of $2f\approx \rm 10-320 \, Hz$.  The induced temperature oscillations $T_{AC}$ were recorded using a chromel-Au/Fe(0.07\%) thermocouple and a SR860 lockin amplifier after being amplified  by a SR554 low-noise transformer preamplifier. Both the heater and thermocouple were mechanically and thermally connected to the samples using GE varnish. The experimental setup is shown in S.M. Fig.~\ref{experimental_setup}. In this configuration, the heat capacity of the system is related to $C\propto P_{AC}/\left\vert T_{AC}\right \vert f$. All measurements were carried out either in a CMR adiabatic demagnetization refrigerator, allowing measurements down to $0.15\, \rm K$, or in a Quantum Design PPMS for sample characterization and measurements above $3\, \rm K$.

The thermodynamic critical temperature $T_{{\rm c},ac}$ of the samples was determined as the average temperature of the sharp rise of $C/T$ occurring when most of the sample becomes superconducting. This sharp rise can be distinguished from the smoother foot above $T_{{\rm c},ac}$ observed for applied $\sigma_{100}$ and $\sigma_{110}$.

Applied stress level was determined using the capacitance versus force calibration curve of the stress cell. Force was then converted into stress by dividing it by the cross section of the sample perpendicular to the applied stress direction. The dimensions of the samples, shown in S.M. Table.~\ref{table1}, were measured using an optical microscope, with a typical resolution of $10\, \mu m $.

As stress is either increased or decreased, we noticed slight changes in the thermal coupling between the heater/thermocouple and the samples. This results in a (small) decrease in the temperature oscillations (with the same excitation current) on measurement carried out consecutively at different stress levels. This effect explains the small increase of the normal state value of $C/T$ with stress observed in Main Fig.~1~(a), for which data were recorded with decreasing $\sigma_{100}$, and the opposite effect for Main Figs.~1~(b)~-~(c), for which data were recorded with increasing $\sigma_{110}$ and $\sigma_{001}$.

\begin{figure}[H]
\centering
\includegraphics[width=0.4\textwidth]{SM_Fig1.pdf}\caption{Photograph of the experimental setup showing the sample (in black and black dashed lines) held with epoxy between the cell clamps (grey). Solid and dashed blue lines indicate the location of the leads for resistivity measurement. Red lined show the constantan heater and its current leads deposited on top of the sample using GE varnish. Yellow and white lines respectively show the Au/Fe(0.07\%) and chromel thermoucouple wires connected underneath the sample using GE varnish.}
\label{experimental_setup}
\end{figure}

\subsection{Unstressed characterization}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{SM_Fig2.pdf}
\caption{ \textbf{(a)} Specific heat data of the three samples showed in the paper and of the additional 100 B sample measured before the samples were mounted in the stress cell. Dashed lines indicated the $T_{\rm c}$ values defined as the average temperature of the superconducting jump. Dotted lines represent the normal state specific heat. \textbf{(b)} Temperature dependence of the resistivity $\rho_{100}$, $\rho_{110}$ and $\rho_{001}$ of our unstressed samples 100, 110 and 001.}
\label{UTe2_Cp_comparison}
\end{figure}

The residual resistivity above $T_{{\rm c},\rho}$, $\rho_0$, has values of the order of $10\, \rm \mu \Omega\, cm$ for $\rho_{100}$ and $\rho_{110}$, and $100\, \rm \mu \Omega\, cm$ for $\rho_{001}$. These values are consistent with previous reports of the resistivity of UTe$_{\rm 2}$~along the different crystal directions \cite{ran2019,eo2021
.  Although we do not find that $\rho_{100}< \rho_{110} $, this is probably due to a fairly large uncertainty of the geometric factors of our samples. In order to measure resistivity over a more homogeneously stressed region of the sample, the voltage leads are centered in the gap between the mounting plates and placed close together ($150-600\, \rm \mu m$).

\subsection{Additional 110 sample}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{SM_Fig3.pdf}
\caption{\textbf{(a)} Temperature dependence of the heat capacity $C/T$ (top panel) and electrical resistivity $\rho_{100}$ (bottom panel) of sample 100 B at the indicated  $\sigma_{001}$ values. Colored dashed lines mark the average temperature of the sharp rise of the jump in $C/T$ at $T_{{\rm c},ac}$  and the temperature $T_{{\rm c},\rho}$ below which the resistivity of the sample is zero.
 \textbf{(b)} Stress dependence of $T_{{\rm c},\rho}$ (circles) and $T_{{\rm c},ac}$ (crosses) for our samples 100 (red) and 100 B (maroon). Circles represent $T_{{\rm c},\rho}$ and crosses $T_{{\rm c},ac}$ extracted from a). Colored dashed lines are linear fits $T_{{\rm c},ac}(\sigma)$, which slopes $\frac{{\rm d}T_{{\rm c},ac}}{{\rm d}\sigma}$ are displayed on the lower left corner of the figure.}
\label{SM_Fig3}
\end{figure}

\subsection{DFT calculation of the elastic constants}
The elastic constants were calculated by using density functional theory (DFT) with the pseudopotential projector-augmented wave method \cite{SM__Kresse1999} implemented in the Vienna Ab initio Simulation Package (VASP) \cite{SM__Kresse1996,SM__Kresse1993}.  A local density approximation to the exchange-correlation functional combined with the correlation correction \cite{SM_Liechtenstein1995} with $U=6 \, \rm eV$ and $J=0.57 \, \rm eV$ to the U-5f orbitals was used. The spin-orbit coupling was incorporated. In the nonmagnetic constraint, the simulations give the fully relaxed lattice constants $a = 4.11\,  \mbox{\normalfont\AA}$, $b = 6.13\,  \mbox{\normalfont\AA}$, $c = 13.93\,  \mbox{\normalfont\AA}$, in good agreement with the experimental measurements at ambient conditions \cite{SM_Boehme1992}. Further relaxation calculations were then performed for varying uniaxial strains applied along $a$- and $c$-axis. For each structure, the elastic moduli were then evaluated. The obtained DFT elastic tensor (in GPa) is:

 \begin{equation*}
 \begin{bmatrix}
\sigma_{xx} \\
\sigma_{yy} \\
\sigma_{zz} \\
\sigma_{xy} \\
\sigma_{yz} \\
\sigma_{zx} \\
\end{bmatrix}
=
 \begin{bmatrix}
96.96782    & 40.94466 &48.58602     & 0 & 0 & 0    \\
40.94466       & 140.74513     &46.67250           & 0 & 0 & 0    \\
48.58602           & 48.58602         &101.32160                & 0 & 0 & 0    \\
0           & 0         &0                & 27.06222      & 0 & 0    \\
0           & 0         &0                & 0      & 19.61040      & 0    \\
0           & 0         &0                & 0      & 0 & 57.25525    \\
\end{bmatrix}
 \begin{bmatrix}
\epsilon_{xx} \\
\epsilon_{yy} \\
\epsilon_{zz} \\
\epsilon_{xy} \\
\epsilon_{yz} \\
\epsilon_{zx} \\
\end{bmatrix}.
\end{equation*}

\subsection{Phenomenological model for the effect of shear strain}

Consider two nearly degenerate superconducting instabilities, with
gaps $\Delta_{\gamma_{1}}$ and $\Delta_{\gamma_{2}}$, that transform
as different one-dimensional irreducible representations $\gamma_{1}$ and $\gamma_{2}$ of the D$_{2h}$ point group.
The Kerr data \cite{hayes2021} indicates that the quantity $i\left(\Delta_{\gamma_{1}}\Delta_{\gamma_{2}}^{*}-\Delta_{\gamma_{1}}^{*}\Delta_{\gamma_{2}}\right)$
has the same symmetry as an out-of-plane magnetic field $H_{z}$,
which transforms as $B_{1g}$ and is odd under time-reversal. As discussed
in Ref. \cite{hayes2021}, this restricts $\gamma_{1}$ and $\gamma_{2}$
to the following combinations: $\left(B_{3u},B_{2u}\right)$, $\left(B_{1u},A_{u}\right)$,
$\left(B_{3g},B_{2g}\right)$, $\left(B_{1g},A_{g}\right)$. Thus,
for any of these combinations, the time-reversal even quantity $\left(\Delta_{\gamma_{1}}\Delta_{\gamma_{2}}^{*}+\Delta_{\gamma_{1}}^{*}\Delta_{\gamma_{2}}\right)$
transforms as $B_{1g}$, i.e. as shear strain $\epsilon_{xy}$. To
quadratic order in $\Delta_{\gamma_{i}}$ and to leading order in
$\epsilon_{xy}$, the Landau free-energy expansion is:

\begin{equation}
F=a_{1}\left|\Delta_{\gamma_{1}}\right|^{2}+a_{2}\left|\Delta_{\gamma_{2}}\right|^{2}+\tilde{\lambda}\epsilon_{xy}\left(\Delta_{\gamma_{1}}\Delta_{\gamma_{2}}^{*}+\Delta_{\gamma_{1}}^{*}\Delta_{\gamma_{2}}\right)
\end{equation}

Here, $\tilde{\lambda}$ is a coupling constant, $a_{i}=a_{i,0}\left(T-T_{c,i}^{(0)}\right)$
with $a_{i,0}>0$ and $i=1,2$, and $T_{c,i}^{(0)}$ are the unstressed
transition temperatures associated with the two superconducting instabilities.
After minimizing the relative phase between the two order parameters,
we can rewrite this expression in the matrix form:

\begin{equation}
F=\left(\begin{array}{cc}
\left|\Delta_{\gamma_{1}}\right| & \left|\Delta_{\gamma_{2}}\right|\end{array}\right)\left(\begin{array}{cc}
a_{1} & -\left|\tilde{\lambda}\epsilon_{xy}\right|\\
-\left|\tilde{\lambda}\epsilon_{xy}\right| & a_{2}
\end{array}\right)\left(\begin{array}{c}
\left|\Delta_{\gamma_{1}}\right|\\
\left|\Delta_{\gamma_{2}}\right|
\end{array}\right)
\end{equation}

Diagonalization of the matrix gives two superconducting transition
temperatures given by the condition $a_{1}a_{2}=\left|\tilde{\lambda}\epsilon_{xy}\right|$.
We find:

\begin{align}
T_{c,1} & =\left(\frac{T_{c,1}^{(0)}+T_{c,2}^{(0)}}{2}\right)+\sqrt{\left(\frac{T_{c,1}^{(0)}-T_{c,2}^{(0)}}{2}\right)^{2}+\frac{\lambda^{2}\epsilon_{xy}^{2}}{4}}\\
T_{c,2} & =\left(\frac{T_{c,1}^{(0)}+T_{c,2}^{(0)}}{2}\right)-\sqrt{\left(\frac{T_{c,1}^{(0)}-T_{c,2}^{(0)}}{2}\right)^{2}+\frac{\lambda^{2}\epsilon_{xy}^{2}}{4}}
\end{align}
where we defined $\lambda\equiv2\tilde{\lambda}/\sqrt{a_{1,0}a_{2,0}}$.
Thus, the splitting $\Delta T_{c}$ between the two transitions depends
on $\epsilon_{xy}$ according to:

\begin{equation}
\Delta T_{c}=\sqrt{\left(\Delta T_{c}^{(0)}\right)^{2}+\lambda^{2}\epsilon_{xy}^{2}}
\end{equation}
where $\Delta T_{c}^{(0)}$ is the splitting between the two transitions
in the absence of applied stress.

\bibliographystyle{apsrev4-2}


\section{Introduction}
\label{introduction}

Machine learning has risen as an important way for modeling and control of complex dynamic systems to accelerate challenging robotics applications~\cite{Nguyen-Tuong:CP:2011,Jordan:Science:2015,Pierson:AR:2017}. This emerging field has witnessed advances, especially based on two frameworks. The first framework exploits reinforcement learning to train robots such that they learn optimal control policies from experiences to accomplish certain goals~\cite{Kober:IJRR:2013}. Despite successes in tackling various tough problems, this framework in general demands hefty amounts of training data and offers difficult generalization to tasks different from specified goals~\cite{Kober:IJRR:2013}. By contrast, the second framework seeks to integrate machine learning-based modeling with optimal control~\cite{Nguyen-Tuong:CP:2011}. Specifically, it extracts data-driven models for robotic systems and then synthesizes optimal control using such explicit models. The corresponding robot control methods would arguably present higher data efficiency and generalizability. Associated with this, a still open question is how to enable control design that well fits with machine learning models, which continually calls for the development of constructive methods.

In this study, we consider the problem of controlling neural network models for dynamic systems and investigate nonlinear model predictive control (NMPC) for them. Neural networks have gained increasing use in modeling robots with complicated dynamics. However, the application of NMPC to them is non-trivial. Conventional NMPC entails online constrained optimization, but numerical optimization of neural network dynamics can be extremely burdensome due to a mix of nonlinearity, nonconvexity, and heavy computation. This is particularly true when a neural network has many hidden layers and nodes. In addition, some robotics applications involve cost functions that have zero gradients, precluding the use of gradient-based optimization.

To address the challenge, we propose to leverage a new NMPC approach, which was developed in our previous study~\cite{Askari:ACC:2021} and named {\em constraint-aware particle filtering/smoothing NMPC} (CAP-NMPC), to control neural network dynamics. Departing from the optimization-based view, the CAP-NMPC interprets NMPC through the lens of Bayesian estimation and designs a constrained particle filtering/smoothing method to achieve it. This approach, at its core, uses sequential Monte Carlo sampling to estimate the optimal control actions from a reference signal that is to be followed by the neural network model over a receding horizon. The sampling-based implementation is arguably advantageous in several ways for control of neural network dynamics. It can make sufficient search for the best control actions across a large control space when using adequate numbers of particles.
The sampling-based computation also obviates the need for iteratively computing gradients and is relatively easy and efficient to achieve.

We further consider the problem of motion planning for autonomous vehicles. We construct a neural network model to capture the vehicle dynamics and then use the CAP-NMPC to perform motion planning in dynamic driving scenarios (e.g., lane changing and moving obstacles). The simulation results demonstrate the effectiveness of the design.

\section{Related Work}
\label{literature}

Neural networks have proven useful for data-driven modeling of many systems that resist analytical modeling or suffer significant modeling uncertainty.
A common approach is using feedforward neural networks to approximate the state dynamics functions. While simple network structures can be enough for some systems~\cite{Draeger:CSM:1995,Piche:NIPS:1999,Williams:ICRA:2017, Williams:TRO:2018}, the recent development of deep neural networks has allowed to capture more complex dynamics~\cite{Broad:arXiv:2018,Nagabandi:ICRA:2018}. Further, recurrent neural networks can effectively learn closed-loop or residual dynamics~\cite{Garimella:2018,Rankovic:IJCCC:2012}.

	Based on the prediction by a neural network model, NMPC can be readily designed to predictively optimize a system's behavior, as is pursued in~\cite{Draeger:CSM:1995,Piche:NIPS:1999,Williams:ICRA:2017, Williams:TRO:2018,Broad:arXiv:2018,Nagabandi:ICRA:2018,Garimella:2018,Rankovic:IJCCC:2012}. Yet, numerical optimization in this setting is recognized as a thorny issue~\cite{Williams:TRO:2018,Broad:arXiv:2018}, as gradient-based methods often find themselves inadequate and computationally expensive in the face of nonconvex optimization of neural network dynamics. Gradient-free NMPC is thus desired, and the sampling-based approach holds significant promise. A simple random-sampling shooting method is used in~\cite{Nagabandi:ICRA:2018} to treat NMPC of a deep neural network model, which generates control action sequences randomly and chooses the sequence that leads to the highest expected cumulative reward. The study in~\cite{Williams:TRO:2018} formulates information-theoretic NMPC for neural network models and then synthesizes an algorithm based on iteratively weighting sampled sequences of control actions. A direct estimation approach to motion planning for autonomous vehicles is presented in~\cite{Berntorp:TIV2019}, where the planning requirements are modeled as measurements and the trajectory of the vehicle is inferred using a particle filter. The CAP-NMPC approach differs from~\cite{Nagabandi:ICRA:2018,Williams:TRO:2018, Berntorp:TIV2019} on three aspects. First, it is developed from the perspective of Bayesian estimation and builds on particle filtering/smoothing to solve the NMPC problem, which is different that the direct filtering approach in~\cite{ Berntorp:TIV2019}. This feature also allows for diverse realizations of the CAP-NMPC since the literature includes a rich set of particle filtering/smoothing techniques. Second, rather than randomly creating all the control action samples at once and then evaluating their competence (e.g., cumulative rewards or weights), this approach generates samples sequentially based on the neural network model and then weights and resamples them. This principled manner has the potential to search the control space more effectively. Third, by design, the CAP-NMPC approach takes generic state and input constraints into account, as needed for robotic systems subject to operating constraints.

A critical function of autonomous driving, motion planning has attracted enormous interest in the past decade. State of the art has three main approaches: input space discretization with collision checking, randomized planning, and NMPC~\cite{Schwarting:ARCRAS:2018,Gonzalez:TITS:2016,Paden:TIV:2016}. Among them, NMPC can optimize vehicle motion using sophisticated vehicle models under practical constraints, thus capable of planning and executing safety-critical, aggressive (near-limit), or emergency maneuvers in complex environmental conditions~\cite{Paden:TIV:2016}. NMPC-based motion planning has gained a growing body of work recently, e.g.,~\cite{Williams:TRO:2018,Chen:TIV:2019,Brunner:CDC:2017,Nolte:IVS:2017,Cardoso:ICRA:2017,Liu:IVS:2017,Gao:DSCC:2010}, just to highlight a few. However, almost all studies, except~\cite{Williams:TRO:2018}, consider physics-based vehicle models to our knowledge, even though
modeling errors or biases are inevitable in the real world. Neural network vehicle models have demonstrated substantial achievements recently~\cite{Spielberg:SR:2019}. This paper will study motion planning based on them, using the CAP-NMPC approach as the enabling tool.

\section{The {CAP-NMPC} Approach}

In this section, we provide an overview of the CAP-NMPC approach. A more detailed description is available in~\cite{Askari:ACC:2021}.

\subsection{NMPC through the Lens of Bayesian Estimation}
\label{NMPC-Bayesian-view}

Consider a nonlinear dynamic system of the form:
\begin{align}\label{state-equation}
x_{k+1} = f(x_{k},u_{k}),
\end{align}
where $x_{k} \in \mathbb{R}^{n_x}$ is the system state, and $u_{k} \in \mathbb{R}^{n_u}$ is the control input. The nonlinear mapping $f:  \mathbb{R}^{n_x} \times \mathbb{R}^{n_u} \rightarrow \mathbb{R}^{n_x}$ characterizes  the state transition. Here, it  represents  a feedforward neural network   model learned from data made on the actual system, or a hybrid model combining  physics with neural networks.
 The system is subject to the following inequality constraints:
\begin{align}\label{inequality-constraints}
g_j({x}_k,u_{k})\leq 0,  \ \forall j=1,\ldots,m,
\end{align}
where $m$ is the total number of constraints.  We consider an   NMPC problem  such that $x_k$ tracks a reference signal $r_k$, which is stated as follows:
\begin{subequations} \label{NMPC-Standard}
\begin{align}
\min_{u_{k:k+H}} \quad &  \sum_{t=k}^{k+H}
  u_t^\top  Q u_t + ( x_t -  r_t)^\top R  ( x_t -  r_t),\\
\mathrm{s.t.} \quad & x_{t+1} = f(x_{t},u_{t}) , \\ \label{ineq-constraint}
& g_j(x_t,u_t) \leq 0 \quad \forall j=1,\ldots,m,\\ \nonumber
& t = k, \ldots, k+H,
\end{align}
\end{subequations}
where $H$ is the length of the upcoming horizon,  $u_{k:k+H} = \left\{u_k, u_{k+1}, \ldots, u_{k+H}\right\}$, and  $Q$ and $R$  are weighting matrices.  The above problem is solved   through time in a receding-horizon manner by computing   the optimal control input sequence $u_{k:k+H}^\ast$.  At every time, the first element, $u_k^\ast$  will be applied to control the system, and the rest  discarded. The same   optimization and control procedure will repeat itself recursively at the future time instants.

While the NMPC problem in~\eqref{NMPC-Standard} is usually solved through numerical optimization in the literature, it can be interpreted as a Bayesian  estimation problem and thus addressed.  The following theorem shows this connection~\cite{Askari:ACC:2021,Stahl:SCL:2011}.
\begin{thm}\label{ML-NMPC-Equivalence}
For the horizon $t=k,\ldots,k+H$, consider the virtual system
\begin{align}\label{Virtual-System}
\left\{
\begin{aligned}
x_{t+1} &= f(x_t, u_t),\\
u_{t+1} &=  w_t,\\
r_t &= x_t+ v_t,
\end{aligned}
\right.
\end{align}
where   $w_t \sim \mathcal{N} \left(0,Q^{-1} \right)$, $v_t \sim \mathcal{N} \left(0, R^{-1} \right)$, and the reference $r_t$ serves as virtual measurements. The maximum likelihood estimation of $x_t$ and $u_t$ via
\begin{align}\label{ML-Formulation}
\max_{x_{k:k+H}, u_{k:k+H}} \log \p(x_{k:k+H}, u_{k:k+H}   |  r_{k:k+H}),
\end{align}
where $x_{k:k+H} = \left\{x_k,x_{k+1}, \ldots, x_{k+H}\right\}$ (similarly for $u_{k:k+H} $ and $r_{k:k+H} $),  is equivalent to the NMPC problem in~\eqref{NMPC-Standard} without the inequality constraints in~\eqref{ineq-constraint}.
\end{thm}

Theorem~\ref{ML-NMPC-Equivalence} indicates that the original   NMPC problem can be converted into the problem of {\em estimating} the optimal control actions along with states given the specified reference. This viewpoint then ushers a new way of treating NMPC through estimation methods.

Before proceeding further, we explain how to incorporate the inequality constraints into a Bayesian estimation procedure, since they are an integral part of the NMPC formulation. Here, we adopt the barrier function method   to create  virtual measurements about the constraint satisfaction:
\begin{align}\label{Virtual-measurement-barrier-function}
z_t = \phi\left(g( {x}_t, u_t)\right) + \eta_t,
\end{align}
where $z$ is the virtual measurement variable, $g$ is the collection of $g_j$ for $j=1,\ldots, m$, $\eta$ is an additive small noise, and $\phi$ is a barrier function.
Here, $\phi$ is chosen to be the softplus function:
\begin{equation}
 \phi{\left( s \right)} = \frac{1}{\alpha}\ln\left(1+\exp(\beta s )\right),
\end{equation}
which is parametrized by two tuning factors $\alpha$ and $\beta$.
Note that $\phi$ is  fully continuous  and   through appropriate parameterization, outputs almost zero at a point within the constraint set and   large values at points outside the set.  Hence,  $\phi$ can be used to quantify the constraint satisfaction or violation, and we let the virtual measurement $z_t$  be $0$. Further, we can include $z_t$ into~\eqref{Virtual-System} to allow estimation with an awareness of constraints.

Now, consider~\eqref{Virtual-System} along with $z_t$ and rewrite it compactly in an augmented form:
\begin{align}\label{Virtual-System-Compact}
\left\{
\begin{aligned}
\bar{x}_{t+1} &= \bar{f}(\bar{x}_t) + \bar{w}_t, \\
\bar r_t &=\bar h(  \bar{x}_t ) +\bar v_t,
\end{aligned}
\right.
\end{align}
for $t = k, \ldots, k+H$,
where $\bar{x}_t=\left[x_{t}^{\top} \ u_{t}^{\top}\right]^{\top}$, $\bar{w}_t=\left[0^\top \ w_t^{\top}\right]^{\top}$,  $\bar{v}_t=\left[v_t^\top \ \eta_t^{\top}\right]^{\top}$, $\bar f $ stems from $f$, and  $\bar h$ results from $r_t$ and $z_t$.  Based on the above, to address the original NMPC problem in~\eqref{NMPC-Standard}, we only need to perform state estimation for the above augmented system via
\begin{align} \label{Bayesian-Estimation-Augmented-System}
\max_{\bar x_{k:k+H}} \log \p( \bar x_{k:k+H}  |  \bar r_{k:k+H}).
\end{align}
 This can be achieved by particle filtering/smoothing, which is known as an effective means of state estimation for even highly nonlinear systems.

\subsection{Development of the {CAP-NMPC} Approach} \label{CAP-NMPC-deployment}

The state estimation problem in~\eqref{Bayesian-Estimation-Augmented-System} involves both filtering and smoothing.
We first look at the forward filtering by considering $\p (\bar x_{k:t} | \bar r_{k:t})$ for $k \leq t \leq k+H$. Because it is practically impossible to obtain an analytical expression  of $\p (\bar x_{k:t} | \bar r_{k:t})$ for nonlinear systems, we  approximate it by using  a sample-based empirical distribution.  A common and useful approach is to do  importance sampling. That is, one draws samples from an alternative known distribution $\q (\bar x_{k:t} | \bar r_{k:t})$, which is called   importance or proposal distribution, and then evaluate  the weights of the samples in relation to $\p (\bar x_{k:t} | \bar r_{k:t})$. Suppose that $N$ samples, $\bar x_{k:t}^i$ for $i=1,\ldots,N$, are drawn from $\q (\bar x_{k:t} | \bar r_{k:t})$. Their importance weights are  given by
\begin{align}\label{Weights}
W_t^i = \frac{\p (\bar x_{k:t} | \bar r_{k:t})}{\q (\bar x_{k:t} | \bar r_{k:t})},
\end{align}
which are then normalized to be between 0 and 1. As such, $\p (\bar x_{k:t} | \bar r_{k:t})$ can be  approximated as
\begin{align*}
\p (\bar x_{k:t} | \bar r_{k:t}) \approx \sum_{i=1}^N W_k^i \delta \left( \bar x_{k:t} - \bar x_{k:t}^i\right).
\end{align*}
Note that~\eqref{Weights} also implies a recurrence relation in the  weight update:
\begin{align*}
W_t^i = \frac{\p (\bar x_{k:t} |\bar r_{k:t})}{\q (\bar x_{k:t} | \bar  r_{k:t})} = \frac{\p ( \bar r_{t} | \bar x_t^i) \p(\bar x_t^i | \bar x_{t-1}^i)}{\q ( \bar x_t^i | \bar x_{t-1}^i, \bar r_{k:t} )} W_{t-1}^i.
\end{align*}

One has different ways to implement the above procedure, with the key lying in choosing the importance distribution $q$. A  straightforward choice is to let $\q ( \bar x_t | \bar x_{t-1}, \bar r_{k:t} ) =  \p (\bar x_t | \bar x_{t-1})$. Given this choice,  we can draw samples  $\bar x_t^i \sim \p(\bar x_t | \bar x_{t-1}^i)$ at time $t$ and compute the associated normalized weights via
\begin{align}\label{const_weight}
W_t^i = \frac{\p(\bar r_t | \bar x_t^i)}{\sum_{j=1}^N \p(\bar r_t | \bar x_t^j)}.
\end{align}
This implementation is called  the {\em bootstrap particle filter}. For a particle filtering run,  a majority of the particles may have zero or almost zero weights after a few time steps. This is known as the issue of particle degeneracy, which lowers the overall quality of the particles   and reduces the estimation performance.  To resolve this issue, resampling can be used to replace  low-weight particles by those with high weights \cite{Sarkka:CUP:2013}.

The backward smoothing follows the forward filtering as we only require $\p(\bar x_k | \bar r_{k:k+H})$, which can lead to a more accurate estimation of $\bar x_k$. There are different particle smoothers, and  the reweighting particle smoother will suffice here. It reweights the samples in a recursive  backward manner via
\begin{align}\label{reweighting}
W_{t|k+H}^i = \sum_{j=1}^N W_{t+1|k+H}^i\frac{W_t^i \p(\bar{x}_{t+1}^j|\bar{x}_t^i)}{ \sum_{l=1}^N W_t^l \p(\bar{x}_{t+1}^j|\bar{x}_t^l)},
\end{align}
where   $W_{k+H|k+H}^i = W_{k+H}^i$. The resultant procedure is called  {\em reweighted particle smoother}. After the smoothing,  the   empirical distribution for $\p(\bar{x}_k| \bar r_{k:k+H})$ is given by
\begin{equation*}
\p(\bar{x}_k| \bar r_{k:k+H})\approx\sum_{i=1}^N W_{k|k+H}^i\delta(\bar{x}_k-\bar{x}_k^i).
\end{equation*}
Finally, the best estimate of $\bar x_k$ from $\bar r_{k:k+H}$ is
\begin{align}\label{estimate}
\hat{\bar x}_k^\ast =  \mathbb{E}\left(\bar x_k \, | \, \bar r_{k:k+H}\right) = \sum_{i=1}^N W_{k|k+H}^i \bar x_k^i,
\end{align}
from which   the optimal control input  $u_k^*$  can be read.

Summarizing the above, the CAP-NMPC approach can be outlined as in Algorithm~\ref{NMPC-CA-PF/S}. It presents a fully sampling-based implementation of NMPC, which is promising for control of neural network dynamics. It should also be noted that this approach can well admit other realizations of particle filtering/smoothing, depending on the needs for sophistication and accuracy of estimation.

\begin{algorithm}[t]
\fontsize{10}{10}
  \caption{CAP-NMPC: NMPC via Constraint-Aware Particle Filtering/Smoothing}	  \label{NMPC-CA-PF/S}
  \begin{algorithmic}[1]
\State {Set up   NMPC}   by specifying the dynamic system~\eqref{state-equation}

\State {Recast NMPC} as particle filtering/smoothing by setting up the virtual system~\eqref{Virtual-System-Compact}

\For{$k=1, \ldots, T$}

\item[]{\hspace{14pt}\bf\em    Forward filtering}

\For{$t = k,\ldots, k+H$}

\If{$t=k$}

\State Draw samples $\bar x_k^i \sim  p(\bar x_k)$, $i=1, \ldots, N$

\Else

\State Draw samples $\bar x_t^i \sim \p(\bar x_t | \bar x_{t-1}^i)$, $i=1,\ldots,N$

\State Evaluate sample weights via~\eqref{const_weight}

\State Do resampling based on the weights

\EndIf

\EndFor

\item[]{\hspace{14pt}\bf\em    Backward smoothing}

\For{$t = k+H,\ldots, k$}

\If{$t=k+H$}

\State Assign $W_{k+H|k+H}^i = W_{k+H}^i$, $i=1,\ldots,N$

\Else

\State Reweight the particles via~\eqref{reweighting}

\EndIf

\EndFor

\State Compute the optimal estimation of $\bar x_k$ via~\eqref{reweighting}

\State Export $u_k^\ast$ from $\bar x_k^\ast$, and apply it to the system~\eqref{state-equation}

\EndFor

\iffalse
\State Draw initial set of particles $\bar{x}_k^i$ for $i=1,\ldots,N$ from the prior $p(\bar{x}_k)$
\State{Repeat the follwoing steps for $j=k,\ldots,k+H$:}
\Statex
\begin{itemize}
\item Draw samples $\bar{x}_j^i\sim p(\bar{x}_{j}|\bar{x}_{j-1})$ for $i=1,\ldots,N$.
\item Compute the VO $c_j$ from Eq. (11) for $i=1,\ldots,N$.
\item Assign the particle weight $\tilde{w}_j^i$ using Eq. (13) for $i=1,\ldots,N$.
\item Normalize weights $w_j^i=\frac{\tilde{w}_j^i}{\sum_{i=1}^Nw_j^i}$.
\item Perform resampling based on the weights $w_j$
\end{itemize}
\State Initialize the smoother weights $w_{k+H|k+H}^i$ by the filtering weight $w_{k+H}^i$ for $i=1,\ldots,N$.
\State Repeat the following for $j=k+H,\ldots,1$
\begin{itemize}
\item Assign the smoother weight using Eq. (10).
\end{itemize}
\State Obtain the control input $u_k=\sum_{i=1}^N w_{k|k+H}u_k^i$, where $u_k^i$ is the second element of the augmented vector $\bar{x}_k^i$.
\Output{$u_k$}
\fi
\end{algorithmic}
\end{algorithm}

\section{Motion Planning Using Neural Network Dynamics \textit{via} CAP-NMPC} \label{motion-planning-formulation}

In this section, we present the autonomous vehicle motion planning problem.

\subsection{Motion Planning Problem Formulation}

We consider an autonomous vehicle whose goal is to arrive at the desired goal state $x^g$ from an initial state $x_0$   while considering its own dynamics~\eqref{state-equation} and constraints imposed by the surrounding~\cite{Chen:TIV:2019}. The objective of planning is to find an optimal trajectory without violating the constraints. To this end, an NMPC problem can be formulated as
\begin{subequations} \label{Motion-planning-NMPC}
\begin{align}
\min_{u_{k:k+H}} \quad &  \sum_{t=k}^{k+H}
  u_t^\top  Q u_t + ( x_t -  x^g)^\top R  ( x_t -  x^g) ,\\ \label{motion-planning-model}
\mathrm{s.t.} \quad & x_{t+1} = f(x_{t},u_{t}) , \\ \label{motion-planning-inequality-constraints}
& g_j(x_t,u_t) \leq 0 \quad \forall j=1,\ldots,m,\\ \nonumber
& t = k, \ldots, k+H,
\end{align}
\end{subequations}
which follows~\eqref{NMPC-Standard} by assigning the goal state to be the reference.  The model in~\eqref{motion-planning-model} captures the vehicle dynamics, and~\eqref{motion-planning-inequality-constraints} encompasses all the constraints due to vehicle operation or planning scenarios.

A mission of motion planning   is to make the autonomous vehicle  perform as well as or even better, than human drivers in a wide array of conditions or situations, and a key is using capable models. Despite the utility, physics-based  vehicle models have practical limitations, which may contain errors or biases, or suit only certain specific conditions. Meanwhile, the increasing abundance of data generated by autonomous vehicles makes it possible to construct precise, broadly applicable data-driven models. Among them, neural networks have shown tremendous merits for model-based control of vehicles~\cite{Spielberg:SR:2019}. Therefore, we will leverage a neural network model for the considered motion planning problem, with more details offered in~\ref{neural-network-dynamics}. Further,  a sampling-based implementation of  NMPC better suits control of neural networks than numerical optimization, as argued in Sections~\ref{introduction}-\ref{literature}. Recent studies also suggest sampling can be more competitive than numerical optimization in dealing with large-scale problems~\cite{Ma-Chen-Flammarion-Jordan:PNAS:2019}. We hence will use the CAP-NMPC approach to address the motion planning problem~\eqref{Motion-planning-NMPC} with neural network dynamics.

\subsection{The Neural Network Model}
\label{neural-network-dynamics}

A feedforward neural network is utilized to represent the vehicle dynamics. We describe the neural-network-based parameterization of  the vehicle system  and then explain the data collection and training process.

\subsubsection{Vehicle Dynamics Parameterization}\label{sec:dynamic-system-parameterization}
There are different approaches to parameterize the vehicle dynamics using neural networks.  A simple approach is to feed the state $x_{k}$ and control $u_{k}$ to a neural network  and make the network predict the next state of the system $x_{k+1}$. However, learning such a function could be difficult due to the small sampling time $\Delta t$. This leads to $x_k$ and $x_{k+1}$ being very similar in making, causing the  network to effectively learn an identity transformation~\cite{Nagabandi:ICRA:2018}.  Hence, we instead focus on parameterizing the state transition function, as it shows the incremental changes in the state. The  parametrization is then given by
\begin{align*}
x_{k+1} = x_{k} + \Delta t\hat{f}(x_{k}, u_{k}; \theta),
\end{align*}
where $ \hat{f}$ represents the neural network to approximate the vehicle dynamics, and $\theta$ collects the weights of the network.

\subsubsection{Data collection and pre-processing}\label{sec:data-collection}

We generate and collect the training data from the considered vehicles model.
This is attained by sampling a set of states from a uniform distribution $x_{k} \sim \mathrm{Uni}(\mathcal{X}_k)$ and applying a random control, which is also sampled from a uniform distribution $u_{k} \sim \mathrm{Uni}(\mathcal{U}_k)$. Here, $\mathcal{X}$ and $\mathcal{U}$ are the respective feasible state and control constraint sets defined by the inequality constraints~\eqref{ineq-constraint}. Then, $x_k$ and $u_k$ are applied to the vehicle model, and the next state $x_{k+1}$ of the system is recorded. This constructs the required transition $\{x_{k}, u_{k}, x_{k+1}\}$ for training. The input to the neural network is the state and control pair $(x_{k}, u_{k})$, and the corresponding output is incremental state change $\Delta x_k=x_{k+1} - x_{k}$. The input and output are also normalized to ensure an almost equal contribution of each element of the state and control to the loss function. The obtained training dataset is $\mathcal{D}$ containing a million data points.

\subsubsection{Neural Network Structure and Training}\label{sec:training-nn}

We use a dense neural network with four hidden layers of size $200, 300, 300, 100$, respectively, with the rectified linear function as the activation function. We then train the neural network by minimizing the mean squared error
\begin{align}
\mathcal{L}(\theta; \mathcal{D}) = \frac{1}{M} \sum_{i=1}^M  \norm{\Delta  x_{k}^{(i)} - \hat{f}(x_{k}^{(i)}, u_{k}^{(i)}; \theta)}^{2},
\label{eq:neural_net_loss}
\end{align}
where $M$ is the number of training samples in  $\mathcal{D}$. The training is done using Adam optimizer~\cite{Kingma:ICLR:2015}. We also use a different test dataset $\mathcal{D}_{\mathrm{test}}$ to evaluate the prediction accuracy of the trained neural network over data points not included in the training set. Finally, we point out that such a dense network is used in order to better validate the CAP-NMPC approach, even though a less dense one could also sufficiently approximate the vehicle dynamics.

\begin{figure}[t]
	  \centering
	  \includegraphics[trim=25 80 35 100 ,clip,width=\linewidth]{straight_traj.pdf}
	  \caption{Planned trajectory for   Scenario 1 using the CAP-NMPC approach.}
	  \label{fig:straight-track-static-traj}
\end{figure}

\begin{figure}[t]
	  \centering
	  \includegraphics[trim = 25 15 30 15,clip,width=\linewidth]{Straight_Track_static_control.pdf}
	  \caption{Velocity, acceleration and steering profiles for Scenario 1.}
	  \label{fig:straight-track-static-control}
\end{figure}
\section{Numerical Simulation}\label{numerical-sim}

Based on Section~\ref{motion-planning-formulation}, we now illustrate the effectiveness of the CAP-NMPC approach for motion planning.

\subsection{Constraints}\label{constraint-def}

The planning vehicle, or ego car, has a state $x_k$  including its 2D position $(x_k^p,y_k^p)$, linear velocity $\nu_k$, and heading angle $\psi_k$, as in the considered single-track vehicle model~\cite{rajmani}. Its control input  $u_k$ includes the acceleration $a_k$ and steering angle $\delta_k$. Three constraints are considered in the simulation study: road boundaries, obstacle collision avoidance, and the limits on the control actions.

\subsubsection{Road Boundary Constraints}\label{road-boundaries}
The  ego car is constrained to stay within the road boundaries. The maximum orthogonal distance $d_B$ from the road boundary is
\begin{align*} \label{distance-to-track}
d_B(x_k) = \max_j \left\{ \norm{x_k-B_j}_2 \right\},
\end{align*}
where $B_j$ denotes the orthogonal point on the $j$-th boundary of the road from the ego car. We require $d_B(x_k)\leq L$, where $L$ is the width of the road, such that  the ego car lies inside the road boundaries.

\subsubsection{Obstacle Avoidance Constraint} \label{obstacle-avoidance}The ego car can encounter many other vehicles in its path and should keep a safe distance from them. The distance between the ego car and the closest obstacles is expressed as
\begin{align*}
d_O(x_k) = \min_l  \norm{x_k-O_l}_2 ,
\end{align*}
where $O_l$ is the $l$-th obstacle. To ensure safety, we enforce  $d_O(x_k) \geq 0$.

\begin{figure}[t]
	  \centering
	  \includegraphics[trim=15 5 15 15,clip,width=\linewidth]{Straight_Track_stat_dist.pdf}
	  \caption{Minimum relative distances between the ego car and static obstacles for Scenario 1.}
	  \label{fig:sine-track-static-dist}
\end{figure}
\subsubsection{Control Constraints}
The maneuverability of the car is limited in real world. Here, we represent these limits by specifying the control  bounds:
\begin{equation*}\label{control-bounds}
\underline{u} \leq u_k \leq \overline{u},
\end{equation*}
where $\underline{u}$ is the lower bound, and $\overline{u}$ is the upper bound.
\begin{figure}[t]
	  \centering
	  \includegraphics[trim = 20 70 40 55 ,clip,width=\linewidth]{sine_traj_ann.pdf}
	  \caption{Planned trajectory for   Scenario 2 using the CAP-NMPC approach.}
	  \label{fig:sine-track-moving-traj}
\end{figure}
\begin{figure}[t]
	  \centering
	  \includegraphics[trim=25 15 30 15,clip,width=\linewidth]{Sine_Track_moving_control.pdf}
	  \caption{Velocity, acceleration and steering profiles for Scenario 2.}
	  \label{fig:sine-track-moving-control}
\end{figure}
\begin{figure}[t]
	  \centering
	  \includegraphics[trim=15 5 15 15,clip,width=\linewidth]{Sine_Track_moving_dist.pdf}
	  \caption{Minimum relative distances between the ego car and other moving vehicles in Scenario 2.}
	  \label{fig:sine-track-moving-dist}
\end{figure}

\subsection{Case Studies}
We consider two scenarios to examine the performance of the CAP-NMPC approach for   motion planning. In Scenario 1, the ego car is to travel on a straight road with static  obstacles. Scenario 2 involves  three moving vehicles on  a curved road and requires the ego car to overtake them. The prediction horizon length in the scenarios is $H=10$. The number of particles used in the CAP-NMPC is $N=300$. All the cars' lengths and widths are $4\: \mathrm{m}$ and $2\: \mathrm{m}$, respectively. The road width is $10\: \mathrm{m}$, and the sampling period $\Delta{t}=0.2\: \mathrm{s}$.

{\em Scenario 1: Straight road with static obstacles.}
We place three stationary obstacles across the road. The ego car's trajectory is shown in Fig.~\ref{fig:straight-track-static-traj}, where the depth of the color from light to dark encodes the temporal information from present to future. The solid black line represents the road boundaries. We can observe that the ego vehicle is able to smoothly maneuver, bypassing the obstacles and keeping a safe distance at all times. Fig.~\ref{fig:straight-track-static-control} depicts the vehicles velocity and control profiles. It is evident that the CAP-NMPC allows the ego car to plan proficient maneuvers from zero velocity until reaching the goal state while satisfying the specified constraints and reference velocity (pink dotted line in Fig.~\ref{fig:straight-track-static-control}). It can be seen that the generated path with CAP-NMPC is indeed feasible for the vehicle to execute. In addition, in Fig.~\ref{fig:sine-track-static-dist}, we depict the minimum relative distance between the ego car and the obstacles, which is kept above zero to ensure collision avoidance.

{\em Scenario 2: Curved road with other moving vehicles.}
In this scenario, we consider a curved road with moving obstacles. The initial speed of all vehicles is $5 \:\mathrm{m/s}$ and the trajectories of the obstacle vehicles are obtained by assigning a goal position for each obstacle vehicle by CAP-NMPC algorithm. The obstacle vehicles are set in motion only when the ego car approaches them in order to ensure interaction with the ego car. This scenario is significantly more challenging given the dynamic environment, making collision avoidance a difficult task. For instance, observe from Figs. 4-6 the instant where the ego car is at $x = 20\:\mathrm{m}$. At this instant, the ego car has two other cars adjacent to it, one in the left lane and the other straight ahead. In addition, the road boundary is on the right of the ego car. The ego car has to carefully maneuver without leaving the road and colliding with the adjacent cars. To this end, the CAP-NMPC could successfully predict the upcoming situation and increase its velocity to smoothly pass between the two vehicles at a safe distance, as shown in Fig.~\ref{fig:sine-track-moving-dist}. Eventually, the ego car formulates a trajectory allowing it to overtake the other vehicles without collision to reach the specified goal state.

With this, we demonstrated the effectiveness of using a sample-based NMPC framework to solve the motion planning problem of a vehicle while considering neural network dynamics through the use of our proposed CAP-NMPC algorithm for realistic driving scenarios.

\section{Conclusion}

A growing  convergence of machine learning and advanced control is driving the frontiers of robotics and  stimulates new research problems. In this paper, we investigated the    problem of NMPC for neural network dynamics. The motivation lies in that, even though  NMPC has become an important method for robotics applications, its popular implementation based on numerical optimization can   meet only limited success when given neural network models. In a departure, we proposed to use a sampling-based NMPC approach, which built upon an Bayesian estimation perspective of NMPC and leveraged particle filtering/smoothing to estimate the best control actions. We then considered the motion planning problem for an autonomous vehicle with neural network dynamics and deployed the proposed approach to solve it under different driving conditions. The results demonstrated the potency of CAP-NMPC to handle constraints and compute   feasible trajectories effectively.
Our future work will pursue integrating   the CAP-NMPC approach with other existing, very successful sampling-based motion planning techniques such as rapidly-exploring random trees or probabilistic roadmap methods, toward further improving the computational efficiency and feasibility in planning. The proposed approach can find prospective use in a wide variety of other robot control problems.

\balance
\bibliographystyle{IEEEtran}