# Isospin-density dependent pairing from infinite nuclear matter to finite
nuclei

Xu Meng  1 School of Physics, Beihang University, Beijing 100191, China
Shisheng Zhang  1 School of Physics, Beihang University, Beijing 100191, China
Lin Guo  1 School of Physics, Beihang University, Beijing 100191, China
Lisheng Geng  mbox1 School of Physics, Beihang University, Beijing 100191,
China
2 Beijing Key Laboratory of Advanced Nuclear Materials and Physics, Beihang
University, Beijing 100191, China
3 Beijing Advanced Innovation Center for Big Data-Based Precision Medicine,
School of Medicine and Engineering, Beihang University, Beijing, 100191
4School of Physics and Microelectronics, Zhengzhou University, Zhengzhou,
Henan 450001, China
Ligang Cao  5 College of Science and Technology, Beijing Normal University,
Beijing 100875, China

###### Abstract

The effective isospin-density dependent pairing interaction (P1) [S. S. Zhang,
U. Lombardo and E. G. Zhao, Sci. Chin. Phys. Mech. Astro. 54, 236 (2011)]
extracted from neutron pairing gaps for 1S0 in asymmetric nuclear matter
calculations [S. S. Zhang, L. G. Cao, U. Lombardo, et al. Phys. Rev. C 81,
044313 (2010)] is employed to study the bulk properties of Ca, Ni, Zr and Sn
isotopes. The odd-even mass (OEM) staggering is calculated by the Skyrme
Hartree-Fock plus BCS method (SHF + BCS) with the SkP interaction. For
comparison, we study two other types of isovector effective pairing
interactions. One is also extracted from pairing gaps of infinite nuclear
matter by the Brueckner-Hartree-Fock (BHF) method but for free spectrum (P2).
The other is obtained by fitting the empirical OEM (P3). An isoscalar
effective pairing interaction (P4) is also adopted which is determined by
fitting the empirical OEM. We find that interaction P1 can better describe the
OEM staggering of Ni, Zr and Sn isotopes by $14.3\%,41\%,30.4\%$ compared with
interaction P2, in terms of root mean square deviations to the empirical OEM,
respectively. On the other hand, the performance of P1 and P2 is comparable
for Ca isotopes. For Ca and Ni isotopes, P1 behaves similarly as P3, but for
Zr isotopes P1 is better than P3 by $\sim 34\%$. Among the four pairing
interactions studied, P4 performs the worst. Therefore, one may conclude that
for neutron pairings in finite nuclei, the isovector pairings are preferred
than the isoscalar one. It is quite interesting to note that the pairing
interaction P1 extracted from nuclear matter calculations can describe pairing
gaps of finite nuclei as well as or even better than the interaction P3
directly fitted to finite nuclei.
Keywords: nuclear matter, pairing correlation, odd-even mass staggering,
finite nuclei

###### pacs:

21.10.Dr, 21.30.Fe, 21.60.-n

##  I Introduction

Pairing correlations play an important role in describing many observables and
processes in nuclear physics, for instance, the odd-even mass (OEM) staggering
for finite nuclei Zhang2011_SCPMA54-236 ; Zhang2010_PRC81-044313 , the
superfluidity and cooling of neutron stars Dean2003_RMP75-607 ;
Dong2016_TAJ817-6 , r-process nucleosynthesis Zhang2015_PRC91-045802 , etc.

Tremendous efforts have been made to extract pairing interactions from
different observables. One way is to fit the OEM staggering of finite nuclei
via the empirical three-point or five-point pairing gap formulae with the
experimental binding energies as inputs Bertulani2009_PRC80-027303 ;
Goriely2009_PRL102-152503 . But pairing gaps can be different for three-point
and five-point formulae, and sometimes cannot reproduce small pairing for
nuclei with (double) magic numbers or truly reflect the experimental OEM
difference since energy density functionals for odd-A systems are not as good
as those for even-even systems Hinohara2016_PRL116-152502 . Recently, there is
a proposal that nucleonic pairing can be extracted from nuclear density
functional theory for pairing rotational bands in even-even nuclei with the
Quasi-particle Random Phase Approximation (QRPA) method
Hinohara2016_PRL116-152502 . A separable force of finite range is widely used
to describe pairing correlations in normal nuclei Tianyuan2009_PLB676-44  and
has recently been applied to provide effective pairing interactions for
hyperons Rong2020_PLB807-135533 . Another alternative is to extract pairing
interactions from asymmetric nuclear matter (ANM) calculations with the
microscopic Brueckner-Hartree-Fock (BHF) method Zhang2011_SCPMA54-236 ;
Margueron2008_PRC77-054309 ; Chamel2008_NPA812-72 ; Zhang2010_PRC81-044313 ;
Zhang2016_PRC93-044329 , and adopt local density approximation to obtain the
isospin-density dependent parameters for finite nuclei. In the last decade,
Margueron, Sagawa and Hagino introduced zero-range isospin-density dependent
effective pairing interactions Margueron2007_PRC76-064316 ;
Margueron2008_PRC77-054309  by fitting to the corresponding pairing gaps of
symmetric nuclear matter (SNM) and pure neutron matter (PNM) obtained by the
BHF method with and without medium polarization effect Cao2006_PRC74-064301 .
But those pairing gaps presented in Ref. Margueron2007_PRC76-064316  are
extracted from free spectrum instead of the mean field spectrum as clarified
in our previous paper Zhang2011_SCPMA54-236 . In that paper, we proposed a new
effective pairing interaction Zhang2011_SCPMA54-236  from the mean field
spectrum based on self consistent calculations, denoted by P1 in the
following.

In this paper, we aim at applying the new pairing interaction P1 to describe
finite nuclei and to see whether one can provide a universal description of
pairing correlations in nuclear matter and finite nuclei. For comparison, we
also study three other pairing interactions. And we label the pairing
interaction of Ref. Margueron2008_PRC77-054309  as P2, which is obtained in
the same way as P1 but fitted to pairing gaps of free spectrum instead of the
mean field. We also consider a second isovector pairing interaction, referred
to as P3, which is extracted from fitting to the experimental OEM using the
empirical three-point formula Bertulani2009_PRC80-027303 . In addition to the
above two isovector types of pairing interactions, we also study an isoscalar
pairing interaction, denoted as P4, to check the impact of isospin dependence.

For the mean field part, we use the `ev8` code Bonche2005_CPC171-49 ;
Ryssens2015_CPC187-175  with the SkP force. The particle-particle channel is
described by the BCS approximation with four different pairing interactions
detailed above. Experimental binding energies are taken from AME2016
Wang2017_ChinPhysC41-030003 .

The paper is organized as follows. In Sec. II, we give a brief introduction of
the Skyrme Hartree-Fock plus BCS method (SHF + BCS) and describe the isospin-
density dependent and isoscalar pairing interactions. Numerical details are
also presented in this section. Then, we take Ca, Ni, Zr and Sn isotopes as
examples to compare the results of four effective pairing interactions in Sec.
III. Discussions on the corresponding potentials are analyzed as well.
Finally, we make a brief summary in Sec. IV.

##  II Theoretical Framework

In this section, we briefly review the SHF + BCS method, describe the isospin-
density dependent pairing interactions, and spell out some numerical details.

###  II.1 Skyrme force

The Skyrme force is widely used in Hartree-Fock calculations. Its energy
density functional sly4  contains eight terms

$\mathcal{H}=\mathcal{K}+\mathcal{H}_{0}+\mathcal{H}_{3}+\mathcal{H}_{\rm
eff}+\mathcal{H}_{\rm fin}+\mathcal{H}_{\rm so}+\mathcal{H}_{\rm
sg}+\mathcal{H}_{\rm Coul}.\\\ $ (1)

They are the kinetic energy term $\mathcal{K}=(\hbar^{2}/2m)\tau$, zero range
term $\mathcal{H}_{0}$, density dependent term $\mathcal{H}_{3}$, effective
mass term $\mathcal{H}_{\rm eff}$, finite range term $\mathcal{H}_{\rm fin}$,
spin-orbit term $\mathcal{H}_{\rm so}$, and $\mathcal{H}_{\rm sg}$ due to the
tensor coupling with spin and its gradient.

These terms have the explicit expressions as follows,

