# Connecting low-energy CP violation, resonant leptogenesis and neutrinoless
double beta decay in a radiative seesaw model

Bikash Thapa and Ng. K. Francis bikash2<EMAIL_ADDRESS>Department of Physics, Tezpur University, Tezpur - 784028

###### Abstract

We present a study of resonant leptogenesis in a radiative seesaw model. We
consider the case where two quasi-degenerate right-handed neutrinos realize
resonant leptogenesis, and the CP violation necessary to achieve leptogenesis
occurs through the CP phases present in the neutrino mixing matrix. A
numerical analysis is performed by taking the best-fit values from the current
global data for three neutrino mixing angles and two mass-squared differences.
We have shown how the predicted value of baryon asymmetry depends on the Dirac
and Majorana CP phases. With the particular choice for the mass parameters,
this model prefers a normal hierarchy of neutrino masses based on the value of
baryon asymmetry predicted. Using the constrained CP phases, we evaluate the
effective neutrino mass, which is relevant to the neutrinoless double beta
decay.

##  1 Introduction

Evidence from neutrino oscillation experiments [1] suggests that the standard
model (SM) is incomplete and cannot incorporate neutrino mass. An attractive
possibility of physics beyond the SM to account for tiny neutrino mass is the
seesaw mechanism [2, 3, 4]. In addition to explaining the tiny neutrino mass,
the baryon asymmetry of the Universe (BAU) can be dynamically generated within
the seesaw mechanism via leptogenesis [5]. In leptogenesis, the CP-violating
and out-of-equilibrium decays of right-handed neutrinos generate lepton
asymmetry, which is then partially converted into baryon asymmetry by the
electroweak sphaleron processes [6].

In standard thermal leptogenesis with hierarchical masses of the right-handed
neutrinos, the observed value of BAU can be explained if their mass scale is
$\mathcal{O}$($10^{9}$) GeV [7]. However, this mass scale may be lower in
cases with nearly degenerate right-handed neutrinos. Such a scenario is known
as Resonant leptogenesis [8, 9]. In resonant leptogenesis scenario, sufficient
amount of BAU may be achieved by assuming the mass degeneracy of the right-
handed neutrinos to be $\leq 10^{-8}$. Such a strong degeneracy is unnatural
and fine-tuned.

Also, since resonant leptogenesis occurs at lower temperatures, the CP phases
present in the neutrino mixing matrix can act as a source of CP violation
required to explain the observed BAU via leptogenesis successfully. The value
of CP phases in the neutrino mass matrix is not as established as the mixing
angles by various neutrino oscillation experiments. Measurements from long-
baseline experiments such as T2K [10] and NO$\nu$A [11] along with reactor
experiments: Daya-Bay [12], RENO [13] and Double-Chooz [14] suggests a
preference for Dirac phase, $\delta\sim 1.5\pi$. Future experiments like DUNE
[15] and T2HK [16] may give a much more precise determination of the Dirac CP-
violating phase. On the other hand, if neutrinos are of Majorana nature, the
standard parameterization of the Pontecorvo-Maki-Naka-Sakata (PMNS) [17, 18,
19] mixing matrix will have two additional Majorana phases. The Majorana
nature of massive neutrinos may be probed in neutrinoless double beta decay
experiments, and these experiments can also provide us with knowledge about
the Majorana phases. These motivates us to consider the investigation of
baryogenesis and explore its dependence on CP phases present in the PMNS
matrix.

Even though the Yukawa coupling matrix, crucial to the generation of CP
asymmetry through the decays of right-handed neutrinos, cannot be exclusively
reconstructed by low-energy neutrino observables. However, bridging the
connection between low-energy CP phases and leptogenesis is interesting.
Numerous studies explore this connection in different frameworks [20, 21, 22,
23, 24, 25, 26, 27]. The authors of reference [20] show the connection between
low-energy CP violation and leptogenesis in the minimal seesaw model with two
right-handed neutrinos. They explored the scenario of baryogenesis through
resonant leptogenesis, considering all the CP violation requirement for
successful leptogenesis is fulfilled by the CP phases present in the neutrino
mixing matrix. Our approach is to explore the parameter space for the low-
energy CP phases such that a successful explanation of the observed BAU may be
given within the framework of the radiative seesaw model with an inert Higgs
doublet.

In this work, we consider a model of radiative neutrino masses having two
right-handed neutrinos with TeV-scale masses. First, we have analyzed the
scenario of resonant leptogenesis by considering a tiny splitting between the
two right-handed neutrinos. An interesting feature of this model is that the
observed BAU can be achieved even if the degeneracy is relaxed [28], making
the model much more natural. Then, we carry out the study of resonant
leptogenesis, assuming that the CP violation required comes exclusively from
the neutrino mixing matrix. In other words, we will show the parameter space
for the Dirac and Majorana phases such that the observed value of BAU
constrains them.

This paper is organized as follows: In section 2, we have presented a brief
description of the framework of the scotogenic model with two right-handed
neutrinos. We have written the Yukawa coupling matrix in Casas-Ibarra type
parameterization and established assumptions that the CP violation comes
exclusively from the neutrino mixing matrix. In section 3, we have introduced
the terms relevant for evaluating baryon asymmetry. We explore the parameter
space of our model using the constraints from the observed BAU in this
section. We further show the results of our numerical evaluation of the
coupled Boltzmann equation, which govern the evolution of the number density
of various particles involved. Section 4 includes the results on neutrinoless
double beta decay obtained using the constrained parameter space of section 3.
We finally conclude our work in section 5.

##  2 Radiative neutrino mass model with inert an Higgs doublet

In this work, we have considered the study of a radiative seesaw model, first
proposed by E. Ma [29]. The model shows an extension of the fermion sector
with two right-handed neutrinos, $N_{i}$, and the scalar sector includes an
additional inert Higgs doublet $\eta$. The particle content of the model under
$SU(2)_{L}\times U(1)_{Y}\times Z_{2}$ may be summarised as

$\displaystyle(\nu^{\alpha}_{L},\alpha_{L})\sim(2$
$\displaystyle,-\frac{1}{2},+),~{}~{}\alpha^{\alpha}_{R}\sim(1,1,+),~{}~{}(\Phi^{+},\Phi^{0})\sim(2,\frac{1}{2},+)$
$\displaystyle
N_{i}\sim(1,0,-),~{}~{}(\eta^{+},\eta^{0})\sim(2,\frac{1}{2},-)$ (1)

The SM fields are $Z_{2}$ even whereas the new fields are assumed to be odd
under the $Z_{2}$ symmetry. The invariant Lagrangian of the new fields may be
written as

$\mathcal{L}\supset-h_{\alpha
i}\bar{l}^{\alpha}_{L}\tilde{\eta}N_{i}+\frac{1}{2}M_{i}\bar{N}_{i}(N^{c})_{i}+h.c.,$
(2)

where $h_{\alpha i}$ denote the Yukawa couplings, $l^{\alpha}_{L}$ (with
$\alpha=e,\mu,\tau$) is the SM lepton doublets,
$\tilde{\eta}=i\sigma_{2}\eta^{*}$, with $\sigma_{2}$ being the 2nd Pauli
matrices and $M_{i}$ is the Majorana masses of the right-handed neutrinos. It
is assumed that the inert Higgs doublet has no vacuum expectation value
($vev$), and we find that no Dirac mass term is generated after electroweak
symmetry breaking (EWSB). The neutrino remains massless at the tree level and
can be generated at the one-loop level. The potential of the model becomes

$\displaystyle
V=-\mu^{2}H^{\dagger}H+m^{2}_{\eta}\eta^{\dagger}\eta+\frac{\lambda_{1}}{2}(H^{\dagger}H)^{2}+\frac{\lambda_{2}}{2}(\eta^{\dagger}\eta)^{2}+\lambda_{3}(H^{\dagger}H)$
$\displaystyle(\eta^{\dagger}\eta)+\lambda_{4}(H^{\dagger}\eta)(\eta^{\dagger}H)$
$\displaystyle\frac{\lambda_{5}}{2}\left[(H^{\dagger}\eta)(H^{\dagger}\eta)+(\eta^{\dagger}H)(\eta^{\dagger}H)\right],$
(3)

where $\lambda_{i}$ represents the quartic couplings and are assumed to be
real without loss of generality. If we take $H=(0,(v+h)/\sqrt{2})^{T}$ and
$\eta=(\eta^{+},(\eta_{R}+i\eta_{I})/\sqrt{2})^{T}$, the masses of the
physical scalar states can be written as

$\displaystyle m^{2}_{h}$ $\displaystyle=\lambda_{1}v^{2},$ $\displaystyle
m^{2}_{\eta^{\pm}}$ $\displaystyle=m_{\eta}^{2}+\frac{v^{2}}{2}\lambda_{3},$
$\displaystyle m^{2}_{\eta_{R}}$
$\displaystyle=m_{\eta}^{2}+\frac{v^{2}}{2}(\lambda_{3}+\lambda_{4}+\lambda_{5}),$
$\displaystyle m^{2}_{\eta_{I}}$
$\displaystyle=m_{\eta}^{2}+\frac{v^{2}}{2}(\lambda_{3}+\lambda_{4}-\lambda_{5}),$
(4)

where $v$ is the $vev$ of the SM Higgs doublet, $m_{h}$ denotes the mass of
the SM Higgs, $m_{\eta_{R}}$ and $m_{\eta_{I}}$ are the mass of the real
scalar $\eta_{R}$ and the real pseudoscalar $\eta_{I}$, respectively. It is
assumed that both $\eta_{R}$ and $\eta_{I}$ are lighter than the complex
scalar $\eta^{\pm}$, whose mass is represented by $m^{2}_{\eta^{\pm}}$ in Eq.
(2). It is clear that the mass difference between $\eta_{R}$ and $\eta_{I}$ is
$m_{\eta_{R}}-m_{\eta_{I}}=v^{2}\lambda_{5}$ and therefore, as
$\lambda_{5}\rightarrow 0$, $\eta_{R}$ and $\eta_{I}$ becomes degenerate. We
note that in this model, $Z_{2}$ symmetry stabilizes the lightest $Z_{2}$ odd
particle, and this particle can play the role of dark matter candidate111Here,
we have restricted our study to baryon asymmetry of the Universe. A systematic
study of the dark matter problem in such a model is carried out in [28, 30,
31, 32].

The radiative neutrino mass model gives rise to a neutrino mass matrix of the
form[29]

$\left(M_{\nu}\right)_{\alpha\beta}=\sum_{i}\frac{M_{i}h^{*}_{\alpha
i}h^{*}_{\beta
i}}{32\pi^{2}}\left[L\left(m_{\eta_{R}}^{2}\right)-L\left(m_{\eta_{I}}^{2}\right)\right],$
(5)