$\begin{split}\mathcal{H}_{0}=&\frac{1}{4}t_{0}\left[\left(2+x_{0}\right)\rho^{2}-\left(2x_{0}+1\right)\left(\rho_{\rm
p}^{2}+\rho_{\rm n}^{2}\right)\right],\\\
\mathcal{H}_{3}=&\frac{1}{24}t_{3}\rho^{\alpha}\left[\left(2+x_{3}\right)\rho^{2}-\left(2x_{3}+1\right)\left(\rho_{\rm
p}^{2}+\rho_{\rm n}^{2}\right)\right],\\\ \mathcal{H}_{\rm
eff}=&\frac{1}{8}\left[t_{1}\left(2+x_{1}\right)+t_{2}\left(2+x_{2}\right)\right]\tau\rho\\\
&+\frac{1}{8}\left[t_{2}\left(2x_{2}+1\right)-t_{1}\left(2x_{1}+1\right)\right]\left(\tau_{\rm
p}\rho_{\rm p}+\tau_{\rm n}\rho_{\rm n}\right),\\\ \mathcal{H}_{\rm
fin}=&\frac{1}{32}\left[3t_{1}\left(2+x_{1}\right)-t_{2}\left(2+x_{2}\right)\right]\left(\nabla\rho\right)^{2}\\\
&-\frac{1}{32}\left[3t_{1}\left(2x_{1}+1\right)+t_{2}\left(2x_{2}+1\right)\right]\left[\left(\nabla\rho_{\rm
p}\right)^{2}+\left(\nabla\rho_{\rm n}\right)^{2}\right],\\\ \mathcal{H}_{\rm
so}=&\frac{1}{2}W_{0}\left[\vec{J}\cdot\nabla\rho+\vec{J}_{\rm
p}\cdot\nabla\rho_{\rm p}+\vec{J}_{\rm n}\cdot\nabla\rho_{\rm n}\right],\\\
\mathcal{H}_{\rm
sg}=&-\frac{1}{16}\left(t_{1}x_{1}+t_{2}x_{2}\right)\vec{J}^{2}+\frac{1}{16}\left(t_{1}-t_{2}\right)\left[\vec{J}_{\rm
p}^{2}+\vec{J}_{\rm n}^{2}\right],\end{split}$ (2)

in which $t_{0}$, $t_{1}$, $t_{2}$, $t_{3}$, $x_{0}$, $x_{1}$, $x_{2}$,
$x_{3}$, and $W_{0}$ are fitted parameters, and $\rho=\rho_{\rm p}+\rho_{\rm
n}$, $\tau=\tau_{\rm p}+\tau_{\rm n}$ are matter and kinetic densities,
$\vec{J}=\vec{J}_{\rm p}+\vec{J}_{\rm n}$ is the spin-orbit current sly4 . We
adopt the SkP force Dobaczewski1984_NPA422-103  in our present study, which is
obtained by paying particular attention to pairing properties and accurate
description of binding energies.

###  II.2 Isospin-density dependent pairing interactions

As a linear interpolation of the particle-particle interaction between
symmetric nuclear matter and pure neutron matter, isospin-density dependent
zero-range effective interactions are derived in Refs.
Margueron2007_PRC76-064316 ; Margueron2008_PRC77-054309 , for neutrons and
protons, respectively,

$\begin{split}g_{\rm
n}(\rho,\beta)&=1-\eta_{s}(\rho/\rho_{0})^{\alpha_{s}}(1-\beta)-\eta_{n}(\rho/\rho_{0})^{\alpha_{n}}\beta,\\\
g_{\rm
p}(\rho,\beta)&=1-\eta_{s}(\rho/\rho_{0})^{\alpha_{s}}(1+\beta)+\eta_{n}(\rho/\rho_{0})^{\alpha_{n}}\beta,\\\
\end{split}$ (3)

where the four parameters $\eta_{s},\eta_{n},\alpha_{s}$ and $\alpha_{n}$ are
adjusted to reproduce the exact values of the pairing gaps in infinite nuclear
matter. The saturation density $\rho_{0}$ of the SNM is $\rho_{0}=0.16\ {\rm
fm}^{-3}$ and asymmetric parameter $\beta$ is defined as $\beta=(N-Z)/A$ or
$\beta=(\rho_{\rm n}-\rho_{\rm p})/(\rho_{\rm n}+\rho_{\rm p})$, in which $N$
($Z$) is neutron (proton) number, $A=N+Z$ is mass number and $\rho_{\rm n}$
($\rho_{\rm p}$) refers to neutron (proton) density.

The isoscalar pairing interaction reads,

$g(\rho,\beta=0)=1-\eta_{s}(\rho/\rho_{0})^{\alpha_{s}},\\\ $ (4)

where $\eta_{s}=1$ and $\alpha_{s}=1$.

In the BCS approximation schul , the pairing matrix element reads,

$\bar{v}^{\rm pair}_{q,k\bar{k}m\bar{m}}=-V_{q}\int{\rm d}^{3}r\
g_{q}(\rho,\beta)\Psi^{\dagger}_{k}(\bm{r})\Psi^{\dagger}_{\bar{k}}(\bm{r})\Psi_{m}(\bm{r})\Psi_{\bar{m}}(\bm{r}),\\\
$ (5)

where $q$ stands for n or p, $V_{q}=V_{0}$ is the pairing strength determined
by the scattering length, which reproduces the phase shift in the low energy
region for a given cutoff energy, $g(\rho,\beta)$ from Eq. (3) or Eq. (4) is
the form factor of the isovector or isoscalar pairing interaction, and
$\Psi_{k}(\bm{r})$ is the wave function of the $k$th HF single-particle (s.p.)
level Krieger1990_NPA517-275 .

###  II.3 Pairing energy, binding energy and pairing gap in the SHF + BCS
method

In the SHF+BCS method, the binding energy $B$ of a nucleus can be written as a
sum of five parts Ryssens2015_CPC187-175 ,

$B=-E_{\rm kin}-E_{\rm Sk}-E_{\rm Coul}-E_{\rm pair}-E_{\rm corr},\\\ $ (6)

where $E_{\rm kin}$ is the kinetic energy, $E_{\rm Sk}$ is the Skyrme energy,
$E_{\rm Coul}$ is the Coulomb energy, $E_{\rm pair}$ is the pairing energy,
and $E_{\rm corr}$ is the center-of-mass correction energy. $E_{\rm pair}$ is
the main contribution of the pairing interaction to the binding energy. Other
parts of the binding energy, such as $E_{\rm Sk}$, are functions of the
density $\rho$, and the density $\rho$ is a function of the occupation
probability $v_{k}^{2}$ of single-particle levels. Therefore, the binding
energy will also be influenced by the pairing interaction through $v_{k}^{2}$.

In particular, the pairing energy $E_{\rm pair}$ can be written in the
canonical basis as

$E_{\rm pair}=\sum_{k,m>0}f_{k}u_{k}v_{k}f_{m}u_{m}v_{m}\bar{v}^{\rm
pair}_{k\bar{k}m\bar{m}},\\\ $ (7)

in which $u_{k}^{2}+v_{k}^{2}=1$, $\bar{v}^{\rm pair}_{k\bar{k}m\bar{m}}$ is
the pairing matrix element, and

$f_{k}=[1+e^{(\epsilon_{k}-\lambda_{q}-E_{C})/\mu_{q}}]^{-1/2}[1+e^{(\lambda_{q}-\epsilon_{k}-E_{C})/\mu_{q}}]^{-1/2},\\\
$ (8)

is the cutoff factor, where $\epsilon_{k}$ is the energy of the $k$th s.p.
level, $\lambda$ is the Fermi energy, $E_{C}$ is the truncation energy of the
pairing interaction, and $\mu_{q}$ is fixed to be 0.5 MeV
Krieger1990_NPA517-275 .

The set of equations that determine the occupation probability of single-
particle states $v_{k}^{2}$ are derived from the variation of

$\frac{\delta}{\delta
v_{j}}\Big{(}2\sum_{k>0}\epsilon_{k}v_{k}^{2}-\lambda_{q}\langle\hat{N}_{q}\rangle+E_{\rm
pair}\Big{)}=0,\\\ $ (9)

where $\lambda_{q}$ is the Lagrange multiplier, which is introduced to obtain
the requested mean number of protons and neutrons.

The probability of the s.p. state $|k\rangle$ and its time-reversal s.p. state
$|\bar{k}\rangle$ being occupied by one pair of neutrons or protons can be
expressed as

$v_{k}^{2}=\big{[}1-(\epsilon_{k}-\lambda_{q})/\sqrt{(\epsilon_{k}-\lambda_{q})^{2}+f_{k}^{2}\Delta_{k\bar{k}}^{2}}\
\big{]}/2,\\\ $ (10)

in which

$\Delta_{k\bar{k}}=\sum_{m>0}f_{m}u_{m}v_{m}\bar{v}^{\rm
pair}_{k\bar{k}m\bar{m}}\\\ $ (11)

is the pairing gap of the $k$th single-particle level.

For odd-$A$ nuclei, we use the blocking method of Ref. Bender2000_EPJA8-59  to
consider the odd particle. When a pair of s.p states $|k\rangle$ or
$|\bar{k}\rangle$ is chosen to be occupied or blocked, the occupied or
unoccupied probability in the BCS state is fixed to be
$u_{k}^{2}=v_{k}^{2}=1/2$, and the pairing gap of that s.p. levels is fixed to
be $\Delta_{k}=0$ MeV.

###  II.4 Numerical details

The ev8 code solves the HF+BCS equations to obtain binding energies
iteratively with the imaginary time step method Davies1980_NPA342-111 . For
even-even and odd-$A$ nuclei, we use the empirical three-point formula to
extract the OEM staggering of isotopes,

$\Delta_{n}^{(3)}(N,Z)=-\frac{\pi_{N}}{2}[B(N+1,Z)-2B(N,Z)+B(N-1,Z)],\\\ $
(12)

where $B(N,Z)$ is the binding energy of a nucleus, and $\pi_{N}$ is the parity
of the isotope with neutron number $N$.