with $L\left(m^{2}\right)$ defined as,

$L\left(m^{2}\right)=\frac{m^{2}}{m^{2}-M_{i}^{2}}\log\left(\frac{m^{2}}{M_{i}^{2}}\right).$
(6)

Analogous to the type-I seesaw formula we write the neutrino mass matrix of
Eq. (5) as,

$M_{\nu}=h^{*}\Lambda^{-1}h^{\dagger},$ (7)

where $\Lambda$ is a diagonal matrix of the form [30],

$\Lambda_{i}=\frac{2\pi^{2}}{\lambda_{5}}\xi_{i}\frac{2M_{i}}{v^{2}}$ (8)

and

$\xi_{i}=\left(\frac{1}{8}\frac{M_{i}^{2}}{m_{\eta_{R}}^{2}-m_{\eta_{I}}^{2}}\left[L\left(m_{\eta_{R}}^{2}\right)-L\left(m_{\eta_{I}}^{2}\right)\right]\right)^{-1}.$
(9)

The active neutrino mass matrix given in Eq. (5) can be diagonalized as
$U^{\dagger}M_{\nu}U^{*}=D_{\nu}=\textrm{diag}(m_{1},m_{2},m_{3})$, where the
$m_{i}$ are the masses of the light neutrino eigenstate. Here, $U$ is a
unitary matrix also known as PMNS matrix and is represented as

$U=\begin{pmatrix}c_{12}c_{13}&s_{12}c_{13}&s_{13}e^{-i\delta_{CP}}\\\
-c_{23}s_{12}-s_{23}c_{12}s_{13}e^{i\delta_{CP}}&c_{23}c_{12}-s_{23}s_{12}s_{13}e^{i\delta_{CP}}&s_{23}c_{13}\\\
s_{23}s_{12}-c_{23}c_{12}s_{13}e^{i\delta_{CP}}&-s_{23}c_{12}-c_{23}s_{12}s_{13}e^{i\delta_{CP}}&c_{23}c_{13}\end{pmatrix}\cdot\textrm{diag}(1,e^{i\sigma},e^{i\rho}),$
(10)

where $s_{ij}=\sin\theta_{ij}$ and $c_{ij}=\cos\theta_{ij}$, $\delta_{CP}$ is
the Dirac phase and $\sigma$, $\rho$ are the Majorana phases. Since, we have
extended the fermion sector with two right-handed neutrinos, only two light
neutrinos achieve non-zero mass eigenvalues i.e., either $m_{1}=0$ (Normal
Hierarchy) or $m_{3}=0$ (Inverted Hierarchy). This also means that one of the
Majorana phases in $U$ becomes not well defined and only one combination of
Majorana phases remains, namely, ($\sigma-\rho$) in NH or $\sigma$ in IH.

In order to study the scenario of resonant leptogenesis in this model it is
convenient to write the Yukawa matrix $h$ in Casas-Ibarra (CI) type
parametrization [33],

$h=U\sqrt{D_{\nu}}R^{\dagger}\sqrt{\Lambda}$ (11)

where $R$ is a complex orthogonal matrix satisfying $RR^{T}=1$, which is of
the form

$R=\begin{pmatrix}\cos\omega_{12}&\sin\omega_{12}&0\\\
-\sin\omega_{12}&\cos\omega_{12}&0\\\
0&0&1\end{pmatrix}\begin{pmatrix}\cos\omega_{13}&0&\sin\omega_{13}\\\ 0&1&0\\\
-\sin\omega_{13}&0&\cos\omega_{13}\end{pmatrix}\begin{pmatrix}1&0&0\\\
0&\cos\omega_{23}&\sin\omega_{23}\\\
0&-\sin\omega_{23}&\cos\omega_{23}\end{pmatrix},$ (12)

where $\omega_{ij}$ are complex parameters. For two right-handed neutrino
case, we have $\omega_{12}=\omega_{13}=0$ in NH and
$\omega_{13}=\omega_{23}=0$ in IH case. Also, in this work we primarily focus
on the possibility that the CP phases in the PMNS matrix of Eq. (10) acts as a
source of the BAU. Under such a scenario we are considering Im$\omega_{23}$ or
Im$\omega_{12}$ to be zero for the NH or IH case, respectively i.e., the CP
violating parameters responsible for successful baryogenesis are $\delta_{CP}$
and $(\sigma-\rho)$ or $\sigma$.

##  3 Resonant leptogenesis

In this section, we discuss the scenario of resonant leptogenesis within the
radiative neutrino mass model with an inert Higgs doublet, popular as the
scotogenic model. In thermal leptogenesis, where the mass spectrum of the
right-handed neutrino is hierarchical, the lower bound on the mass of the
right-handed neutrino is $\mathcal{O}$($10^{9}$) GeV [7]. This high mass scale
may be lowered in the case of resonant leptogenesis, where the masses of the
right-handed neutrinos are nearly degenerate. In the type-I seesaw model with
nearly degenerate right-handed neutrinos, the observed BAU can be explained
through CP violating, out-of-equilibrium decays of right-handed neutrinos with
TeV-scale masses. However, this comes at the cost of strong, fine-tuned
degeneracy. Such a strong degeneracy is relaxed ($>10^{-8}$) in this model of
a radiative seesaw with an inert Higgs doublet.

The interference of the tree-level decay of right-handed neutrino with the
one-loop self-energy and vertex diagrams gives CP violation and hence produces
non-zero lepton asymmetry. In resonant leptogenesis, the self-energy
correction is resonantly enhanced, and the flavor-dependent asymmetry
parameter is given by [34, 35]

$\varepsilon_{\alpha i}=\sum_{j\neq
i}\frac{\textrm{Im}\left[h^{\dagger}_{i\alpha}h_{\alpha
j}\left(h^{\dagger}h\right)_{ij}\right]+\frac{M_{i}}{M_{j}}\textrm{Im}\left[h^{\dagger}_{i\alpha}h_{\alpha
j}\left(h^{\dagger}h\right)_{ji}\right]}{\left(h^{\dagger}h\right)_{ii}\left(h^{\dagger}h\right)_{jj}}\cdot\frac{\left(M_{i}^{2}-M_{j}^{2}\right)M_{i}\Gamma_{j}}{\left(M_{i}^{2}-M_{j}^{2}\right)^{2}+M_{i}^{2}\Gamma_{j}^{2}},$
(13)

where the decay width $\Gamma_{i}$ is defined as

$\Gamma_{i}=\frac{M_{i}}{8\pi}\left(h^{\dagger}h\right)_{ii}\left(1-\eta_{i}\right)^{2},~{}\textrm{with}~{}\eta_{i}=\frac{m_{\eta}^{2}}{M_{i}^{2}}$
(14)

The CP asymmetry generated from the decay of nearly degenerate right-handed
neutrinos generates the lepton asymmetry, whose value can be evaluated by
analysis of the Boltzmann equation. Here, we solve the coupled Boltzmann
equations (Eq. (15)) which describes the evolution of the right-handed
neutrino density, N${}_{N_{i}}$ and lepton number density,
N${}_{N_{\alpha\alpha}}$ ($\alpha=e,\mu,\tau$) [8, 36].

$\displaystyle\frac{d\textrm{N}_{N_{i}}}{dz}$ $\displaystyle=-D_{i}\
\left(\textrm{N}_{N_{i}}-\textrm{N}_{N_{i}}^{\textrm{eq}}\right)$ (15)
$\displaystyle\frac{d\textrm{N}_{\alpha\alpha}}{dz}$
$\displaystyle=-\sum_{i=1}^{2}\varepsilon_{i\alpha}D_{i}\left(\textrm{N}_{N_{i}}-\textrm{N}_{N_{i}}^{\textrm{eq}}\right)$
$\displaystyle\qquad-\frac{1}{4}\left\\{\sum_{i=1}^{2}(rz)^{2}D_{i}\mathcal{K}_{2}(rz)+W_{\Delta
L=2}\right\\}\textrm{N}_{\alpha\alpha},$ (16)

The above Boltzmann equations are flavor-diagonal, and we have considered the
decays, inverse decays of the right-handed neutrinos, and the $\Delta L=2$
washout processes. The equilibrium number density of the right-handed
neutrinos is given by,

$\textrm{N}_{N_{i}}^{\textrm{eq}}=\frac{45}{2\pi^{2}g_{*}}z^{2}\mathcal{K}_{2}(z),$
(17)

where $\mathcal{K}_{2}(z)$ and $g_{*}$ are the modified Bessel function of the
second kind and the number of relativistic degrees of freedom, respectively.
The decay parameter, $D_{i}$ is defined as

$D_{i}=\frac{z}{H(z=1)}\frac{\Gamma_{N_{i}}}{\textrm{N}_{N_{i}}^{\textrm{eq}}},$
(18)

with $H$ being the Hubble parameter. The $\Delta L=2$ washout processes
consist of scattering such as $l\eta\leftrightarrow\bar{l}\eta^{*}$ and
$ll\leftrightarrow\eta^{*}\eta^{*}$ and can be written as [30]

$\Delta W=\frac{\Gamma_{\Delta
L=2}}{Hz}=\frac{36\sqrt{5}M_{Pl}}{\pi^{\frac{1}{2}}g_{l}\sqrt{g_{*}}v^{4}}\frac{1}{z^{2}}\frac{1}{\lambda_{5}^{2}}M\bar{m}_{\xi}^{2}$
(19)

where the $\bar{m}_{\xi}$ is the effective mass parameter and is defined as

$\displaystyle\bar{m}_{\xi}^{2}$
$\displaystyle=\sum_{i,j}\xi_{i}\xi_{j}\textrm{Re}\left[\left(RD_{\nu}R^{\dagger}\right)_{ij}^{2}\right]$
$\displaystyle\approx
4\xi_{1}^{2}m_{1}^{2}+\xi_{2}^{2}m_{2}^{2}+\xi_{3}^{2}m_{3}^{2}$ (20)

To examine the possibility of CP phases present in the PMNS matrix as the sole
source of CP violation necessary to achieve successful leptogenesis, we take
Re$\omega_{23}$ or Re$\omega_{12}$ to be $\frac{\pi}{4}$ for NH or IH cases,
respectively. We fix, the mass of right-handed neutrinos, $M_{1}=5$ TeV, the
mass splitting parameter, $\Delta M=10^{-7}$, $m_{\eta_{R}}=1.5$ TeV, and
$\lambda_{5}=10^{-3}$. We observe that with these choices of parameters, in
the case of IH, the value of baryon asymmetry, $\eta_{B}$ is quite small,
$\mathcal{O}$($10^{-20}$), and hence, we will confine our study to the NH
case. Since the BAU is generated at the TeV scale, we take into consideration
the flavor effects of leptogenesis. This makes the CP violation in the PMNS
matrix relevant to the production of the BAU. Here, the primary focus is to
show how the generated baryon asymmetry depends on the CP phases $\delta_{CP}$
and $\alpha=\sigma-\rho$.

To estimate the baryon asymmetry, $\eta_{B}$, we numerically solve the
Boltzmann equations in Eqns. (15) and (3). We scan the $\delta_{CP}-\alpha$
plane by taking the whole range of $[0,2\pi]$ for both $\delta_{CP}$ as well
as $\alpha$ and evaluate $\eta_{B}$ by taking the best-fit values for the
three mixing angles and two mass-squared differences from the global results
of latest neutrino oscillation experiments [37]. To compare the numerical
results for $\eta_{B}$ with the experimental results given by Planck [38], we
make a $\chi^{2}$ analysis and minimize the function,

$\chi^{2}=\sum_{i}\frac{\left(\lambda^{model}-\lambda^{expt}\right)^{2}}{\Delta\lambda^{2}},$
(21)

where $\lambda^{model}$ represents the $\eta_{B}$ predicted by the model,
$\lambda^{expt}$ is the experimentally measured value of $\eta_{B}$, and
$\Delta\lambda$ denote the $1\sigma$ range.

Figure 1: The left panel shows the contour plot of $\eta_{B}$ in the
$\delta_{CP}-\alpha$ plane. The red cross mark indicates the
$\chi^{2}$-minimum value for ($\delta_{CP},\alpha$). $\eta_{B}$ has negative
values in regions with dotted lines. The right panel shows the allowed region
of parameter space constrained by the Planck bound, $\eta_{B}=(6.10\pm
0.04)\times 10^{-10}$.

Figure 1 shows the results of the parameter scan. The contour plot of
$\eta_{B}$ in the $\delta_{CP}-\alpha$ plane is shown in the left panel. It is
clear that the value of $\eta_{B}$ calculated in this scenario depends on both
the CP phases. The right panel of figure 1 represents the region in
$\delta_{CP}-\alpha$ space that is constrained by the observed value of BAU,
$\eta_{B}=(6.10\pm 0.04)\times 10^{-10}$. The best-fit value in the parameter
space of our model is evaluated using Eq. (21), which is denoted by a red star
mark and corresponds to $\delta_{CP}=0.46\pi$ and $\alpha=0.62\pi$. Using the
best-fit point, we show in figure 2 the evolution of lepton number density for
three flavors and baryon asymmetry as a function of $z$. The final value of
baryon asymmetry is observed at a high value of $z$ and is found to be
$\eta_{B}\approx 6.1\times 10^{-10}$.

Figure 2: Variation of lepton number density Nαα and $\eta_{B}$ as a function
of $z$.

In this section, we demonstrate how BAU could be explained via resonant
leptogenesis within the scotogenic model such that the required CP violation
comes exclusively from the PMNS matrix. By confronting the predictions of our
model with experimentally measured $\eta_{B}$, we obtained the allowed region
for the CP phases. In the next section, we will study the process of
neutrinoless double beta decay. We evaluate the effective light neutrino mass
on which the amplitude of neutrinoless double beta decay depends.

##  4 Neutrinoless double beta decay

If the neutrinos are of Majorana nature it can mediate the lepton number
violating neutrinoless double beta decay: $(A,Z)\rightarrow(A,Z+2)+2e^{-}$
[39, 40]. The decay rate of such a process is proportional to the effective
light neutrino mass, $|\langle m_{ee}\rangle|$ and is given by 222The
contribution from heavy right-handed leptons is negligible, and hence we take
into consideration the effect of only the light neutrinos.

$|\langle m_{ee}\rangle|=\lvert\sum_{i}m_{i}U_{ei}^{2}\rvert.$ (22)

We take the allowed region of parameter space obtained through the requirement
of successful resonant leptogenesis and show how the predicted $|\langle
m_{ee}\rangle|$ is constrained. The left panel of figure 3 shows the
calculated value of effective light neutrino mass, $|\langle m_{ee}\rangle|$
by scanning $\delta_{CP}$ and $\alpha$ over the range $[0,2\pi]$. It is clear
that the value of $|\langle m_{ee}\rangle|$ depends significantly on
$\delta_{CP}$ as well. In the right panel of figure 3, we show the predicted
values of $|\langle m_{ee}\rangle|$ constrained by the successful generation
of BAU via resonant leptogenesis. The predicted values range from 22.9 meV to
32.8 meV. The minimal value of $|\langle m_{ee}\rangle|$ is obtained when
$\alpha=0.8\pi$ and the maximal value is observed at $\alpha=1.7\pi$.

Figure 3: The left panel shows the contour plot of $|\langle m_{ee}\rangle|$
in $\delta_{CP}-\alpha$ plane. The red star mark indicate the best-fit point
corresponding to $\chi^{2}$-minimum value. The right panel shows the variation
of effective light neutrino mass with the Majorana phase, $\alpha$,
constrained by the requirement of successful resonant leptogenesis.

##  5 Conclusions

In this work, we studied a radiative seesaw model of neutrino masses with an
additional inert Higgs doublet. Such a model is fascinating as it can explain
the origin of neutrino masses, BAU, and Dark Matter. Here, we have considered
the extension of the fermion sector by two nearly degenerate right-handed
neutrinos with TeV-scale masses. Considering the quasi-degenerate nature of
the right-handed neutrinos, we investigate the possibility of generating the
BAU through the scenario of resonant leptogenesis. In this case leptogenesis
occurs at TeV-scale temperature through the out-of-equilibrium, CP violating
decay of $Z_{2}$ odd right-handed neutrino into $Z_{2}$ even SM leptons and
inert Higgs doublet: $N_{i}\rightarrow l\eta$,
$N_{i}\rightarrow\bar{l}\eta^{*}$.

We have considered that the BAU can be generated with the CP violation, which
may be measured at neutrino oscillation and neutrinoless double beta decay
experiments. In other words, we have assumed that the CP violation necessary
for successful leptogenesis comes from the phases present in the neutrino
mixing matrix (PMNS matrix). Furthermore, we found that the model prefers NH
of neutrino masses for the choices of mass parameters and the quartic coupling
presented in this work. Thus, the parameters upon which the estimation of BAU
depends are the Dirac phase, $\delta_{CP}$, and the Majorana phase,
$\alpha=\sigma-\rho$.

We numerically solved the coupled Boltzmann equations and found that this
model can explain the observed BAU even when the source of CP violation are
the CP phases, $\delta_{CP}$, and $\alpha$. Thus, we found the region of
parameter space that successfully produces the observed baryon asymmetry. The
best-fit point values for the parameters corresponding to $\chi^{2}$-minimum
value are $\delta_{CP}=0.46\pi$ and $\alpha=0.62\pi$. Next, we discussed how
the requirement of successful resonant leptogenesis constrains the effective
light neutrino mass, which is relevant for neutrinoless double beta decay.
Using the constrained parameter space of the model, we calculated the
$\lvert\langle m_{ee}\rangle\rvert$. We found that significantly, the
$\lvert\langle m_{ee}\rangle\rvert$ depends on both the CP phases. Future
generation neutrinoless double beta decay experiments can probe such a model.

## Acknowledgement

BT acknowledges the Department of Science and Technology (DST), Government of
India for INSPIRE Fellowship vide Grant No. DST/INSPIRE/2018/IF180588. The
research of NKF is funded by DST-SERB, India under Grant No. EMR/2015/001683.