We adopt three isospin-density dependent pairing interactions: P1
Zhang2011_SCPMA54-236 , P2 Margueron2007_PRC76-064316 ;
Margueron2008_PRC77-054309 , and P3 Bertulani2009_PRC80-027303 . For
comparison, the isoscalar interaction P4 Bertulani2009_PRC80-027303  is also
considered. We list the parameters of the pairing interactions in Table 1. One
can see that P1 and P2 have the same cut-off energy $E_{C}=40$ MeV, while P3
and P4 have a smaller cut-off energy $E_{C}=5$ MeV. As we know, effective
pairing interactions are sensitive to the values of the energy (or momentum)
cut-off (see, e.g., Fig.4 of Ref. Zhang2010_PRC81-044313  for PNM). Therefore,
the parameters are quite different for P2 and P3, especially the potential
strength $V_{0}$. This will be further analyzed in the following section.

Before large-scale calculations, we check the convergence of OEM staggering
with respect to the number of neutron wave functions. This is necessary
because the calculated OEM staggering of finite nuclei might not be right if
the basis space of wave functions is not large enough. In the ev8 code, the
number of neutron wave functions,“nwaven”, is an input parameter. We choose
neutron-rich Ca, Ni, Zr and Sn isotopes with magic or sub-magic proton number,
to calculate OEM staggering using the three-point formulas with the P1 pairing
interaction. We can see from Table 2 that neutron pairing gaps $\Delta_{n}$
converge to a certain value with the number of neutron wave functions
increasing. The absolute value of the differences for $\Delta_{n}$ from $N=90$
to $N=100$ represented by $\delta$ is smaller than 0.03 MeV for all the cases,
which guarantees the convergence. In our later study, $N=100$ is used.

In most calculations, an absolute accuracy of $10^{-7}$ MeV can be achieved
for binding energies. Since the binding energies are at the order of 1000 MeV,
the relative accuracy is better than $10^{-10}$. However, for some odd-$A$
nuclei, it is difficult to achieve such an accuracy. Therefore, the criterion
of convergence is fixed to be $10^{-5}$ MeV, which corresponds to an relative
accuracy $10^{-8}$.

Table 1: Parameters of the pairing interactions studied in the present work, in which P1, P2, and P3 are of isovector type, and P4 is of isoscalar type. $E_{C}$ is the truncated energy, $V_{0}$ is the pairing strength, $\eta_{s}$, $\alpha_{s}$, $\eta_{n}$, and $\alpha_{n}$ are the parameters of Eq. (3). Parameters | P1 | P2 | P3 | P4
---|---|---|---|---
$E_{C}$ (MeV) | 40 | 40 | 5 | 5
$V_{0}$ (MeV fm3) | 542 | 542 | 824 | 1400
$\eta_{s}$ | 0.729 | 0.664 | 0.677 | 1
$\alpha_{s}$ | 0.522 | 0.522 | 0.365 | 1
$\eta_{n}$ | 1.010 | 1.010 | 0.931 | -
$\alpha_{n}$ | 0.525 | 0.525 | 0.378 | -
Table 2: Neutron pairing gaps $\Delta_{\rm n}$ for Ca, Ni, Zr and Sn isotopes as a function of the number of neutron wave functions for the pairing interaction P1. The absolute values of the difference for $\Delta_{n}$ from $N=90$ to $N=100$ are represented by $\delta$. All the quantities are in units of MeV. | $N=80$ | $N=90$ | $N=100$ | $\delta$
---|---|---|---|---
$\Delta_{\rm n}(^{52}{\rm Ca})$ | 1.03 | 0.95 | 0.95 | 0.00
$\Delta_{\rm n}(^{74}{\rm Ni})$ | 1.23 | 1.24 | 1.22 | 0.02
$\Delta_{\rm n}(^{110}{\rm Zr})$ | 1.40 | 1.37 | 1.40 | 0.03
$\Delta_{\rm n}(^{116}{\rm Sn})$ | 1.36 | 1.38 | 1.38 | 0.00
Figure 1: Odd-even mass staggering as a function of mass number A for Ca, Ni,
Zr and Sn isotopes with the SkP force. P1 (solid red lines with triangles), P2
(solid dark blue lines with squares), and P3 (solid green lines with diamonds)
are isovector pairing interactions and P4 (solid light blue lines with stars)
is isoscalar one. The black solid circles with error bars correspond to
experimental data.

##  III Results and discussions

We calculate the binding energies of even-even and odd-$A$ Ca, Ni, Zr and Sn
isotopes with the SkP force, then derive the OEM staggering via the three-
point formula of Eq. (12). It is known that pairing gaps obtained from the OEM
staggering for those isotopes with magic or semi-magic number $N=20$, 28, 40,
50 and 82 are not reliable, such as 40Ca, 48Ca, 56Ni, 78Ni, 90Zr, 100Sn and
132Sn. Therefor, we ignore the neutron pairing gaps $\Delta_{n}$ for these
nuclei in the later analysis.

In Fig. 1, we show the OEM staggering as a function of mass number A for Ca,
Ni, Zr and Sn isotopes. The solid red (dark blue) lines with triangles
(squares) refer to the neutron pairing gaps obtained with the isovector
pairing interactions P1 (P2), which is extracted from the BHF + BCS
calculations for SNM and PNM with the mean field (free) spectrum. The solid
green lines with diamonds correspond to the predictions of the isovector
pairing interaction P3, which is obtained by fitting to the experimental OEM
staggering. While the solid light blue lines with stars denote the results
obtained with the isoscalar pairing interaction P4. The black solid circles
with error bars labelled by Exp. represent the experimental data. Generally
speaking, the neutron pairing gaps $\Delta_{\rm n}^{(3)}$ obtained with the
isovector pairing interactions P1, P2 and P3, are much closer to the
experimental data than those from the isoscalar pairing interaction P4.
Therefore, the isospin-dependence of the pairing interaction seems to be
crucial to reproduce the experimental OEM staggering.

To quantitatively evaluate the deviations of the theoretical predictions from
the experimental OEM staggering, we calculate the Root Mean Square Error
(RMSE) of the neutron pairing gaps for all the isovector pairing interactions,
except for magic or semi-magic nuclei as mentioned above. The RMSE is defined
by

$\sigma=\sqrt{\frac{\sum\limits_{i=1}^{N}(\Delta_{i,{\rm cal}}-\Delta_{i,{\rm
exp}})^{2}}{N}},\quad i=1,2,\cdots,N\\\ $ (13)

where $\Delta_{i,{\rm cal}}$ and $\Delta_{i,{\rm exp}}$ are the calculated and
the experimental OEM staggering of an isotope labeled by $i$, and $N$ is the
number of isotopes considered. We list the RMSEs in Table 3 for a better
understanding of the predictions of different pairing interactions in
comparison with the experimental data. As mentioned above, the neutron pairing
gaps for magic nuclei 40Ca, 48Ca, 56Ni, 78Ni, 90Zr, 100Sn and 132Sn are
omitted in calculating the RMSEs.

For Ni and Zr isotopes, P1 turns out to be the best among the isovector
pairing interactions. The RMSEs of the OEM staggering for Ni, Zr and Sn
isotopes are reduced by $14.3\%$, $41.0\%$, and $30.4\%$ respectively,
compared with those of P2, which is also extracted from infinite nuclear
matter with the BHF method but for free spectrum, instead of the mean field
spectrum. Moreover, the predictions of P1 for Ca isotopes are almost the same
as those of P2, and are better than those of P3 by $\sim 5\%$. One should note
that P3 is obtained by fitting to the experimental OEM staggering. In
addition, for Zr isotopes P1 can reduce the RMSEs of OEM staggering by $\sim
34.3\%$ compared with P3. For Sn isotopes, P1 can also give a better
description of OEM staggering than P2 by $30\%$, but not as good as P3. It
seems to be consistent with the conclusion of Ref. Bertulani2009_PRC80-027303
that P3 is particularly good for Sn isotopes together with the SkP force. As
for P4, it is much worse than the other three isovector pairings as we noted
earlier. Again, we stress that that neutron pairing gaps from infinite nuclear
matter can be a good constraint for constructing neutron pairing interactions
for finite nuclei.

To further investigate the difference among the three isovector pairing
interactions, we simply choose 116Sn as an example to plot in Fig. 2 the
neutron pairing interaction $V_{\rm n}=V_{0}\ g_{\rm n}(\rho,\beta)$ as a
function of density $\rho$ with fixed asymmetric parameter
$\beta=(66-50)/116\approx 0.138$. For increasing $\rho$, $V_{\rm n}$ decreases
monotonically. The curve of P3 is much different from those of P1 and P2, due
to the different energy cutoff. On the other hand, P1 and P2 look similar.
Suppose that finite nuclei have the saturation density $\rho=0.16$ fm -3, then
$V_{\rm n}$ of P2 is larger than that of P1. Correspondingly, the pairing gaps
obtained with P2 should also be larger than those obtained with P1 in Fig. 1.
The main difference between P1 and P2 is the value of $\eta_{s}$, see Table 1,
$\eta_{s}$=0.664 for P2 is slightly smaller than $\eta_{s}$=0.729 for P1.
Since $\eta_{s}$ is negative in $g_{\rm n}(\rho,\beta)$ [see Eq. (3)], smaller
$\eta_{s}$ of P2 results in a larger $V_{\rm n}$ compared with P1. Quite
remarkably, it seems that pairing gaps in ANM impose strong constraint on
those in finite nuclei, to the extent that the same pairing interactions can
be used for both cases. More studies are needed to further corroborate such a
conclusion.