## References

  * [1] Yoshiyuki Fukuda, T Hayakawa, E Ichihara, K Inoue, K Ishihara, Hirokazu Ishino, Y Itow, T Kajita, J Kameda, S Kasuga, et al.  Evidence for oscillation of atmospheric neutrinos.  Physical review letters, 81(8):1562, 1998.
  * [2] Peter Minkowski.  $\mu$→ e$\gamma$ at a rate of one out of 109 muon decays?  Physics Letters B, 67(4):421–428, 1977.
  * [3] Tsutomu Yanagida.  Horizontal symmetry and masses of neutrinos.  Progress of Theoretical Physics, 64(3):1103–1105, 1980.
  * [4] Rabindra N Mohapatra.  Mechanism for understanding small neutrino mass in superstring theories.  Physical review letters, 56(6):561, 1986.
  * [5] Masataka Fukugita and Tsutomu Yanagida.  Barygenesis without grand unification.  Physics Letters B, 174(1):45–47, 1986.
  * [6] Vadim A Kuzmin, Valery A Rubakov, and Mikhail E Shaposhnikov.  On anomalous electroweak baryon-number non-conservation in the early universe.  Physics Letters B, 155(1-2):36–42, 1985.
  * [7] Sacha Davidson and Alejandro Ibarra.  A lower bound on the right-handed neutrino mass from leptogenesis.  Physics Letters B, 535(1-4):25–32, 2002.
  * [8] Apostolos Pilaftsis and Thomas EJ Underwood.  Resonant leptogenesis.  Nuclear Physics B, 692(3):303–345, 2004.
  * [9] Apostolos Pilaftsis and Thomas EJ Underwood.  Electroweak-scale resonant leptogenesis.  Physical Review D, 72(11):113001, 2005.
  * [10] K Abe, N Abgrall, H Aihara, Y Ajima, JB Albert, D Allan, P-A Amaudruz, C Andreopoulos, B Andrieu, MD Anerella, et al.  The T2K experiment.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 659(1):106–135, 2011.
  * [11] DS Ayres, NO$\nu$A Collaboration, et al.  Proposal to build a 30 kiloton off-axis detector to study $\nu$$\mu$→ $\nu$e oscillations in the numi beamline.  arXiv preprint hep-ex/0503053, 30, 2005.
  * [12] FP An, JZ Bai, AB Balantekin, HR Band, D Beavis, W Beriguete, M Bishai, S Blyth, K Boddy, RL Brown, et al.  Observation of electron-antineutrino disappearance at Daya bay.  Physical Review Letters, 108(17):171803, 2012.
  * [13] Jung Keun Ahn, S Chebotaryov, JH Choi, S Choi, W Choi, Y Choi, HI Jang, JS Jang, EJ Jeon, IS Jeong, et al.  Observation of reactor electron antineutrinos disappearance in the RENO experiment.  Physical Review Letters, 108(19):191802, 2012.
  * [14] Y Abe, JC Dos Anjos, JC Barriere, E Baussan, I Bekman, M Bergevin, TJC Bezerra, L Bezrukov, E Blucher, C Buck, et al.  Improved measurements of the neutrino mixing angle $\theta_{13}$ with the Double Chooz detector.  Journal of High Energy Physics, 2014(10):1–44, 2014.
  * [15] Dune Collaboration et al.  Long-baseline neutrino facility (LBNF) and deep underground neutrino experiment (DUNE) conceptual design report volume 1: The LBNF and DUNE projects.  arxiv. org, 2016.
  * [16] Hyper-Kamiokande Proto-Collaboration, K Abe, H Aihara, C Andreopoulos, I Anghel, A Ariga, T Ariga, R Asfandiyarov, M Askins, JJ Back, et al.  Physics potential of a long-baseline neutrino oscillation experiment using a J-PARC neutrino beam and Hyper-Kamiokande.  Progress of theoretical and experimental physics, 2015(5):053C02, 2015.
  * [17] Bruno Pontecorvo.  Inverse $beta$ processes and nonconservation of lepton charge.  Zhur. Eksptl’. i Teoret. Fiz., 34, 1958.
  * [18] Ziro Maki, Masami Nakagawa, and Shoichi Sakata.  Remarks on the unified model of elementary particles.  Progress of Theoretical Physics, 28(5):870–880, 1962.
  * [19] R. L. Workman and Others.  Review of Particle Physics.  PTEP, 2022:083C01, 2022.
  * [20] Takehiko Asaka and Takahiro Yoshida.  Resonant leptogenesis at TeV-scale and neutrinoless double beta decay.  Journal of High Energy Physics, 2019(9):1–16, 2019.
  * [21] Xinyi Zhang, Jiang-Hao Yu, and Bo-Qiang Ma.  Leptogenesis from low-energy CP violation in minimal left-right symmetric model.  Nuclear Physics B, 976:115670, 2022.
  * [22] Zhi-zhong Xing and Di Zhang.  Bridging resonant leptogenesis and low-energy CP violation with an RGE-modified seesaw relation.  Physics Letters B, 804:135397, 2020.
  * [23] Kristian Moffat, S Pascoli, ST Petcov, and J Turner.  Leptogenesis from low energy CP violation.  Journal of High Energy Physics, 2019(3):1–53, 2019.
  * [24] Silvia Pascoli, Serguey T Petcov, and Antonio Riotto.  Connecting low energy leptonic CP violation to leptogenesis.  Physical Review D, 75(8):083511, 2007.
  * [25] Sin Kyu Kang.  Low-energy CP violation and leptogenesis in a minimal seesaw model.  Journal of the Korean Physical Society, 78(9):743–749, 2021.
  * [26] Gustavo C Branco, Takuya Morozumi, BM Nobre, and MN Rebelo.  A bridge between CP violation at low energies and leptogenesis.  Nuclear Physics B, 617(1-3):475–492, 2001.
  * [27] Gustavo Castello Branco, R Gonzalez Felipe, FR Joaquim, I Masina, MN Rebelo, and Carlos A Savoy.  Minimal scenarios for leptogenesis and CP violation.  Physical Review D, 67(7):073025, 2003.
  * [28] Shoichi Kashiwase and Daijiro Suematsu.  Baryon number asymmetry and dark matter in the neutrino mass model with an inert doublet.  Physical Review D, 86(5):053001, 2012.
  * [29] Ernest Ma.  Verifiable radiative seesaw mechanism of neutrino mass and dark matter.  Physical Review D, 73(7):077301, 2006.
  * [30] Thomas Hugle, Moritz Platscher, and Kai Schmitz.  Low-scale leptogenesis in the scotogenic neutrino mass model.  Physical Review D, 98(2):023020, 2018.
  * [31] Laura Lopez Honorez, Emmanuel Nezri, Josep F Oliver, and Michel HG Tytgat.  The inert doublet model: an archetype for dark matter.  Journal of Cosmology and Astroparticle Physics, 2007(02):028, 2007\.
  * [32] Ethan M Dolle and Shufang Su.  Inert dark matter.  Physical Review D, 80(5):055012, 2009.
  * [33] JA Casas and Alejandro Ibarra.  Oscillating neutrinos and $\mu$→ e, $\gamma$.  Nuclear Physics B, 618(1-2):171–204, 2001.
  * [34] Apostolos Pilaftsis.  CP violation and baryogenesis due to heavy majorana neutrinos.  Physical Review D, 56(9):5431, 1997.
  * [35] A Anisimov, A Broncano, and M Plümacher.  The CP-asymmetry in resonant leptogenesis.  Nuclear Physics B, 737(1-2):176–189, 2006.
  * [36] Andrea De Simone and Antonio Riotto.  Quantum boltzmann equations and leptogenesis.  Journal of Cosmology and Astroparticle Physics, 2007(08):002, 2007\.
  * [37] Pablo F de Salas, DV Forero, S Gariazzo, P Martínez-Miravé, O Mena, CA Ternes, M Tórtola, and JWF Valle.  2020 global reassessment of the neutrino oscillation picture.  Journal of High Energy Physics, 2021(2):1–36, 2021.
  * [38] Peter AR Ade, N Aghanim, M Arnaud, Mark Ashdown, J Aumont, C Baccigalupi, AJ Banday, RB Barreiro, JG Bartlett, N Bartolo, et al.  Planck 2015 results-xiii. cosmological parameters.  Astronomy & Astrophysics, 594:A13, 2016.
  * [39] JD Vergados, H Ejiri, and F Šimkovic.  Theory of neutrinoless double-beta decay.  Reports on Progress in Physics, 75(10):106301, 2012.
  * [40] Heinrich Päs and Werner Rodejohann.  Neutrinoless double beta decay.  New Journal of Physics, 17(11):115010, 2015.


# Bounding Optimality Gaps for Non-Convex Optimization Problems:
Applications to Nonlinear Safety-Critical Systems

Prithvi Akella and Aaron D. Ames1 *This work was supported by the AFOSR Test
and Evaluation Program, grant FA9550-19-1-03021All authors are with the
California Institute of Technology {pakella<EMAIL_ADDRESS>
###### Abstract

Efficient methods to provide sub-optimal solutions to non-convex optimization
problems with knowledge of the solution’s sub-optimality would facilitate the
widespread application of nonlinear optimal control algorithms. To that end,
leveraging recent work in risk-aware verification, we provide two algorithms
to (1) probabilistically bound the optimality gaps of solutions reported by
novel percentile optimization techniques, and (2) probabilistically bound the
maximum optimality gap reported by percentile approaches for repetitive
applications, e.g. Model Predictive Control (MPC). Notably, our results work
for a large class of optimization problems. We showcase the efficacy and
repeatability of our results on a few, benchmark non-convex optimization
problems and the utility of our results for controls in a Nonlinear MPC
setting.

##  I INTRODUCTION

Optimal controllers have emerged in recent years as the preeminent choice for
controller synthesis as they offer a natural way of expressing (perhaps)
disparate control objectives [1, 2, 3], e.g. as in Model Predictive Control
(MPC) [4, 5, 6], control-barrier-function based quadratic programs [7, 8, 9],
and optimal path planning [10, 11, 12], among others. Specifying to the first
two examples, optimal controllers in this vein leverage existing system models
to inform controller choice at the current or successive time-steps, e.g. the
predictive aspect of MPC references the utilization of a system model to
predict how a given input will affect system evolution and uses that
prediction to determine optimal input choice for the next few time-steps.
Notably, one-step control-barrier-function-based quadratic programs ensure
convexity for nonlinear system models, though predicting over these models
remains a non-convex problem [13, 14, 15, 16].

The utility of these controllers inspired efforts to develop and exploit rapid
solution techniques for these non-convex optimization problems, though some
improvement avenues still exist. For example, typical approaches in the
literature include Model Predictive Path Integral Control (MPPI) [17, 18],
Collocation methods [19, 20], Bayesian Optimization [21, 22], Learning-based
methods [23, 24], and more recently, percentile-based methods [25, 26]. In
specific for the latter three learning-based methods, both bayesian
optimization and other learning-based methods guarantee eventual convergence
to an optimal solution but do not provide sample complexity bounds. Percentile
methods offer robust sample complexity bounds and have seen considerable
success at efficiently providing sub-optimal solutions to constrained
nonlinear optimal controllers. However, unlike the first two methods, they
fail to provide even probabilistic bounds on the optimality gap between the
reported percentile solution and the true optimal value.

Figure 1: An overview of our proposed method. Sampling methods employed in
risk-aware verification provide decent results to non-convex optimization
problems, though their optimality gaps are not well-understood. By defining a
variance function around the previously generated sample set, we can
implicitly identify decisions (orange on the bottom) whose variances exceed
the optimality gap of our reported sampling solution and identify such an
optimality gap via a second application of the same sampling methods.

Our Contribution: We aim to provide probabilistic guarantees on the optimality
gap of percentile solutions — the deviance in cost evaluation of the
percentile solution and the minimum function value — for their respective
optimization problems. In doing so, we hope to develop the first in a line of
algorithms that provide robust sample complexity bounds for efficient
identification of sub-optimal solutions with knowledge of the solution’s sub-
optimality. To that end, our contributions are twofold:

  * •

Subject to a fairness assumption, we provide sample-complexity bounds for the
identification of the optimality gap for percentile solutions to a large class
of optimization problems.

  * •

For optimization problems to be solved repeatedly, e.g. MPC, we provide a
method to determine the probability of percentile solutions exhibiting
optimality gaps lower than a determined upper bound.

We develop such a procedure to offer efficient solutions to constrained
optimal control problems with known sample complexity and sub-optimality
bounds. This procedure would deviate from existing approaches insofar as it
would consider general (perhaps) nonlinear constraints for optimal
controllers, and provide sample-complexity results on the sub-optimality of
the provided solution — a discrepancy with existing methods, e.g. MPPI,
Collocation, Bayesian Optimization, and other, learning-based approaches.

Structure: Section II reviews percentile optimization. Section III-A outlines
our solution to bound the optimality gaps of percentile solutions. Section
III-B extends these results for optimization problems to be solved repeatedly,
e.g. Nonlinear MPC. Finally, Section IV provides examples to validate our
theoretical results. Specifically, Section IV-A validates the probabilistic
bounding results of Section III-A, and Section IV-C illustrates the utility of
the results in Section III-B in a Nonlinear MPC setting.

##  II A Review of Percentile Optimization

To facilitate the statement of our contributions, we first review percentile
optimization in this section. Consider the following optimization problem:

$J^{*}=\min_{\operatorname{s}\in\operatorname{\mathbb{S}}}~{}J(\operatorname{s}),$
(1)

subject to the following assumption:

######  Assumption 1.

The decision space $\operatorname{\mathbb{S}}$ is a set with bounded volume,
i.e.
$\int_{\operatorname{\mathbb{S}}}~{}1~{}ds=V_{\operatorname{\mathbb{S}}}<\infty$
or $\operatorname{\mathbb{S}}$ has a finite number of elements. Furthermore,
the cost function $J$ is bounded over $\operatorname{\mathbb{S}}$, i.e.
$\exists~{}m,M\in\mathbb{R},\operatorname*{~{}\mathrm{s.t.}~{}}m\leq
J(\operatorname{s})\leq
M,~{}\forall~{}\operatorname{s}\in\operatorname{\mathbb{S}}$.

This assumption permits us to define two functions, $\mathcal{V}$, which
corresponds to the volume fraction occupied by a subset $A$ of
$\operatorname{\mathbb{S}}$, and $F$, which corresponds to the set of strictly
better decisions for a provided decision
$\operatorname{s}^{\prime}\in\operatorname{\mathbb{S}}$ (we assume continuous
spaces for notational ease):

$\displaystyle\mathcal{V}(A)=\frac{\int_{A}~{}1~{}ds}{\int_{\operatorname{\mathbb{S}}}~{}1~{}ds},$
(2) $\displaystyle
F(\operatorname{s}^{\prime})=\\{\operatorname{s}\in\operatorname{\mathbb{S}}~{}|~{}J(\operatorname{s})<J(\operatorname{s}^{\prime})\\}.$
(3)

Naturally then, for a given decision
$\operatorname{s}^{\prime}\in\operatorname{\mathbb{S}}$, were
$\mathcal{V}(F(\operatorname{s}^{\prime}))\leq\epsilon$ for some
$\epsilon\in(0,1]$, i.e. $\operatorname{s}^{\prime}$ is such that the volume
fraction of strictly better decisions is no more than $\epsilon$, then
$\operatorname{s}^{\prime}$ would be in the $100(1-\epsilon)\%$-ile with
respect to minimizing $J$. Likewise, the associated minimum cost of such a
decision $J(\operatorname{s}^{\prime})$ should also be a probabilistic lower
bound on achievable costs. Both of these notions are expressed formally in the
theorem below, which combines similar results from [25, 27]. To facilitate the
statement of results to follow, we will first state the assumption, then the
theorem.

######  Assumption 2.

Let $I_{1}=\\{(\operatorname{s}_{i},J(\operatorname{s}_{i}))\\}_{i=1}^{N_{p}}$
be a set of $N_{p}$ decisions and costs for decisions $\operatorname{s}_{i}$
sampled independently via
$\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]$, with
$\zeta^{*}_{N_{p}}$ the minimum sampled cost and
$\operatorname{s}^{*}_{N_{p}}$ the (perhaps) non-unique decision with minimum
cost.