To show the predictive power of the pairing interactions, we extend our
calculations to neutron rich nuclei 58Ca, 80Ni, 112Zr, 140Sn, which are
plotted in Fig. 1 as well. We have also calculated the OEM scattering of Ca
and Sn isotopes with the SLy4 force sly4 , and found that it does not predict
the pairing gaps as well as the SkP force. Generally speaking, it
underestimates the pairing gaps compared with measured data and doubles the
RMSEs than the SkP force for the three isovector pairing interactions.
Therefore, it justifies the use of SkP force in the SHF+BCS model to check the
validity of pairing interactions and it confirms the claim of Ref.
Dobaczewski1984_NPA422-103  that SkP can give better descriptions of pairing
gaps.

Table 3: RMSE of the OEM staggering of Ca, Ni, Zr and Sn isotopes in units of MeV and the relative errors. The OEM staggering for magic nuclei 40Ca, 48Ca, 56Ni, 78Ni, 90Zr, 100Sn and 132Sn are neglected in the calculation of RMSE. | P1 | P2 | $\dfrac{P2-P1}{P2}$ | P3 | $\dfrac{P3-P1}{P3}$ | P4 | $\dfrac{P4-P1}{P4}$
---|---|---|---|---|---|---|---
Ca | 0.35 | 0.33 | -6.1% | 0.37 | 5.4% | 1.17 | 70.1%
Ni | 0.30 | 0.35 | 14.3% | 0.30 | 0.0% | 0.94 | 68.1%
Zr | 0.23 | 0.39 | 41.0% | 0.35 | 34.3% | 1.00 | 77.0%
Sn | 0.16 | 0.23 | 30.4% | 0.12 | -33.3% | 0.89 | 82.0%
Figure 2: Neutron pairing interaction $V_{\rm n}=V_{0}\ g_{\rm n}(\rho,\beta)$
as a function of density $\rho$. $g_{\rm n}(\rho,\beta)$ is the neutron
pairing form factor calculated from Eq. (3). Parameters of the pairing
interactions are listed in Table 1.

##  IV Summary

We studied the OEM staggering of Ca, Ni, Zr and Sn isotopes with the effective
pairing interaction P1, together with other three types of pairing
interactions for comparison, using the Skyrme Hartree-Fock plus BCS method
with the SkP force. We showed that P1 is suitable for the description of OEM
staggering in these isotopes, especially, much better than the other pairing
interactions for Ni and Zr isotopes. For example, the pairing gaps for Zr
isotopes are $41\%$ ($34\%$) better than those obtained with P2 (P3),
accordingly. Our predictions for Ca isotopes are comparable with those of P2
and P3, which is reasonable since the parameters of P3 are obtained by fitting
to the experimental OEM staggering Bertulani2009_PRC80-027303 . For Sn
isotopes, the predictions of P1 are almost $\sim 30\%$ better than those of
P2, and are comparable with the results of P3. It is quite interesting to find
that although the isospin-density dependent pairing interaction P1 is
extracted from the bare interaction for nuclear matter within the framework of
the BHF method, it can give reasonable description of the OEM staggering for
Ca, Ni and Zr isotopes ($30\%$ better than the results obtained with P2). The
isoscalar pairing P4 is also considered and turns out to be not very good,
which shows the importance of the isospin effect in pairing interactions.

One should note that the pairing interaction P1 is obtained by fitting to the
pairing gaps of SNM and PNM with the bare interaction in the framework of the
microscopic BHF method. It is interesting to see that without any tunable
parameters, it can give pretty good description of the OEM staggering compared
with measured data. For some isotopes, e.g. Ca and Zr isotopes (or Ni
isotopes), the predictions are even better (or comparable) than those of P3,
extracted by fitting to the experimental pairing gaps. From this point of
view, it is reasonable to say that neutron pairing gaps of infinite nuclear
matter can be a good constraint of the neutron pairing interaction in finite
nuclei. In the future, our ansatz could also be tested for proton pairing
gaps.

##  V Acknowledgements

Discussions with Prof. Hiroshima Watanabe and Prof. Hagino are gratefully
acknowledged. Authors show their great thanks to Prof. Shan-Gui Zhou for his
careful reading of this article. This work was supported partially by the
National Natural Science Foundation of China under Grant No. 11775014, No.
11975096, No. 11735003, No. 11975041, and No. 11961141004.

## References

  * (1) S. S. Zhang, U. Lombardo and E. G. Zhao, Sci. Chin. Phys. Mech. Astro. 54, 236 (2011).
  * (2) S. S. Zhang, L. G. Cao, U. Lombardo, et al. Phys. Rev. C 81, 044313 (2010).
  * (3) D. J. Dean and M. Hjorth-Jensen, Rev. Mod. Phys., 75, 607 C656 (2003).
  * (4) J. M. Dong, U. Lombardo, H. F. Zhang, and W. Zuo, The Astrophysical Journal 817, 6 (2016).
  * (5) S. S. Zhang, J. P. Peng, M. S. Smith, et al. Phys. Rev. C 91, 045802 (2015).
  * (6) C. A. Bertulani, H. F. Lu and H. Sagawa, Phys. Rev. C 80, 027303 (2009).
  * (7) S. Goriely, N. Chamel and J. M. Pearson, Phys. Rev. Lett., 102, 152503 (2009).
  * (8) N. Hinohara and W. Nazarewicz, Phys. Rev. Lett., 116, 152502 (2016).
  * (9) Y. Tian, Z. Y. Ma and P. Ring, Phys. Lett. B 676, 44 (2009).
  * (10) Y. T. Rong, P. W. Zhao and S. G. Zhou, Phys. Lett. B 807, 135533 (2020).
  * (11) J. Margueron, H. Sagawa and K. Hagino, Phys. Rev. C 77, 054309 (2008).
  * (12) N. Chamel, S. Goriely and J. Pearson, Nucl. Phys. A 812, 72 (2008).
  * (13) S. S. Zhang, L. G. Cao, U. Lombardo, et al. Phys. Rev. C 93, 044329 (2016).
  * (14) J. Margueron, H. Sagawa and K. Hagino, Phys. Rev. C 76, 064316 (2007).
  * (15) L. G. Cao, U. Lombardo and P. Schuck, Phys. Rev. C 74, 064301 (2006).
  * (16) P. Bonche, H. Flocard and P. Heenen, Comput. Phys. Commun. 171, 49 (2005).
  * (17) W. Ryssens, V. Hellemans, M. Bender, et al. Comput. Phys. Commun. 187, 175 (2015).
  * (18) M. Wang, G. Audi, F. G. Kondev, et al. The AME2016 atomic mass evaluation (II). Tables, graphs and references. Chin. Phys. C 41, 030003 (2017).
  * (19) E. Chabanat, P.Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A 635 (1998) 231, Nucl. Phys. A 643 (1998) 441.
  * (20) J. Dobaczewski, H. Flocard and J. Treiner, Nucl. Phys. A 422, 103 (1984).
  * (21) U. Lombardo and H.-J. Schulze, “Superfluidity in Neutron Star Matter” in Physics of Neutron Star Interiors, Lecture Notes in Physics vol.578, pp .30-54, Eds. D. Blaschke, N. K. Glendenning and A. Sedrakian (Springer Verlag, 2001).
  * (22) S. Krieger, P. Bonche, H. Flocard, et al. Nucl. Phys. A 517, 275 (1990).
  * (23) M. Bender, K.Rutz, P.-G. Reinhard, and J. A. Maruhn, Eur. Phys. J. A 8, 59-75 (2000).
  * (24) K. T. R. Davies, H. Flocard, S. Kireger, and M. S. Weiss, Nucl. Phys. A 342, 111 (1980).


equ[][]

# The Dilemma of Choice:
Addressing Constraint Selection for Autonomous Robotic Agents

Hardik Parwana1,2, Ruiyang Wang1,2, and Dimitra Panagou2
1Both the authors have equal contribution.2Hardik Parwana, Ruiyang Wang, and
Dimitra Panagou are with the Department of Robotics, University of Michigan,
Ann Arbor, MI 48109, USA<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>work was partially sponsored by the Office of Naval
Research (ONR), under grant number N00014-20-1-2395. The views and conclusions
contained herein are those of the authors only and should not be interpreted
as representing those of ONR, the U.S. Navy or the U.S. Government.

###### Abstract

The tasks that an autonomous agent is expected to perform are often optional
or are incompatible with each other owing to the agent’s limited actuation
capabilities, specifically the dynamics and control input bounds. We encode
tasks as time-dependent state constraints and leverage the advances in multi-
objective optimization to formulate the problem of choosing tasks as selection
of a feasible subset of constraints that can be satisfied for all time and
maximizes a performance metric. We show that this problem, although amenable
to reachability or mixed integer model predictive control-based analysis in
the offline phase, is NP-Hard in general and therefore requires heuristics to
be solved efficiently. When incompatibility in constraints is observed under a
given policy that imposes task constraints at each time step in an
optimization problem, we assign a Lagrange score to each of these constraints
based on the variation in the corresponding Lagrange multipliers over the
compatible time horizon. These scores are then used to decide the order in
which constraints are dropped in a greedy strategy. We further employ a
genetic algorithm to improve upon the greedy strategy. We evaluate our method
on a robot waypoint following task when the low-level controllers that impose
state constraints are described by Control Barrier Function-based Quadratic
Programs and provide a comparison with waypoint selection based on knowledge
of backward reachable sets.