######  Theorem 1.

Let Assumptions 1 and 2 hold. Then $\forall~{}\epsilon\in[0,1]$, the
probability of sampling a decision whose cost is at-least $\zeta^{*}_{N_{p}}$
is at minimum $1-\epsilon$ with confidence $1-(1-\epsilon)^{N_{p}}$, i.e.

$\operatorname{\mathbb{P}}^{N_{p}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\operatorname{\mathbb{P}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[J(\operatorname{s})\geq\zeta^{*}_{N_{p}}\right]\geq
1-\epsilon\right]\geq 1-(1-\epsilon)^{N_{p}}.$ (4)

Furthermore, $\forall~{}\epsilon\in(0,1]$, $\operatorname{s}^{*}_{N_{p}}$ is
in the $100(1-\epsilon)\%$-ile with minimum confidence
$1-(1-\epsilon)^{N_{p}}$, i.e.

$\operatorname{\mathbb{P}}^{N_{p}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\mathcal{V}(F(\operatorname{s}^{*}_{N_{p}}))\leq\epsilon\right]\geq
1-(1-\epsilon)^{N_{p}}.$ (5)

##  III Bounding the Optimality Gap

To formally bound the optimality gap, we must first define it. To that end,
consider the same general optimization problem as per (1). The optimality gap
is defined as follows.

######  Definition 1.

For general optimization problems of the form in (1), the optimality gap of a
decision $\operatorname{s}\in\operatorname{\mathbb{S}}$, denoted as
$\operatorname{G}(\operatorname{s})$, is the deviance between the decision and
optimal values, i.e.
$\operatorname{G}(\operatorname{s})=J(\operatorname{s})-J^{*}$.

Then our formal problem statement follows.

######  Problem 1.

Let Assumption 2 hold, and let the optimality gap $\operatorname{G}$ be as per
Definition 1. Identify an upper bound to
$\operatorname{G}(\operatorname{s}^{*}_{N_{p}})$ and the probability with
which this bound holds.

###  III-A Results for Bounded Optimization Problems

First, we note that if Assumption 2 holds, we can define the following
variance function $\operatorname{\mathbb{V}}$ over the set $I_{1}$:

$\operatorname{\mathbb{V}}:\operatorname{\mathbb{S}}\to\mathbb{R}\operatorname*{~{}\mathrm{s.t.}~{}}\operatorname{\mathbb{V}}(\operatorname{s})=\min_{\operatorname{s}_{i}\in
D\subseteq I_{1}}~{}|J(\operatorname{s})-J(\operatorname{s}_{i})|.$ (6)

The restriction of $\operatorname{s}\in D\subseteq I_{1}$ in the definition of
$\operatorname{\mathbb{V}}$ above is purely for practical purposes. From a
theoretical standpoint, one could use the entire information set $I_{1}$ to
define $\mathbb{V}$, but this tends to increase sample requirements as will be
discussed in sections to follow (e.g. Figure 3)

Intuitively then, we aim to maximize $\operatorname{\mathbb{V}}$ via a
percentile method to identify a variance that supersedes the optimality gap
$\operatorname{G}(\operatorname{s}^{*}_{N_{p}})$ of our chosen decision. To do
so, we first require the following fairness assumption — that it is possible
to sample variances at least as large as the optimality gap, as otherwise, it
would be impossible to take a percentile approach. To formally state this
assumption, we first define $\Omega_{r}$ to be the $r$-level set of
$\operatorname{\mathbb{V}}$:

$\Omega_{r}=\\{\operatorname{s}\in\operatorname{\mathbb{S}}~{}|~{}\operatorname{\mathbb{V}}(s)\leq
r\\}.$ (7)

######  Assumption 3.

Let the variance function $\operatorname{\mathbb{V}}$ be as per (6), let the
optimality gap $\operatorname{G}$ be as per Definition 1, let $\mathcal{V}$ be
as per (2), let $\Omega_{r}$ be as per (7), and let Assumption 2 hold. The
level set $\Omega_{\operatorname{G}(s^{*}_{N_{p}})}$ of decisions whose
variance is at most the optimality gap of $\operatorname{s}^{*}_{N_{p}}$ does
not encompass $\operatorname{\mathbb{S}}$, i.e.

$\mathcal{V}\left(\operatorname{\mathbb{S}}\setminus\Omega_{\operatorname{G}(\operatorname{s}^{*}_{N_{p}})}\right)>0.$
(8)

Second, we have the following result regarding level sets $\Omega$.

######  Lemma 1.

Let $\Omega_{r}$ be as defined in (7) and let $\mathcal{V}$ be as per (2). The
following statements are all equivalent.

  1. (a)

$r\geq s$,

  2. (b)

$\Omega_{r}\supseteq\Omega_{s}$,

  3. (c)

$\mathcal{V}(\Omega_{r})\geq\mathcal{V}(\Omega_{s})$,

  4. (d)

$\operatorname{\mathbb{P}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}[\Omega_{r}]\geq\operatorname{\mathbb{P}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}[\Omega_{s}]$.

Proof: The equivalency between (a) and (b) stems via the definition of
$\Omega_{r}$ in (7). The equivalency between (b) and (c) stems via the
definition of the volume fraction function in (2). Finally, the equivalency
between (c) and (d) stems via the uniform distribution assigning probabilistic
weight to subsets of $A\subseteq\operatorname{\mathbb{S}}$ equivalent to
$\mathcal{V}(A)$. $\mathbin{{\rule{5.38193pt}{5.38193pt}}}$

Third, based on our fairness assumption, we know there exists a non-zero
probability of sampling decisions such that their variances are at least the
optimality gap.

######  Lemma 2.

Let Assumption 3 hold, then there exists $p>0$ corresponding to the
probability of uniformly sampling a decision
$\operatorname{s}\in\operatorname{\mathbb{S}}\setminus\Omega_{\operatorname{G}(\operatorname{s}^{*}_{N_{p}})}$,
i.e.

$\operatorname{\mathbb{P}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\operatorname{\mathbb{S}}\setminus\Omega_{\operatorname{G}(\operatorname{s}^{*}_{N_{p}})}\right]=p>0.$
(9)

Proof: The result holds by definition of the uniform distribution over
$\operatorname{\mathbb{S}}$ and equation (8).
$\mathbin{{\rule{5.38193pt}{5.38193pt}}}$

Then the main result in the utilization of a percentile approach to upper
bound the optimality gap stems from the prior lemmas and assumption. Formally,
we aim to take a percentile approach to the following optimization problem:

$\max_{\operatorname{s}\in\operatorname{\mathbb{S}}}~{}\operatorname{\mathbb{V}}(\operatorname{s}),$
(10)

which results in the following theorem.

######  Theorem 2.

Let Assumptions 1 and 3 hold, let $p$ satisfy (9), let
$I_{2}=\\{(\operatorname{s}_{i},\operatorname{\mathbb{V}}(\operatorname{s}_{i}))\\}_{i=1}^{N_{v}}$
be a set of $N_{v}$ decisions $\operatorname{s}_{i}$ independently sampled
from $\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]$ with their
corresponding variances as per (6) and with
$\operatorname{\mathbb{V}}^{*}_{N_{v}}$ the maximum sampled variance. Then,
$\forall\epsilon\in[0,p]$, $\operatorname{\mathbb{V}}^{*}_{N_{v}}$ exceeds the
optimality gap of the percentile solution with confidence
$1-(1-\epsilon)^{N_{v}}$, i.e.

$\operatorname{\mathbb{P}}^{N_{v}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\operatorname{\mathbb{V}}^{*}_{N_{v}}\geq\operatorname{G}(\operatorname{s}^{*}_{N_{p}})\right]\geq
1-(1-\epsilon)^{N_{v}}.$ (11)