##  I Introduction

Autonomous systems are often deployed with the expectation to complete
multiple, possibly time-critical, tasks. For example, these tasks could be to
follow time-stamped waypoints returned by a high-level path planner or to
simultaneously perform environment exploration and object manipulation. The
successful completion of these tasks is dependent on several factors: (1) the
actuation capabilities of the autonomous agent, namely the dynamics and
control input bounds, (2) compatibility of the tasks, that is, whether it is
physically possible to perform both the tasks at the same time, (3) safety,
that is, whether a task can be completed while respecting state constraints
such as collision avoidance. The type of reasoning we aim to develop for
robots in this paper is: ”How should I choose which tasks to perform?”. We
argue that this objective is hard to solve in real-time and therefore
necessitates the use of heuristics. In this work, we pose tasks as, possibly
time-dependent, state constraints and define a performance metric that is
dependent on these tasks. The objective is to navigate the robotic agent in a
constrained environment, for example in presence of obstacles, to maximize the
performance metric while satisfying the strict safety constraints. A good
performance score will require the satisfaction of optional task-based
constraints and thus promotes the completion of tasks.

Several path planning algorithms such as Model Predictive Control(MPC)[1] and
RRT[2] are frequently employed to plan paths for robots in constrained
environments. However, these methods often ignore the nonlinear dynamics and
do not have the modality to specify time-dependent tasks or to choose among
them. Task or resource allocation[3, 4], mostly studied in the context of
multi-agent systems, assigns each agent a task according to its capabilities.
These methods also often ignore the dependence of dynamics and assume a user-
given cost, representative of time to completion, to assign tasks[5, 4]. An
Optimal allocation scheme considering actuation capabilities is usually
intractable due to the resulting mixed-integer nonlinear optimization problems
and therefore hierarchical approaches that decouple task allocation from path
planning have been proposed[6, 7].

Control barrier function(CBF)[8] is a mathematical tool to ensure the
satisfaction of state constraints by enforcing a bound on the rate of change
of constraint specifying barrier functions. Signal Temporal Logic (STL)[9]
specifying the user tasks have been transformed into CBFs for implementation
by a controller for successful completion of tasks[10, 11]. However, the
aforementioned papers do not provide a reasoning mechanism for the agent to
drop constraints for the sake of higher priority constraints or when tasks are
incompatible with each other. STL specifications have also been enforced
through Mixed Integer MPC formulations[12] albeit for linear dynamics only to
keep it computationally tractable.

In this work, we take a multiobjective optimization-inspired[13, 14] viewpoint
on task selection. When an optimization problem is infeasible, that is does
not accept a solution, several methods exist to find the maximal subset of
constraints that can be satisfied[15]. In general, finding this set has
exponential complexity in the number of constraints and therefore heuristics
are often proposed[16, 17]. Lagrange multipliers[18] also offer insight into
the satisfiability of the constraints and are often used to select candidate
constraints that could be made feasible, for example, in the popular active
set method[19] and interior point methods[20] for solving quadratic programs.
A related problem is to find a subset of the given cardinality of a set that
optimizes a performance metric. This problem is also known to be NP-Hard[21,
22].

We pose the task constraint selection problem in the framework of multi-
objective optimization. In general, our framework can also be solved by
reachability analysis or mixed integer MPC formulations at the expense of high
computational load. We also show that selecting constraints that help maximize
a performance metric for dynamical systems is NP-Hard. When incompatibility in
constraints is observed under a given policy that imposes task constraints at
each time step in an optimization problem, we assign a Lagrange score to each
of these constraints based on the variation in the corresponding Lagrange
multipliers over the compatible time horizon. These scores are then used to
decide the order in which constraints are dropped in a greedy strategy. We
further employ a genetic algorithm to improve upon the estimates from the
greedy strategy.

Some of the recent works address similar problems. [23] allows a robot to
execute tasks described by Linear Temporal Logic specifications. However, it
assumes holonomic robots exhibiting straight-line motions and uses a maximum
velocity-based time-to-go estimate as a measure of the reachability of the
target region. [24] and [25] perform reactive task allocation through relaxed
CBFs constraints in a CBF-QP controller by weighing the slack variables
proportional to a user-given priority of the constraint. All of the
aforementioned approaches share the characteristics of a myopic controller,
that is, they may enter regions from where the robot may not be able to
recover and finish the high-priority tasks. They also do not provide an
approach to choose between optional task constraints when they are
incompatible with each other or with the other hard constraints imposed by the
user.

The rest of the paper is structured as follows. Section IV provides our
problem formulation and shows that it is NP-Hard. In Section V, we show how
our problem can be solved using backward reachable sets or mixed-integer
nonlinear MPC. We then propose our heuristic method based on the Lagrange
multiplier in Section V-C and present simulation results in Section VI.

##  II Preliminaries

###  II-A Notation

The set of real numbers is denoted as $\mathbb{R}$ and the non-negative real
numbers as $\mathbb{R}^{+}$. Given $x\in\mathbb{R}^{n}$, $||x||$ denotes
$L_{2}$ norm of $x$. Let $\mathbb{I}(y),y\in\mathbb{R}$ be step identity
function such that $\mathbb{I}(y)=1$ if $y\geq 0$ and $\mathbb{I}(y)=0$
otherwise. The interior and boundary of a set $\mathcal{C}$ are denoted by
$\textrm{Int}(\mathcal{C})$ and $\partial\mathcal{C}$. For
$a\in\mathbb{R}^{+}$, a continuous function
$\alpha:[0,a)\rightarrow[0,\infty)$ is a class-$\mathcal{K}$ function if it is
strictly increasing and $\alpha(0)=0$. Furthermore, if $a=\infty$ and
$\lim_{r\rightarrow\infty}\alpha(r)=\infty$, then it is called
class-$\mathcal{K}_{\infty}$.

##  III Control Barrier Functions

###  III-A Lagrange multiplier

Consider the optimization problem

$\displaystyle\min_{x\in\mathbb{R}^{n}}\quad$ $\displaystyle J(x)$ (1) s.t.
$\displaystyle g(x)\geq 0$ (2)

where $J:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and
$g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{q}$ are the cost and constraint
functions respectively. The Hamiltonian[18] is given by

$\displaystyle H=J+\lambda^{T}g$ (3)

where $\lambda\in\mathbb{R}^{q}$ is the Lagrange multiplier. A necessary
condition for a solution $x^{*}$ to be optimal is

$\displaystyle\lambda^{*T}g(x^{*})=0\textup{ and }\lambda^{*}>0\textup{ if
}g(x^{*})>0$ (4)

$\lambda$ is also the sensitivity of the cost function with respect to
constraint bounds. A larger value of $\lambda^{*}$ for a constraint would
therefore lead to a larger change in cost function if that constraint were
relaxed.

###  III-B Genetic Algorithm (GA)

We give a brief overview of the Genetic Algorithm(GA)[26] that we use in
Section V-C. Consider an unconstrained optimization problem

$\displaystyle\min_{y\in\mathbb{R}^{n}}\quad$ $\displaystyle F(y)$ (5)

where $F:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is the objective function and
y$\in\mathbb{R}^{n}$ is a chromosome defined by n state variables. In order to
minimize $F(x)$, GA keeps a population of $n^{\prime}$ potential chromosomes:

$\displaystyle Y_{GA}=y_{1},y_{2},...,y_{n^{\prime}}$ (6)

Each chromosome has an associated fitness score given by the function
$F(y_{i}):\mathbb{R}^{n}\rightarrow\mathbb{R}$. In every iteration of GA, new
chromosomes are first sampled from the population with probability
proportional to their fitness scores. Next, in the crossover process, sampled
chromosomes exchange part of their states randomly with other chromosomes
inside the sampled population giving rise to $n^{\prime}$ new chromosomes
$Y_{GA}^{{}^{\prime}}$. Finally, in the last process of mutation, The above
process is repeated until the fitness score converges or the maximum number of
iterations(GA steps) is met.

###  III-C System Description

Consider a discrete-time nonlinear dynamical system

$\displaystyle{x}_{t+1}=f(x_{t})+g(x_{t})u,$ (7)

where $x_{t}\in\mathcal{X}\subset\mathbb{R}^{n}$ and
$u_{t}\in\mathcal{U}\subset\mathbb{R}^{m}$ represent the state and control
input, and $f:\mathcal{X}\rightarrow\mathbb{R}^{n}$ and
$g:\mathcal{X}\rightarrow\mathbb{R}^{m}$ are locally Lipschitz continuous
functions. The set $\mathcal{S}(t)$ of allowable states at time $t$ is
specified as an intersection of $N$ sets
$\mathcal{S}_{i}(t),i\in\\{1,2,..,N\\}$, each of which is defined as the
0-superlevel set of a continuously differentiable function
$h_{i}:\mathbb{R}^{+}\times\mathcal{X}\rightarrow\mathbb{R}$ as:

$\displaystyle\mathcal{S}_{i}(t)$
$\displaystyle\triangleq\\{x\in\mathcal{X}:h_{i}(t,x)\geq 0\\},$ (8a)
$\displaystyle\partial\mathcal{S}_{i}(t)$
$\displaystyle\triangleq\\{x\in\mathcal{X}:h_{i}(t,x)=0\\},$ (8b)
$\displaystyle\text{Int}(\mathcal{S}_{i})(t)$
$\displaystyle\triangleq\\{x\in\mathcal{X}:h_{i}(t,x)>0\\}$ (8c)

Let $h_{i},i\in\\{1,..,n_{H}\\}$ be the hard constraints corresponding to, for
example, safety specifications or high-priority tasks that must be completed.
Let $h_{j},j\in\\{1,..,n_{S}\\}$ be the soft constraints that an agent is
allowed to drop and $N=n_{H}+n_{S}$. Let

$\mathcal{S}_{H}=\cap_{i}\mathcal{S}_{i},\quad\mathcal{S}_{S}=\cap_{j}\mathcal{S}_{j}$
(9)

be the intersection of safe sets corresponding to hard and soft constraint
barrier functions, $h_{i}$ and $h_{j}$, respectively.

Finally, the low-level controller that we consider in this work to illustrate
our algorithms is based on the notion of control barrier function[8]. The
following condition on $h_{i}$, called the CBF derivative condition, is
sufficient for the invariance of the set $\mathcal{S}_{i}$ if it can be
satisfied for all $x_{t}\in\mathcal{X}$

$\displaystyle h_{i}(t+1,x_{t+1})-h_{i}(t,x_{t})\geq\alpha_{i}h_{i}(t,x_{t})$
(10)

where $\alpha\in\mathbb{R}^{+}$ is a fixed constant that is usually treated as
a tuning parameter of the controller. The following quadratic program(QP)
based controller[27, 28] is used to enforce satisfaction of all the hard and
soft constraints $h_{i},h_{j}$.

$\displaystyle\min_{u_{t}\in\mathcal{U}}\quad$ $\displaystyle
J(u)=(u_{t}-u_{d_{t}})^{T}P(u_{t}-u_{d_{t}})$ s.t. $\displaystyle
h_{i}(t+1,x_{t+1})\geq(1-\alpha_{i})h_{i}(t,x_{t})~{}\forall
i\in\\{1,..,n_{H}\\}$ $\displaystyle
h_{j}(t+1,x_{t+1})\geq(1-\alpha_{j})h_{j}(t,x_{t})~{}\forall
j\in\\{1,..,n_{S}\\}$ (11)

where $u_{d_{t}}\in\mathbb{R}^{m}$ is a desired or nominal control input
provided by the user.

##  IV Problem Formulation

Let
$r(t,x_{t},u_{t}):\mathbb{R}^{+}\times\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$
be the reward function that determines the cost of choosing $u_{t}$ at state
$x_{t}$ and time $t$. In general, it may be designed to promote the completion
of optional tasks by depending on functions $h_{j},j\in\\{1,..,n_{S}\\}$.

######  Objective 1.

Consider the dynamical system (7). Suppose the initial state
$x_{1}\in\mathcal{S}_{H}\cap\mathcal{S}_{S}$. For a given time horizon $T$,
design control input $u_{t}:\mathcal{X}\rightarrow\mathcal{U},t\in[1,2,..,T]$
such that

  1. 1.

$x_{t+1}\in\mathcal{S}_{H}(t),t\in\\{1,...,T\\}$

  2. 2.

$J=\sum_{t=0}^{t=T}r(t,x_{t},u_{t})$ is maximized

We now provide an alternate formulation of Objective 1 that is relevant when
robots are deployed with a structured controller such as (11). These control
architectures are usually modular enough to consider an arbitrary number of
constraints in their structure. Suppose each combination of constraints is
represented as a different policy
$\pi_{k}:\mathcal{X}\rightarrow\mathcal{U},k\in\\{1,...,N_{\pi}\\}$. Let the
set of all policies that can be used be denoted as

$\displaystyle\Pi=\\{\pi_{1},\pi_{2},...,\pi_{N_{\pi}}\\}$ (12)

and $\Pi(x)$ denotes a multi-valued mapping. In this work, for illustration,
we consider $\pi_{k}$ to be given by (11) but with different combination of
constraints $h_{j},j\in\\{1,..,n_{S}\\}$.

######  Objective 2.

Consider the dynamical system (7) and the policy set $\Pi$. Suppose the
initial state $x_{1}\in\mathcal{S}_{H}\cap\mathcal{S}_{S}$. Design control
input $u_{t}:\mathcal{X}\rightarrow\mathcal{U},t\in[1,2,..,T]$ to solve the
following problem

$\displaystyle u_{1},..,u_{T}=\quad$
$\displaystyle\operatorname*{arg\,max}_{u_{1},..,u_{T}}J=\sum_{t=1}^{T}r(t,x_{t},u_{t})$
(13) s.t. $\displaystyle x_{t+1}=f(x_{t})+g(x_{t})u_{t}$ (14) $\displaystyle
u_{\tau}\in\Pi(x_{\tau})$ (15) $\displaystyle x_{t+1}\in\mathcal{S}_{H}(t)$
(16)

To ensure the problem is well-posed when the QP (11) does not admit a
solution, $r$ may be designed to approach infinity whenever (11) is
infeasible. Note that policy switching is required because each individual
policy $\pi_{k}$ may not be optimal. If $\pi_{k}$ were MPC controllers,
corresponding to a different combination of optional constraints, with
infinite time horizons, then a single policy can be chosen. However, policies
such as CBF-QP(11) or PID are prone to sub-optimality due to their myopic
nature[28, 29] and sensitivity to the controller gains, such as
$\alpha_{i},\alpha_{j}$ in (11) or PID gains in a PID controller. In
principle, we may have multiple policies in set $\Pi$ that try to enforce the
same constraint but with different $\alpha_{i},\alpha_{j}$. We now show that
the above-mentioned objectives are hard to solve in practice.

######  Theorem 1.

Objective 1 is NP-Hard.

###### Proof.

We provide two proofs.

Approach 1: We show that Problem 1 is NP-Hard in general by considering a
special case of the problem. Consider the scenario when $T=1$ and the problem
is infeasible. Further, let the cost function
$r(t,x,u)=-\sum_{i}\frac{1}{\mathbb{I}(h_{i}(t,x))}+\sum_{j}\mathbb{I}(h_{j}(t,x))$.
In this case, $r$ maximizes the number of soft constraints that can be
satisfied simultaneously while satisfying all hard constraints. Therefore,
this special case is an instance of the maximal feasible subset selection
problem that is known to be NP-Hard[15, 17].

Approach 2: Consider again $T=1$ and the index set $V=\\{1,2,..,n_{S}\\}$. Let
$E\subset V$. Consider the cost function $r(t,x_{t},u_{t})=\sum_{i\in
E}\frac{1}{\mathbb{I}(g_{i}(t,x))}+m(t,x_{t},u_{t})$ where $m$ is any
function. and the optimization problem

$\displaystyle\min_{E\subset V}f(E)$ (17)

This cost function approaches infinity if a constraint included in the index
set $E$ is not satisfied by x. Therefore, the above optimization chooses a
feasible subset if one exists and if it optimizes the cost $m(t,x_{t},u_{t})$.
(17) is an instance of a subset selection problem that is known to be NP-
Hard[21, 22]. ∎

Note that the Theorem 1 proves that Objective 1 is NP-Hard in general. For
some criterion functions $r$, exact solutions may be computable in sub-
polynomial time. Exploiting structures in dynamical systems that achieve this
complexity will be considered in future work.

##  V Methodology

In this section, we show different methods to solve Problem 1. Exact solutions
can be obtained by computing backward reachable sets or through model
predictive control. Owing to the computationally prohibitive nature of both
methods even for offline computation, we propose heuristics that can be used
for real-time but sub-optimal decision-making subject to an agent’s actuation
capabilities. Moreover, for some classes of systems such as decentralized
multi-agent systems, reachability and MPC and not viable options due to the
inability to control other agents in the environment.

###  V-A Reachability based Solution

Suppose the backward reachable set of each unsafe region $h_{i}(t,x)<0$ is
$\mathcal{UBR}_{i}(t,x)$. This computation is usually only amenable to offline
computation even for simple systems[30]. Let
$\mathcal{X}_{j}(t)=\\{x|h(t,x)\geq 0\\}$. The agent chooses to satisfy the
constraint $h_{j}(t,x)\geq 0$ if

$\displaystyle X_{j}(t)\cap\mathcal{UBR}_{i}(t)=\emptyset~{}~{}\forall
t>0,\forall i\in\\{1,..,n_{H}\\}$ (18)

###  V-B Model Predictive Control

Here, we show that Objective 1 can be formulated as an MPC problem using the
big-M method[31]. Let $M>>1$ and $\mathbb{I}$ be a binary variable. Consider
the constraint

$\displaystyle h_{j}(t,x)+M\mathbb{I}_{j}\geq 0$ (19)