Proof: First, we know there exists a $p>0$ satisfying (9) via Lemma 2. Second,
we know that for the optimization problem (10), the decision space
$\operatorname{\mathbb{S}}$ and objective function $\operatorname{\mathbb{V}}$
are bounded — this stems via Assumption 1. This permits us to take a
percentile approach to solving (10). Via Theorem 1 and looking at (10) as a
minimization, we have the following result substituting terms for (4) and for
all $\epsilon\in[0,1]$:

$\operatorname{\mathbb{P}}^{N_{v}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\operatorname{\mathbb{P}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\operatorname{\mathbb{V}}(\operatorname{s})\leq\operatorname{\mathbb{V}}^{*}_{N_{v}}\right]\geq
1-\epsilon\right]\geq 1-(1-\epsilon)^{N_{v}}.$ (12)

The set in the interior probability corresponds to the level set of
$\operatorname{\mathbb{V}}^{*}_{N_{v}}$, i.e.

$\operatorname{\mathbb{P}}^{N_{v}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\operatorname{\mathbb{P}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\Omega_{\operatorname{\mathbb{V}}^{*}_{N_{v}}}\right]\geq
1-\epsilon\right]\geq 1-(1-\epsilon)^{N_{v}}.$ (13)

Finally, via (9), we know that the probability of sampling a decision in the
level set corresponding to the optimality gap is $1-p$. Restricting to
$\epsilon\in[0,p]$ which implies that $1-\epsilon\geq 1-p$ and substituting
terms in the inequality above, we have the following result.

$\operatorname{\mathbb{P}}^{N_{v}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\operatorname{\mathbb{P}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\Omega_{\operatorname{\mathbb{V}}^{*}_{N_{v}}}\right]\geq\operatorname{\mathbb{P}}_{\operatorname{\mathrm{U}}[\operatorname{\mathbb{S}}]}\left[\Omega_{\operatorname{G}(\operatorname{s}^{*}_{N_{p}})}\right]\right]\geq
1-(1-\epsilon)^{N_{v}}.$ (14)

Then, the final result holds due to Lemma 1.
$\mathbin{{\rule{5.38193pt}{5.38193pt}}}$

In summary, Theorem 2 tells us that if we wish to bound the optimality gap of
a percentile solution, we need to evaluate the variance of $N_{v}$ uniform
samples $\operatorname{s}$ from $\operatorname{\mathbb{S}}$ with respect to a
subset of the information set $I_{1}$ utilized to generate the percentile
solution. In practice, however, the exact probability $p$ of sampling
decisions with large enough variance will be unknown to the practitioner
apriori. In these cases, it suffices to assume a small enough value for
$\epsilon$, i.e. $10^{-2}$ or smaller, is smaller than $p$. Examples along
this vein will be provided in the following section. Notably, this result
implies that we can utilize percentile methods to both identify decisions that
outperform a large fraction of the decision space and also determine their
optimality gap. Indeed, this result holds even for non-convex optimization
problems, provided they satisfy Assumption 1.

###  III-B Producing Solutions with Maximum Optimality Gaps

The prior section provides a method to determine the upper bound on the
optimality gap of a provided solution via a secondary sampling scheme. This
section provides a method to remove the secondary sampling requirement for
similar optimization problems to be solved successively. In other words,
consider the following general form of optimization problems, where each
instance $l$ satisfies Assumption 1:

$J^{l*}=\min_{\operatorname{s}\in\operatorname{\mathbb{S}}^{l}}~{}J^{l}(\operatorname{s}),~{}(\operatorname{\mathbb{S}}^{l},J^{l})\in\operatorname{\mathbb{O}},~{}l\in
L.$ (15)

Furthermore, we assume that it is possible to randomly sample indices $l$ from
$L$, e.g. via the uniform distribution.

Example Setting: Here, $\operatorname{\mathbb{O}}$ is a set containing pairs
of objective functions and decision spaces. To provide an example consistent
with the sections to follow, consider the following nonlinear dynamical system
with state $x$ and input $u$:

$x_{k+1}=f(x_{k},u_{k}),~{}x\in\mathcal{X},~{}u\in\mathcal{U}$ (16)

Provided a cost function over states and inputs, state constraints, and torque
bounds, we can construct the following finite-time optimal controller with
horizon $H$ for the aforementioned system. Here, all state constraints are
projected to input constraints through prediction over the model (16):

$\displaystyle\operatorname*{argmin}_{\mathbf{u}=(u^{0},u^{1},\dots,u^{H-1})\in\mathcal{U}^{H}}~{}$
$\displaystyle J(\mathbf{u},x_{k}),$ (FTOCP)
$\displaystyle~{}\qquad\mathrm{subject~{}to~{}}$
$\displaystyle\mathbf{u}\in\mathbb{U}(x_{k})\subseteq\mathcal{U}^{H}.$ (17)

The aforementioned finite-time optimal controller (FTOCP) collapses to the
form in (15) if we consider an optimality set $\operatorname{\mathbb{O}}$
indexed by states $x\in\mathcal{X}$ — a specific form of indexing more
generally referred to via $l\in L$ in (15).

Key Insight: The critical insight for this section then is as follows. If we
were to randomly sample via a distribution $\pi$ over
$\operatorname{\mathbb{O}}$, optimization problems of the form in (15) and
calculate the optimality gap $\operatorname{G}^{l}$ of percentile solutions
for that problem, the corresponding gap is a sample of some real-valued random
variable. By taking multiple independent samples of this random variable, we
can provide, using recent results on risk-measure estimation, a probabilistic
upper bound on this random variable, i.e. a probabilistic upper bound on
achievable optimality gaps. Phrased formally, an existing theorem in the
literature we will utilize is as follows.

######  Theorem 3.

Let $X$ be a real-valued random variable with (perhaps) unknown distribution
$\pi$, let $\\{x_{i}\\}_{i=1}^{N}$ be a set of $N$ independent samples of $X$,
and let $x^{*}_{N}$ be the maximum such sample. $\forall~{}\epsilon\in[0,1]$,
the following statement holds:

$\operatorname{\mathbb{P}}^{N}_{\pi}\left[\operatorname{\mathbb{P}}_{\pi}[x\leq
x^{*}_{N}]\geq 1-\epsilon\right]\geq 1-(1-\epsilon)^{N}.$ (18)

Proof: This follows from a scenario argument, see [28], and is a re-phrasing
of Theorem 4 in [25]. $\mathbin{{\rule{5.38193pt}{5.38193pt}}}$

To employ Theorem 3, we require the following definition.

######  Definition 2.

Let $\operatorname{\mathcal{G}}(N_{p})$ be a real-valued random variable with
distribution $\pi_{\operatorname{G}}$ and samples $\operatorname{\mathbf{g}}$
defined as follows: (1) Uniformly sample an index $l\in L$, (2) Take a
percentile approach to solve (15) corresponding to this sampled index $l$,
producing the solution $\operatorname{s}^{l*}_{N_{p}}$, (3) Calculate and
report as a sample $\operatorname{\mathbf{g}}$, the optimality gap
$\operatorname{G}^{l}\left(\operatorname{s}^{l*}_{N_{p}}\right)$.

Then, we can upper bound the optimality gaps of percentile solutions to all
optimization problems formed in the set $\operatorname{\mathbb{O}}$, to some
minimum probability. The formal statement will follow.

######  Theorem 4.

Let $\operatorname{\mathcal{G}}(N_{p})$ be as per Definition 2, and let
$\\{\operatorname{\mathbf{g}}_{i}\\}_{i=1}^{R}$ be a set of $R$ independent
samples of $\operatorname{\mathcal{G}}(N_{p})$ with
$\operatorname{\mathbf{g}}^{*}_{R}$ the maximum such sample. Then
$\forall~{}\epsilon\in[0,1]$, percentile solutions to optimization problems of
the form in (15) will exhibit optimality gaps less than
$\operatorname{\mathbf{g}}^{*}_{R}$ with minimum probability $1-\epsilon$ and
confidence $1-(1-\epsilon)^{R}$, i.e.

$\operatorname{\mathbb{P}}^{R}_{\pi_{\operatorname{G}}}\left[\operatorname{\mathbb{P}}_{\pi_{\operatorname{G}}}[\operatorname{\mathbf{g}}\leq\operatorname{\mathbf{g}}^{*}_{R}]\geq
1-\epsilon\right]\geq 1-(1-\epsilon)^{R}.$ (19)

Proof: Stems via Theorem 3. $\mathbin{{\rule{5.38193pt}{5.38193pt}}}$

In short, if we wish to remove the secondary sampling requirement for the
determination of optimality gaps, we need to be able to calculate the
optimality gap for at least $R$ independently chosen optimization problems of
the form (15). Doing so permits us to make a general statement on the maximum
achievable optimality gaps, to some minimum probability. Now, we will provide
a few examples.

##  IV Examples

Figure 2: Validation Data for Section IV-A corresponding to Theorem 2. (Top)
$100$ reported upper bounds $\operatorname{\mathbb{V}}^{*}_{N_{v}}$ using
Theorem 2 with desired confidence for (11) equal to $0.7$. (Middle) $100$
reported upper bounds $\operatorname{\mathbb{V}}^{*}_{N_{v}}$ with confidence
$0.999$. (Bottom) Running fraction over $2000$ trials of reported upper bounds
$\operatorname{\mathbb{V}}^{*}_{N_{v}}$ exceeding the true optimality gap
$\operatorname{G}(p^{*}_{N_{p}})$ at confidence level $0.999$. Notice how the
fraction of upper bounds exceeding the optimality gap increases as we increase
confidence (top to middle), and the running fraction of upper bounds exceeding
the optimality gap converges to our desired confidence (bottom), corroborating
Theorem 2.

###  IV-A Validating Theorem 2 — Traveling Salesman

The Traveling Salesman Problem (TSP) is a classic non-convex path planning
problem referencing the identification of the path of the shortest distance
traversing each node in a set once. Mathematically, consider a set of
waypoints $W=\\{w_{1},w_{2},\dots,w_{|W|}\\}~{}w_{i}\in\mathbb{R}^{2}$ and the
set of all paths over these waypoints
$P=\\{(i_{1},i_{2},\dots,i_{|W|})~{}|~{}i_{j}\in\\{1,2,\dots,|W|\\},~{}i_{j}\neq
i_{k},~{}\forall~{}j\neq k\\}$. Then formally, the Traveling Salesman Problem
is to

$\min_{p\in P}~{}\sum_{i=0}^{|W|-1}~{}\|p_{i}-p_{i+1}\|+\|p_{0}-p_{|W|}\|.$
(20)

For a graph with $10$ nodes and $3628800$ possible paths, evaluating
$N_{p}=5000$ paths and picking the best one identifies a path $p^{*}_{N_{p}}$
in the $99.9\%$-ile with at least $99\%$ confidence according to Theorem 1.
Using a subset $D$ of the corresponding information set $I_{1}$ for the
determination of such a percentile solution (see Assumption 2 for the
definition of $I_{1}$), we define a variance function
$\operatorname{\mathbb{V}}$ as per equation (6). Finally, to validate the
probabilistic results of Theorem 2, we can also calculate the true probability
$p$ of uniformly sampling paths that exhibit a variance higher than the
optimality gap of our proposed solution $\operatorname{G}(p^{*}_{N_{p}})$ —
this is the minimum probability assumed to exist via Assumption 3 and defined
in equation (9). For this particular node set and percentile solution
$p^{*}_{N_{p}}$ the probability $p=0.1083$.

Name |  $N_{p}$ |  $N_{v}$ |  expected |  |
---|---|---|---|---|---
success |  |  |  |  |
probability (11) |  true success |  |  |  |
probability |  average |  |  |  |
runtime (ms) |  |  |  |  |
R-2 |  $300$ |  $300$ |  $0.95$ |  $\approx 1$ |  5.18
R-10 |  $300$ |  $300$ |  $0.95$ |  $\approx 1$ |  5.17
Ack |  $300$ |  $300$ |  $0.95$ |  $\approx 1$ |  5.28
Ble |  $300$ |  $300$ |  $0.95$ |  $\approx 1$ |  5.17
Levi |  $300$ |  $300$ |  $0.95$ |  $\approx 1$ |  5.31
Himm |  $300$ |  $300$ |  $0.95$ |  $\approx 1$ |  5.30

TABLE I: Data for Section IV-B for the (R-2) Rastigrin 2-D, (R-10) Rastigrin
10-D, (Ack) Ackley, (Ble) Beale, Levi, and (Himm) Himmelblau benchmark
optimization problems.

Figure 3: Validation data for Section IV-B. We claim that by varying the
amount of information used to generate the variance function
$\operatorname{\mathbb{V}}$ in equation (6), we can change the baseline
probability $p$ of sampling a decision whose variance exceeds the optimality
gap of a given solution $\operatorname{G}(s^{*}_{N})$ (such decisions are
highlighted in orange). Notice that as the volume fractions $\chi$ occupied by
the chosen information set $D$ decreases, we see a corresponding increase in
the baseline probability $p$. Section IV-B discusses why this inverse
relationship holds.

Figure Analysis: To validate the results of Theorem 2, we solved for the
minimum number of samples $N_{v}$ required to determine an upper bound
$\operatorname{\mathbb{V}}^{*}_{N_{v}}$ to $\operatorname{G}(p^{*}_{N_{p}})$
with minimum confidence $0.7$ — top figure in Figure 2 requiring $11$ samples—
and minimum confidence $0.999$ — middle figure in Figure 2 requiring $61$
samples. Both of these minimum sample requirements were identified by setting
$\epsilon=p$ in (11) and solving for the minimum integer $N_{v}$ required to
make the right-hand side greater than or equal to our desired confidence. For
each confidence level, we performed $100$ separate trials and reported the
maximum variance $\operatorname{\mathbb{V}}^{*}_{N_{v}}$ produced by each
trial according to Theorem 2. As can be seen in the corresponding data,
increasing the confidence increases the likelihood that the corresponding
reported result $\operatorname{\mathbb{V}}^{*}_{N_{v}}$ exceeds the true
optimality gap $\operatorname{G}(p^{*}_{N_{p}})$ — the red $x$’es in the
middle figure are all above the black, dashed line, whereas a few dip below
the same line in the top figure when we report solutions with lower
confidence. Furthermore, by repeating the procedure once more at confidence
$0.999$, taking $2000$ separate verification runs, and recording whether
$\operatorname{\mathbb{V}}^{*}_{N_{v}}\geq\operatorname{G}(p^{*}_{N_{p}})$ per
run, we can get a sense of the true, running probability that
$\operatorname{\mathbb{V}}^{*}_{N_{v}}\geq\operatorname{G}(p^{*}_{N_{p}})$. As
can be seen in the bottom figure, this probability converges to $0.999$ — the
lower bound expected by Theorem 2. Notably, though, this result implies that
we were able to identify a path in the $99\%$-ile that was no more than $1.12$
times the length of the optimal path, by only evaluating $5061$ uniformly
chosen path samples, less than $0.14\%$ of the overall decision space.

###  IV-B Increasing Success Probabilities — Benchmark Functions

In a brief remark after defining the variance function in equation (6), we
mentioned that by specific choice of a subset $D$, one could increase the
baseline probability $p$ of uniformly choosing samples that exhibit a higher
variance than the optimality gap of the reported percentile solution. This
section provides evidence in support of that statement for a few benchmark
optimization problems. The one referenced in Figure 3, the 2-d Rastigrin
function [29], is as follows:

$\min_{x\in[-5.12,5.12]^{2}}~{}20+\sum_{i=1}^{2}(x_{i}^{2}-10\cos(2\pi
x_{i})).$ (21)

Both the decision space and objective function are bounded, permitting a
percentile solution to (21). Following Theorem 1 and taking $N_{p}=100$
samples for such an approach, we generate the information set $I_{1}$ and
percentile solution $x^{*}_{N_{p}}$. Then, by choice of a subset $D$ of
$I_{1}$ for the definition of the variance function in (6), we claim we can
vary the baseline probability $p$ of sampling decisions that exhibit a higher
variance than the true optimality gap. To show this, define $\chi$ to be the
volume fraction the subset $D$ in (6) occupies of $I_{1}$, i.e. $\chi=1$
implies $D=I_{1}$ and $\chi=0.05$ implies that $D$ contains $1/20$-th as many
elements as $I_{1}$.

Figure Analysis: Figure 3 portrays the results of varying the volume fraction
$\chi$ of decisions $D$ utilized to generate the variance function
$\operatorname{\mathbb{V}}$ per equation (6). The decisions highlighted in
orange are those that exhibit a higher variance than the optimality gap of the
reported solution. Notice that as $\chi$ decreases, $p$ increases as indicated
in the titles. This inverse relationship arises as if we consider two sets
$D_{1}\subset D_{2}$ utilized for generation of the variance functions
$\operatorname{\mathbb{V}}_{1},\operatorname{\mathbb{V}}_{2}$, respectively,
then $\forall~{}\operatorname{s}\in\operatorname{\mathbb{S}}$,
$\operatorname{\mathbb{V}}_{1}(s)\geq\operatorname{\mathbb{V}}_{2}(s)$ by
definition of $\operatorname{\mathbb{V}}$ as a minimization problem. In other
words, decreasing the amount of information provided to the variance function
increases the corresponding conservativeness of the resulting function, which
increases the likelihood of sampling a, now more conservative, upper bound on
the optimality gap. In practice, and in the examples provided in the prior
section, using a dilation $\chi=0.1$ proved most effective, though studying if
there exists an optimal volume fraction remains an open problem and the
subject of future work.

Table Description: By increasing the success probability $p$ of sampling
decisions whose variance exceeds the optimality gap, we can “blindly” use
Theorem 2 to bound the optimality gap of percentile solutions to benchmark
optimization problems. Table IV-A shows our data in this vein. For each
benchmark optimization problem, we produced a percentile solution using
$N_{p}=300$ samples, constructed a variance function
$\operatorname{\mathbb{V}}$ using $10\%$ of the information set $I_{1}$
generated via the percentile method, and took $N_{v}=300$ samples to identify
the probabilistic maximum variance $\operatorname{\mathbb{V}}^{*}_{N_{v}}$.
Under the assumption that $p\geq 0.01$, Theorem 1 tells us that
$\operatorname{\mathbb{V}}^{*}_{N_{v}}\geq\operatorname{G}(x^{*}_{N_{p}})$
with $95\%$ probability — column 4 in Table IV-A. Indeed, over $5000$ trials
following the above procedure for each optimization problem, we were
successfully able to identify a valid upper bound every time. Additionally, as
the sampling method only requires the evaluation of sampled points, which is
relatively quick, the procedure takes very little time to implement, as seen
in the rightmost column.

Figure 4: Validation data for Section IV-C in support of Theorem 4. We claim
that we can upper bound the optimality gap of successive applications of
percentile methods to solve optimization problems of the form in (15). Shown
above in red are the calculated upper bounds for the black lines corresponding
to the $99\%$ cutoff value of optimality gaps for percentile solutions to
(NMPC-B). For the three separate percentile methods shown, we’re able to upper
bound the true value every time, corroborating Theorem 4. This information is
further discussed in Section IV-C.

###  IV-C Validating Theorem 4 — Nonlinear MPC

To validate the results of Theorem 4, we require a series of optimization
problems of the form in (15). Keeping with the example provided in Section
III-B, we aim to bound the optimality gap of percentile solutions for a
Nonlinear MPC controller steering the Robotarium robots [30] — a collection of
agents modelable via unicycle dynamics. Formally, let
$x\in\mathcal{X}=[-1.6,1.6]\times[-1.2,1.2]\times[0,2\pi]$ be the system
state, and let $u\in\mathcal{U}=[-0.2,0.2]\times[-\pi,\pi]$ be the control
input. Then, our model is as follows:

$\displaystyle x_{k+1}$
$\displaystyle=\underbrace{x_{k}+\begin{bmatrix}\cos\left(x_{k}[3]\right)&0\\\
\sin\left(x_{k}[3]\right)&0\\\ 0&1\end{bmatrix}u_{k}(\Delta
t=0.033)}_{f(x_{k},u_{k})}.$ (22)

Furthermore, each agent has a Lyapunov controller $U$ steering it to a
provided waypoint $w\in\mathcal{W}$:

$\displaystyle
U:\mathcal{X}\times\mathcal{W}\triangleq[-1.6,1.6]\times[-1.2,1.2]\to\mathcal{U}.$
(23)

We will use $U$ to construct our MPC algorithm, which provides waypoints to
steer the system around static and moving obstacles toward at least one goal.
Formally, we overlay an $8\times 5$ grid on $\mathcal{W}$ and define the space
of operating environments $\mathcal{D}$ as those environments that: (1) have
$8$ static obstacles (SO) and $3$ goals ($g$), (2) have controlled ($x_{a}$)
and uncontrolled ($x_{o}$) agent starting positions outside static obstacles
and goals, and (3) have at least one feasible path from the controlled agent’s
starting location to one of the three goals. A vector
$d=[x_{a},x_{o},\mathrm{SO},g]\in\mathcal{D}$ corresponds to one such
environmental setup. To account for collisions, consider the following barrier
function $h$, where $P=[\mathbf{I}_{2\times 2}~{}\mathbf{0}_{2\times 1}]$
projects system states to the plane [7]:

$h(x_{a},x_{o},d)=\begin{cases}-5&\mbox{in~{}SO~{}cell},\\\
\|P(x_{a}-x_{o})\|-0.18&\mbox{else}.\end{cases}$ (24)

Then the nominal NMPC algorithm minimizes $S:\mathcal{W}\to\mathbb{R}$ — a
function outputting the shortest path distance from a waypoint to the closest
goal — while ensuring that the existing Lyapunov controller $U$ renders the
barrier $h$ positive for the next $5$ time steps, i.e.,
$\forall~{}j=1,2,\dots,5$,

$\displaystyle w^{*}_{k}=$
$\displaystyle~{}\operatorname*{argmin}_{w\in\mathcal{W}}~{}$ $\displaystyle
S(w),$ (NMPC-A) $\displaystyle\mathrm{subject~{}to~{}}$ $\displaystyle
x^{j}_{k}=f(x^{j-1}_{k},u^{j-1}),$ (a) $\displaystyle
x^{0}_{k}=x_{k},{\color[rgb]{1,1,1}\eqref{eq:constr_2},\eqref{eq:constr_3}}$
(b) $\displaystyle h(x^{j}_{k,a},x_{o},d)\geq 0$ (c) $\displaystyle
u^{j-1}=U\left(x^{j-1}_{k},w\right),$ (d) $\displaystyle
0.05\leq\|w-x_{k}\|\leq 0.2.$ (25)

To ease sampling requirements, we consider an augmented cost $J$ that outputs
$100$ whenever a waypoint $w$ fails to satisfy constraints (a)-(d) in
(NMPC-A), yielding the following:

$\displaystyle w^{d*}=$
$\displaystyle~{}\operatorname*{argmin}_{w\in\mathcal{W}}~{}$ $\displaystyle
J(w,d),$ (NMPC-B) $\displaystyle\mathrm{subject~{}to~{}}$ $\displaystyle
0.05\leq\|w-x\|\leq 0.2.$ (26)

Finally, we note that (NMPC-B) is equivalent to (15), if we consider as index
set $L$, the environment set $\mathcal{D}$.

Figure Analysis: To validate the results of Theorem 4 we uniformly sample
$R=459$ different environments $d\in\mathcal{D}$ and calculate a percentile
solution to (NMPC-B) with $N_{p}=200$ samples, Figure 4-(a); $N_{p}=300$
samples, Figure 4-(b); and $N_{p}=500$ samples, Figure 4-(c). We calculate the
optimality gap for each solution by performing gradient descent on the best
out of $2000$ uniformly chosen samples and reporting the final value as the
true optimal value. According to Theorem 4, in each case, we should produce an
upper bound on optimality gaps $\operatorname{\mathbf{g}}^{*}_{R}$ that
exceeds sample-able optimality gaps $\operatorname{\mathbf{g}}$ with at least
$99\%$ probability and $99\%$ confidence. To validate this statement, we
uniformly sampled $50000$ more environments $d\in\mathcal{D}$ and followed the
prior scheme for each percentile case to determine the distribution of sample-
able optimality gaps. This data is reflected as the distributional data you
see in each subfigure in Figure 4. As can be seen, in each case the reported
upper bound $\operatorname{\mathbf{g}}^{*}_{R}$ exceeds the true $99\%$
cutoff. Furthermore, as the number of samples taken for the percentile
solution increases, the upper bound decreases. This is expected as we are
providing a solution in a higher percentile each time. To emphasize the
utility of this result for controls, say we wished to implement a percentile
procedure with $N_{p}=300$ samples to provide “good” waypoints optimizing for
(NMPC-B). Offline calculation of this optimality gap would provide confidence
that in practice, we would, with $99\%$ probability and with $99\%$
confidence, be choosing waypoints within $0.03m$ of the optimal waypoint at
every iteration. Notably, we would not have to solve the non-convex program
and repetitively sample variances at each time step to make this statement.

##  V Conclusion

We provide a method to bound the optimality gaps of percentile methods to a
wide class of non-convex optimization problems and explore the utility of such
approaches for control in an NMPC setting. We also verify our mathematical
results on several benchmark optimization problems including the traveling
salesman problem. In future work, we hope to decrease the conservatism of our
bounds and formally prove the existence of a baseline probability for certain
classes of optimization problems, i.e. remove the requirement on our fairness
assumption preceding our theoretical statements.

## References

  * [1] F. L. Lewis, D. Vrabie, and V. L. Syrmos, _Optimal control_. John Wiley & Sons, 2012.
  * [2] R. E. Kalman _et al._ , “Contributions to the theory of optimal control,” _Bol. soc. mat. mexicana_ , vol. 5, no. 2, pp. 102–119, 1960.
  * [3] A. Locatelli and S. Sieniutycz, “Optimal control: An introduction,” _Appl. Mech. Rev._ , vol. 55, no. 3, pp. B48–B49, 2002.
  * [4] E. F. Camacho and C. B. Alba, _Model predictive control_. Springer science & business media, 2013.
  * [5] J. B. Rawlings, “Tutorial overview of model predictive control,” _IEEE control systems magazine_ , vol. 20, no. 3, pp. 38–52, 2000.
  * [6] C. E. Garcia, D. M. Prett, and M. Morari, “Model predictive control: Theory and practice—a survey,” _Automatica_ , vol. 25, no. 3, pp. 335–348, 1989\.
  * [7] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 8, pp. 3861–3876, 2016.
  * [8] X. Xu, P. Tabuada, J. W. Grizzle, and A. D. Ames, “Robustness of control barrier functions for safety critical control,” _IFAC-PapersOnLine_ , vol. 48, no. 27, pp. 54–61, 2015.
  * [9] R. Grandia, A. J. Taylor, A. D. Ames, and M. Hutter, “Multi-layered safety for legged robots via control barrier functions and model predictive control,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 8352–8358.
  * [10] P. Raja and S. Pugazhenthi, “Optimal path planning of mobile robots: A review,” _International journal of physical sciences_ , vol. 7, no. 9, pp. 1314–1320, 2012.
  * [11] I. Noreen, A. Khan, and Z. Habib, “Optimal path planning using rrt* based approaches: a survey and future directions,” _International Journal of Advanced Computer Science and Applications_ , vol. 7, no. 11, 2016.
  * [12] B. Riviere, W. Hönig, Y. Yue, and S.-J. Chung, “Glas: Global-to-local safe autonomy synthesis for multi-robot motion planning with end-to-end learning,” _IEEE robotics and automation letters_ , vol. 5, no. 3, pp. 4249–4256, 2020.
  * [13] H. Ma, J. Chen, S. Eben, Z. Lin, Y. Guan, Y. Ren, and S. Zheng, “Model-based constrained reinforcement learning using generalized control barrier function,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 4552–4559.
  * [14] A. Agrawal and K. Sreenath, “Discrete control barrier functions for safety-critical control of discrete systems with application to bipedal robot navigation.” in _Robotics: Science and Systems_ , vol. 13. Cambridge, MA, USA, 2017.
  * [15] J. Zeng, Z. Li, and K. Sreenath, “Enhancing feasibility and safety of nonlinear model predictive control with discrete-time control barrier functions,” in _2021 60th IEEE Conference on Decision and Control (CDC)_. IEEE, 2021, pp. 6137–6144.
  * [16] J. Zeng, B. Zhang, and K. Sreenath, “Safety-critical model predictive control with discrete-time control barrier function,” in _2021 American Control Conference (ACC)_. IEEE, 2021, pp. 3882–3889.
  * [17] G. Williams, A. Aldrich, and E. A. Theodorou, “Model predictive path integral control: From theory to parallel computation,” _Journal of Guidance, Control, and Dynamics_ , vol. 40, no. 2, pp. 344–357, 2017.
  * [18] M. S. Gandhi, B. Vlahov, J. Gibson, G. Williams, and E. A. Theodorou, “Robust model predictive path integral control: Analysis and performance guarantees,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 1423–1430, 2021.
  * [19] J. Tamimi and P. Li, “Nonlinear model predictive control using multiple shooting combined with collocation on finite elements,” _IFAC Proceedings Volumes_ , vol. 42, no. 11, pp. 703–708, 2009.
  * [20] ——, “A combined approach to nonlinear model predictive control of fast systems,” _Journal of Process Control_ , vol. 20, no. 9, pp. 1092–1102, 2010\.
  * [21] O. Andersson, M. Wzorek, P. Rudol, and P. Doherty, “Model-predictive control with stochastic collision avoidance using bayesian policy optimization,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2016, pp. 4597–4604.
  * [22] A. Jain and M. Morari, “Computing the racing line using bayesian optimization,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 6192–6197.
  * [23] Y. Vaupel, N. C. Hamacher, A. Caspari, A. Mhamdi, I. G. Kevrekidis, and A. Mitsos, “Accelerating nonlinear model predictive control through machine learning,” _Journal of process control_ , vol. 92, pp. 261–270, 2020.
  * [24] D. Limon, J. Calliess, and J. M. Maciejowski, “Learning-based nonlinear model predictive control,” _IFAC-PapersOnLine_ , vol. 50, no. 1, pp. 7769–7776, 2017.
  * [25] P. Akella, A. Dixit, M. Ahmadi, J. W. Burdick, and A. D. Ames, “Sample-based bounds for coherent risk measures: Applications to policy synthesis and verification,” _arXiv preprint arXiv:2204.09833_ , 2022.
  * [26] P. Akella, W. Ubellacker, and A. D. Ames, “Probabilistic guarantees for nonlinear safety-critical optimal control,” _arXiv preprint arXiv:2303.06258_ , 2023.
  * [27] P. Akella, M. Ahmadi, and A. D. Ames, “A scenario approach to risk-aware safety-critical system verification,” _arXiv preprint arXiv:2203.02595_ , 2022.
  * [28] M. C. Campi and S. Garatti, “The exact feasibility of randomized solutions of uncertain convex programs,” _SIAM Journal on Optimization_ , vol. 19, no. 3, pp. 1211–1230, 2008.
  * [29] L. A. Rastrigin, “Systems of extremal control,” _Nauka_ , 1974.
  * [30] S. Wilson, P. Glotfelter, L. Wang, S. Mayya, G. Notomista, M. Mote, and M. Egerstedt, “The robotarium: Globally impactful opportunities, challenges, and lessons learned in remote-access, distributed control of multirobot systems,” _IEEE Control Systems Magazine_ , vol. 40, no. 1, pp. 26–44, 2020\.