In an optimization problem, for $\mathbb{I}_{j}=0$, $h_{j}(t,x)\geq 0$ is
imposed as a hard constraint, and for $\mathbb{I}_{j}=1$, it is, effectively,
ignored. Therefore the following optimization problem solves Objective 1.

$\displaystyle\min_{u_{1},..,u_{T}}\quad$ $\displaystyle
J=\sum_{t=1}^{T}r(t,x,u)$ (20) s.t. $\displaystyle x_{t+1}=f(x_{t},u_{t})$
(21) $\displaystyle h_{i}(t,x_{t})\geq 0$ (22) $\displaystyle
h_{j}(t,x_{t})+M\mathbb{I}_{j}\geq 0$ (23)

###  V-C Heuristics

In this section, we present our heuristic approach to solving the
combinatorial NP-Hard Objectives 1 and 2 Consider the point-wise optimal QP
controller (11) that enforces $N$ state constraints. Let the trajectory
resulting from forward propagating the closed-loop system (7) be given by
$(x_{t},u_{t}),t\in[0,..,T^{\prime}]$ where $T^{\prime}\in\\{0,..,\infty\\}$
is such that a solution to (11) exists for all $t<T^{\prime}$ and does not
exist for $t=T^{\prime}$.

Let $\lambda_{j}$ be the Lagrange multiplier of the inequality constraint
$\dot{h}_{j}+\alpha_{j}h_{j}\geq 0$ resulting from solving (11) with a given
combination of soft constraints at time $t$. We assign a Lagrange score
$l_{j}$ to each constraint $j$ as follows

$\displaystyle l_{j}=\sum_{\tau=1}^{T}\lambda_{j}(\tau)$ (24)

1:$x_{1}$ $\triangleright$ Current state of the system

2:$\mathcal{S}_{S}^{{}^{\prime}}=\\{1,..,n_{S}\\}$ $\triangleright$ Constraint
Index to consider

3:for $j\in\\{1,..,n_{S}\\}$ do

4: $l_{j}\leftarrow 0$

5:end for

6:for $t$ in $1,...,T$ do

7: $u_{t},\lambda\leftarrow\pi_{k}(t,x_{t})$

8: $x_{t+1}=f(x_{t})+g(x_{t})u_{t}$

9: if $\pi_{k}(t,x_{t})$ is Infeasible then

10: $j^{\prime}$ = $\operatorname*{arg\,max}_{j}(l_{j})$

11: $\mathcal{S}^{{}^{\prime}}_{S}=\mathcal{S}^{{}^{\prime}}_{S}\setminus
j^{\prime}$

12: Go back to line 3

13: end if

14: for $j\in\mathcal{S}^{{}^{\prime}}_{S}$ do

15: $l_{j}\leftarrow l_{j}+\lambda_{j}$

16: end for

17:end for

18:Return $\mathcal{S}^{{}^{\prime}}_{S}$

Algorithm 1 Greedy Algorithm

Algorithm 18 presents a simple greedy algorithm to select a subset of
constraints by dropping them in order of decreasing the Lagrange score. Since
greedy approaches are prone to producing sub-optimal solutions and are not
guaranteed to give a feasible sub, we also present Algorithm 17 which uses a
genetic algorithm introduction to improve the estimate of an optimal feasible
subset of constraints. Each chromosome is a binary tuple of constraints
$(b_{1},b_{2},..,b_{n_{s}}),b_{j}\in\\{0,1\\}$ such that $b_{j}=1$ implies
$j^{th}$ soft constraint is imposed and vice versa. In line 5, we first
initialize a population of chromosomes in a similar fashion as Algorithm 18 by
successively removing constraints based on their Lagrange scores. The
population of constraint subsets is now passed through an iterative genetic
evolution. In line 10, the fitness score, $F(s)$ is derived from the reward
function $r$ and the Lagrange score of the imposed soft constraints for each
subset selection in the genetic population. The genetic evolution then is
guided by the fitness score of subsets in the population to update the new
population in line 16.

Note that language variables have been used extensively in, for example,
active set[19] and interior point methods[20] to find a feasible solution to
an optimization problem. Our novelty lies in utilizing Lagrange variables to
evaluate the satisfaction of a constraint over a time horizon for an evolving
system such as (7) and using it to find a feasible subset of constraints.

1:$x_{1}$ $\triangleright$ Current state of the system

2:$\mathcal{S}_{S}^{{}^{\prime}}=\\{1,..,n_{S}\\}$ $\triangleright$ Constraint
Index to consider

3:$K_{sample}$ $\triangleright$ Population size

4:$N_{steps}$ $\triangleright$ Number of evolutionary iterations

5:$\mathcal{S}^{sample}_{S}\leftarrow\textup{sample\\_init}(K_{sample},\mathcal{S}^{{}^{\prime}}_{S})$

6:$\mathcal{S}_{S}^{best}\leftarrow\emptyset$ $\triangleright$ Initialize the
optimal subset

7:$F_{min}=\infty$ $\triangleright$ Initialize the least observed Fitness
score

8:for $step$ in $N_{steps}$ do

9: for $s$ in $\mathcal{S}^{sample}_{S}$ do

10: $F(s)\leftarrow\textup{roll\\_out}(\pi_{s},x_{1}))$

11: if $F(s)<F_{min}$ then

12: $F_{min}\leftarrow F(s)$

13: $\mathcal{S}_{S}^{best}\leftarrow s$

14: end if

15: end for

16:
$\mathcal{S}^{sample}_{S}\leftarrow\textup{genetic\\_evolve}(\mathcal{S}^{sample}_{S},F)$

17:end for

Algorithm 2 Evolutionary Soft Constraints Selection

##  VI Simulation Results

In this section, we illustrate the efficacy of our algorithm through
simulations. Consider the scenario shown in Fig. 1 in which a robot is asked
to track a series of time-stamped waypoints provided by the user or any high-
level planner. The waypoints are formatted in CBF constraints as:

$\displaystyle
h_{i}(x_{t},t)=B_{i}(t)\left(d_{max}^{2}-(x_{t}-x_{i})^{T}(x_{t}-x_{i})\right)$
(25)

where $B_{i}(t)$ is the activation function for the waypoint, $x_{t}$ is the
current state of the robot, $x_{i}$ is the waypoint’s position, and $d_{max}$
is the maximum allowable distance from the waypoint. The collision avoidance
with static obstacles, approximated as series of circles, are formulated using
the following barrier function

$\displaystyle h_{j}(x_{t})=(x_{t}-x_{j})^{T}(x_{t}-x_{j})-d_{min}^{2}$ (26)

where $x_{t}$ is the current state of the robot, $x_{o}$ is the obstacle’s
position, and $d_{min}$ is the radius of the obstacle.

Figure 1: A constrained environment with black static obstacles and known
disturbances in light blue. The robot, starting at waypoint 1, needs to track
the series of waypoints in blue in order of increasing number, and finally,
reach the red waypoint.

The QP controller (11) is used to generate control inputs for the robot with a
single integrator dynamics. Each blue waypoint is encoded as a time-dependent
soft constraint such that the robot should be in the circular regions for a
finite amount of time. The final waypoint, red in color, and collision
avoidance with black obstacles are encoded as hard constraints that must
always be satisfied.

The environment is affected by certain known disturbances, whose magnitude is
represented with light blue color in Fig. 1. The nonlinear disturbance may
lead to a collision in presence of input bounds if the robot navigates to one
of the shaded area. Thus, the robot may have to drop some of the constraints
corresponding to blue waypoints to be able to reach the last red waypoint
while remaining safe.

We present results for four cases in Fig. 2. In Case 1, the optimal solution
is obtained by computing backward reachable sets as mentioned in Section V-A
and shown in the first row of Fig.2 for three different magnitudes of
disturbance. In case 2, with reward $r$ designed to maximize the number of
regions visited, we solve Objective 2 with an exhaustive search by evaluating
all possible QP controllers of the form (11) with different combinations of
soft constraints but without switching controllers at each time step. The
results are shown in the second row of Fig.2. In case 3, we use our proposed
heuristic method from Section V-C. We find the best fixed QP controller using
Algorithm 17 for a time horizon $T$ large enough to complete navigation to all
waypoints and the results are shown in the third row of Fig. 2. Finally, we
apply Algorithm 17 in a receding horizon fashion, executing every $T/10$ time
steps for a small time horizon $T/5$ for efficient online implementation. The
results are shown in the last row of Fig. 2.

|  No. of Waypoints Arrived |  No. of Iterations
---|---|---
No. Waypoints |  Distur-bance Level |  Case 1 |  Case 2 |  Case 3 |  Case 2 |  Case 3
7 |  Low |  7 |  7 |  6.5 |  64 |  48
Medium |  7 |  7 |  6.4 |  64 |  48
High |  7 |  4 |  3.8 |  64 |  48
9 |  Low |  9 |  8 |  7.6 |  256 |  64
Medium |  8 |  7 |  6.6 |  256 |  64
High |  7 |  5 |  4.3 |  256 |  64
15 |  Low |  15 |  13 |  10.3 |  16384 |  224
Medium |  14 |  9 |  7.8 |  16384 |  224
High |  13 |  6 |  5.4 |  16384 |  224
TABLE I: Test Results from three different test environments to comparing
solutions from Case 1: reachability analysis, Case 2: exhaustive CBF search,
and Case 3: Algorithm 17.

Table I analyses the performance of different methods with respect to the
number of waypoints in the environment, which determines the number of soft
constraints. The number of iterations required in the exhaustive search for
finding the best CBF controller grows exponentially in the number of
constraints. Algorithm 17 on the other hand, with the maximum number of
allowed mutations reported in the table, still resulted in a similar
performance (number of waypoints covered) as the exhaustive search. All the
GA-based results report the average performance over 100 runs due to
randomness involved in the mutation process. We used 4 chromosomes for the
data reported in Table I.

(

a)

(

b)

(

c)

(

d)

Figure 2: The trajectory calculated from (a) Case 1: Reachability, (b) Case 2:
Offline CBF, (c) Case 3: Algorithm 2,
(d) Case 4: Algorithm 2 with the receding horizon. Each column corresponds to
a different maximum magnitude of disturbance, increasing from left to right.
The shaded blue regions show disturbances affecting the robot. The waypoints
that were reached are shown in green. The colormap represents the timestamp.

Table. II compares the performance of Algorithm 18, Algorithm 17 with zero
mutations, Algorithm 17 with 16 mutations levels, and finally with Algorithm
17 applied in receding horizon fashion with small time horizon. Based on the
number of waypoints covered, we see that mutations help improve the estimates
of a pure greedy strategy. The receding horizon, although does not show the
best performance, can be implemented for real-time implementation. Based on
all the results.

Based on the results reported in Tables I and II, we observe that Algorithm 17
drops fewer constraints and significantly outperforms Algorithm 18 when
disturbance level is high. Although the receding horizon Algorithm 17 is
myopic as it considers a shorter time horizon, it results in reasonably well-
performing trajectories.

|  No. of Waypoints Arrived
---|---
No. Sets |  Disturb Level |  Case 1 |  Case 2 |  Case 3 |  Case 4
9 |  Low |  8.0 |  5.6 |  7.6 |  6.3
Medium |  2.0 |  4.3 |  6.6 |  5.1
High |  2.0 |  1.6 |  4.3 |  2.9
TABLE II: Test Result comparing Case 1: Algorithm 1, Case 2: Algorithm 2 with
0 mutations, Case 3: Algorithm 2 with infinite time horizon, Case 4: Algorithm
2 with receding time horizon.

##  VII Conclusion

We formulated robots tasks specifications as time-dependent state constraints.
We then introduced and addressed a new decision-making problem in which the
autonomous robot is allowed to drop soft constraints permanently to ensure
that hard constraints can be satisfied. We showed that this problem is NP-Hard
and leveraged ideas from multi-objective optimization literature to propose
heuristics based on a variation of Lagrange variables of each constraint over
a time horizon. Evaluation of our algorithm is then done for a test case of
robot waypoint following in presence of obstacles and disturbance. Future work
will involve the evaluation of our method on multi-agent systems and formal
analysis of the convergence guarantees of the proposed heuristic.

## References

  * [1] C. Liu, S. Lee, S. Varnhagen, and H. E. Tseng, “Path planning for autonomous vehicles using model predictive control,” in _2017 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2017, pp. 174–179.
  * [2] S. M. LaValle and J. J. Kuffner, “Rapidly-exploring random trees: Progress and prospects: Steven m. lavalle, iowa state university, a james j. kuffner, jr., university of tokyo, tokyo, japan,” _Algorithmic and computational robotics_ , pp. 303–307, 2001.
  * [3] M.-L. Lee, S. Behdad, X. Liang, and M. Zheng, “Task allocation and planning for product disassembly with human–robot collaboration,” _Robotics and Computer-Integrated Manufacturing_ , vol. 76, p. 102306, 2022.
  * [4] A. Khamis, A. Hussein, and A. Elmogy, “Multi-robot task allocation: A review of the state-of-the-art,” _Cooperative robots and sensor networks 2015_ , pp. 31–51, 2015.
  * [5] M. Darrah, W. Niland, and B. Stolarik, “Multiple uav dynamic task allocation using mixed integer linear programming in a sead mission,” in _Infotech@ Aerospace_ , 2005, p. 7164.
  * [6] S. Choudhury, J. K. Gupta, M. J. Kochenderfer, D. Sadigh, and J. Bohg, “Dynamic multi-robot task allocation under uncertainty and temporal constraints,” _Autonomous Robots_ , vol. 46, no. 1, pp. 231–247, 2022.
  * [7] A. Tika, N. Gafur, V. Yfantis, and N. Bajcinca, “Optimal scheduling and model predictive control for trajectory planning of cooperative robot manipulators,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 9080–9086, 2020\.
  * [8] 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.
  * [9] O. Maler and D. Nickovic, “Monitoring temporal properties of continuous signals,” in _Formal Techniques, Modelling and Analysis of Timed and Fault-Tolerant Systems: Joint International Conferences on Formal Modeling and Analysis of Timed Systmes, FORMATS 2004, and Formal Techniques in Real-Time and Fault-Tolerant Systems, FTRTFT 2004, Grenoble, France, September 22-24, 2004. Proceedings_. Springer, 2004, pp. 152–166.
  * [10] L. Lindemann and D. V. Dimarogonas, “Control barrier functions for signal temporal logic tasks,” _IEEE control systems letters_ , vol. 3, no. 1, pp. 96–101, 2018.
  * [11] ——, “Control barrier functions for multi-agent systems under conflicting local signal temporal logic tasks,” _IEEE control systems letters_ , vol. 3, no. 3, pp. 757–762, 2019.
  * [12] V. Raman, A. Donzé, M. Maasoumy, R. M. Murray, A. Sangiovanni-Vincentelli, and S. A. Seshia, “Model predictive control with signal temporal logic specifications,” in _53rd IEEE Conference on Decision and Control_. IEEE, 2014, pp. 81–87.
  * [13] R. T. Marler and J. S. Arora, “Survey of multi-objective optimization methods for engineering,” _Structural and multidisciplinary optimization_ , vol. 26, pp. 369–395, 2004.
  * [14] K. Deb and K. Deb, “Multi-objective optimization,” in _Search methodologies: Introductory tutorials in optimization and decision support techniques_. Springer, 2013, pp. 403–449.
  * [15] J. W. Chinneck, “The maximum feasible subset problem (maxfs) and applications,” _INFOR: Information Systems and Operational Research_ , vol. 57, no. 4, pp. 496–516, 2019.
  * [16] P. Sadegh, “A maximum feasible subset algorithm with application to radiation therapy,” in _Proceedings of the 1999 American Control Conference (Cat. No. 99CH36251)_ , vol. 1. IEEE, 1999, pp. 405–408.
  * [17] J. W. Chinneck, _Feasibility and Infeasibility in Optimization:: Algorithms and Computational Methods_. Springer Science & Business Media, 2007, vol. 118.
  * [18] D. E. Kirk, _Optimal control theory: an introduction_. Courier Corporation, 2004.
  * [19] E. Wong, _Active-set methods for quadratic programming_. University of California, San Diego, 2011.
  * [20] F. A. Potra and S. J. Wright, “Interior-point methods,” _Journal of computational and applied mathematics_ , vol. 124, no. 1-2, pp. 281–302, 2000\.
  * [21] C. Qian, Y. Yu, and Z.-H. Zhou, “Subset selection by pareto optimization,” _Advances in neural information processing systems_ , vol. 28, 2015.
  * [22] G. Davis, S. Mallat, and M. Avellaneda, “Adaptive greedy approximations,” _Constructive approximation_ , vol. 13, pp. 57–98, 1997.
  * [23] D. Gundana and H. Kress-Gazit, “Event-based signal temporal logic tasks: Execution and feedback in complex environments,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 4, pp. 10 001–10 008, 2022.
  * [24] G. Notomista, S. Mayya, S. Hutchinson, and M. Egerstedt, “An optimal task allocation strategy for heterogeneous multi-robot systems,” in _2019 18th European Control Conference (ECC)_. IEEE, 2019, pp. 2071–2076.
  * [25] C. Wang, S. Zhu, B. Li, L. Song, and X. Guan, “Time-varying constraint-driven optimal task execution for multiple autonomous underwater vehicles,” _IEEE Robotics and Automation Letters_ , 2022.
  * [26] S. J. Russell, _Artificial intelligence a modern approach_. Pearson Education, Inc., 2010.
  * [27] 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.
  * [28] H. Parwana and D. Panagou, “Recursive feasibility guided optimal parameter adaptation of differential convex optimization policies for safety-critical systems,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 6807–6813.
  * [29] M. H. Cohen and C. Belta, “Approximate optimal control for safety-critical systems with control barrier functions,” in _2020 59th IEEE conference on decision and control (CDC)_. IEEE, 2020, pp. 2062–2067.
  * [30] S. Bansal, M. Chen, S. Herbert, and C. J. Tomlin, “Hamilton-jacobi reachability: A brief overview and recent advances,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_. IEEE, 2017, pp. 2242–2253.
  * [31] M. Cococcioni and L. Fiaschi, “The big-m method with the numerical infinite m,” _Optimization Letters_ , vol. 15, no. 7, pp. 2455–2468, 2021.