# The First TESS Self-Lensing Pulses: Revisiting KIC 12254688

Nicholas M. Sorabella Lowell Center for Space Science and Technology, Lowell,
MA 01854  Department of Physics and Applied Physics, University of
Massachusetts Lowell, Lowell, MA 01854  Silas G.T. Laycock Lowell Center for
Space Science and Technology, Lowell, MA 01854  Department of Physics and
Applied Physics, University of Massachusetts Lowell, Lowell, MA 01854
Dimitris M. Christodoulou Lowell Center for Space Science and Technology,
Lowell, MA 01854  Sayantan Bhattacharya Lowell Center for Space Science and
Technology, Lowell, MA 01854  Department of Physics and Applied Physics,
University of Massachusetts Lowell, Lowell, MA 01854

###### Abstract

We report the observations of two self-lensing pulses from KIC 12254688 in
Transiting Exoplanet Survey Satellite (TESS) light curves. This system,
containing a F2V star and white-dwarf companion, was amongst the first self-
lensing binary systems discovered by the Kepler Space Telescope over the past
decade. Each observed pulse occurs when the white dwarf transits in front of
its companion star, gravitationally lensing the star’s surface, thus making it
appear brighter to a distant observer. These two pulses are the very first
self-lensing events discovered in TESS observations. We describe the methods
by which the data were acquired and detrended, as well as the best-fit binary
parameters deduced from our self-lensing+radial velocity model. We highlight
the difficulties of finding new self-lensing systems with TESS, and we discuss
the types of self-lensing systems that TESS may be more likely to discover in
the future.

Compact binary stars (283), Gravitational microlensing (672), White Dwarf
Stars (1799)

††journal: ApJL††software: EMCEE (Goodman & Weare, 2010; Foreman-Mackey et
al., 2013), Python module corner.py (Foreman-Mackey, 2016; Hogg & Foreman-
Mackey, 2018), Python SciPy module optimize (https://scipy.org, Van Rossum &
Drake, 1995), Python package Lightkurve (Cardoso et al., 2018), Mathematica,
v. 13.3 (https://www.wolfram.com/mathematica), MATLAB, v. R2023b
(https://www.mathworks.com).

##  1 Introduction

Self-lensing binary systems are a subset of eclipsing binary systems in which
at least one component is a compact object, such as a white dwarf (WD),
neutron star (NS), or black hole (BH). When the compact object transits in
front of its (typically noncompact) companion, it gravitationally lenses the
companion’s surface, resulting in a symmetric flare in the system’s light
curve (Maeder, 1973; Witt & Mao, 1994; Marsh, 2001; Rahvar et al., 2011;
Sorabella et al., 2022a). Self-lensing pulses can be very useful in estimating
binary parameters, as the pulse profile depends on the masses of the compact
object and star, their separation during the transit, and the orbital
inclination of the system. Given that these pulses will occur every time the
compact object transits its companion, self-lensing pulses are periodic, thus,
they provide a way of also measuring the orbital period of the binary without
the need for spectroscopic measurements. The profile of the lensing pulse can
also provide information regarding the star being lensed, such as its radius
and limb-darkening coefficient. However, studies of self-lensing systems
remain scarce, as only five stellar-mass self-lensing systems (and a sixth
unconfirmed candidate) have been discovered to date (Kruse & Agol, 2014;
Kawahara et al., 2018; Masuda et al., 2019). All of these systems were
discovered with the Kepler Space Telescope and all were found to be quite
similar; each binary system contains a WD that gravitationally lenses a main-
sequence companion not too dissimilar to that of our Sun.

As a spiritual successor to the Kepler Space Telescope, the Transiting
Exoplanet Survey Satellite (TESS) (Ricker et al., 2010) is a very important
tool in the search for new self-lensing pulses. While Kepler had the advantage
of obtaining more data for each target, as it was always observing the same
patch of sky for 4 years, TESS has the distinct advantage of searching far
more area of the sky. TESS images a $24^{\circ}\times 96^{\circ}$ sector of
the sky for 27 days before moving on to another sector. This drastically
increases the number of stars observed and the potential for discovery of more
self-lensing systems (as well as extrasolar systems). It has been predicted by
Wiktorowicz et al. (2021) that TESS could find anywhere from 60-170 self-
lensing sources throughout its operational lifetime, increasing the number of
known self-lensing systems by a factor of 10-30. Combined with other surveys
such as the Zwicky Transient Facility (ZTF) and the upcoming Legacy Survey of
Space and Time (LSST), the number of detectable self-lensing systems could be
of the same order as the population of X-ray binaries in the Milky Way.

We have begun our own search for new self-lensing pulses in the existing TESS
data, although we have yet to find any new systems. Instead, we have found an
abundance of asteroids crossing in front of background stars, giving false-
positive signals that mimic self-lensing pulses Sorabella et al. (2023). We
also checked whether TESS has observed any of the known self-lensing WD
systems (Kawahara et al., 2018, hereafter KMM18) at times when new pulses were
expected to occur. This search highlighted a disadvantage of TESS
observations. As TESS observes a sector for only 27 days, if a binary were to
have a much longer orbital period, there would be a much lower probability of
observing it at the right time. The four KMM18 self-lensing systems have
orbital periods of $>$ 400 days, and three of them were not observed at times
when their pulses were anticipated. One system (KIC 12254688), however, was
not only observed at the right time, but rather fortuitously, TESS detected
two pulses in Sectors 41 and 56, respectively. These pulses mark the first two
self-lensing events reported from TESS archival data, despite the fact that
KIC 12254688 exhibits by far the smallest of the four self-lensing flares in
the Kepler data.

Furthermore, after the end of the primary mission of TESS, the longest-
exposure Full Frame Images (FFIs) have cadences of 10 minutes compared to the
30-minute cadence of both the TESS primary mission and Kepler. In the cases of
Sectors 41 and 56, the cadence was even lower (see Section 2), implying a
signal-to-noise ratio (SNR) much lower than that achieved in the original
Kepler detection of the KIC 12254688 pulse (KMM18). While the low SNR
complicated the extraction of the TESS pulses, we were still able to identify
them and measure them in the noisy light curves of both sectors.

In Section 2, we detail the methods by which we acquired, detrended, binned,
and folded the two light curves from Sectors 41 and 56. In Section 3, we
describe the methodology and the results of fitting our self-lensing model to
the folded TESS pulse of KIC 12254688 before comparing our best-fit model
parameters to those obtained by KMM18 and to those our model produced from the
original Kepler data set. In Section 4, we discuss the implications of this
type of modeling, including the difficulties in discovering such lensing
events in TESS data. In Section 5, we summarize this work and our conclusions.

##  2 Acquiring and Detrending the Pulses

Light curves for Sectors 41 and 56 were acquired using the Python-based time
series analysis package Lightkurve (Cardoso et al., 2018). Specifically, we
used Lightkurve’s integrated TESScut (Brasseur et al., 2019) function
search_tesscut to get the data products for both sectors. The individual frame
exposure times for Sectors 41 and 56 were $475$ and $158$ seconds,
respectively.

We used different methods to remove TESS systematics from the two light
curves. For sector 41, we used Lightkurve’s integrated Cotrending Basis
Vectors function known as CBVCorrector, which removes systematic noise. Due to
unknown errors from CBVCorrector when applied to Sector 56, we relied instead
on Lightkurve’s RegressionCorrector function in this sector.
RegressionCorrector removes scattered light background and spacecraft motion
from TESS light curves using linear regression with a design matrix
constructed from the regressors.

Figure 1: Detrended light curves of KIC 12254688 from TESS Sectors 41 and 56.
The projected pulse midpoints are expected to lie at 2424.734 BTJD and
2843.449 BTJD in Sectors 41 and 56, respectively.

The detrended light curves are shown in Figure 1. One may notice that the
pulses are not inherently visible in these detrended light curves
(particularly in Sector 56), as the expected magnification from the lensing
events (about 1.0005 in relative flux) is completely washed out by the noise
intrinsic to each light curve.

Figure 2: Binned light curves of KIC 12254688 for TESS Sectors 41 and 56,
respectively. We have zoomed in to the locations of the pulses (near the
middle of each plot) so that they are more easily noticeable. The midpoint of
the Sector 41 pulse is at 2424.734 BTJD and Sector 56 is at 2843.449 BTJD.
Sector 56 has a small downward trend to the right of the pulse that we
subsequently remove using a second-order polynomial fit, in which we exclude
the pulse itself.

In hopes of extracting the pulses from these light curves, we binned each data
set with the number of bins chosen to best reveal the pulses. For Sector 41,
this ended up being 11 data points per bin, while for Sector 56, we assigned
45 points per bin. The error bars were calculated using the error of the mean
in each bin. The resulting light curves are shown in Figure 2. At times
2424.734 BTJD and 2843.449 BTJD pulse shape profiles now appear in the data.
We note that we used only the second half of Sector 56 data, in which the
pulse was projected to appear. A downward trend in seen throughout this light
curve. As this trend is not explained by self-lensing or any other relevant
effect (e.g., Doppler-boosting kinematics or ellipsoidal variations due to
tidal forces), we removed it by applying a second-order polynomial fit to the
data on the sides of the pulse, after excluding the pulse itself. This
flattened the light curve on the sides, leaving only the pulse shape centered
at 2843.449 BTJD to rise above the flat background.

We normalized the light curves about the means of their respective data sets,
excluding the pulses in both sets. We then folded the two light curves using
the orbital period from KMM18 ($418.715$ days), combining the two pulses into
one better-sampled pulse. Still, this combined TESS pulse shows much more
variation than the original Kepler pulse due to the worse SNR in the TESS
light curves. This fact inexorably leads to the question of whether a self-
lensing model fitted to these TESS data would return best-fit parameters
similar to those determined by KMM18 from the Kepler light curve of KIC
12254688.

##  3 Fitting the Folded Pulse

As described in Sorabella et al. (2022b), we use our post-Newtonian orbital
motion code derived from the Paczynski-Wiita potential to find the velocities
and binary separation of the two stars in time. In this modeling effort, we
exclude additional effects analyzed in Sorabella et al. (2022a), such as
relativistic Doppler boosting and ellipsoidal variations. While these effects
are expected in some stellar-compact binary systems, they become significant
only in systems with shorter orbital periods (smaller separations) than that
of KIC 1225488 (Masuda & Hotokezaka, 2019).

###  3.1 Basic Equations

The magnification caused by a self-lensing event (e.g., Paczyński, 1986; Witt
& Mao, 1994; Marsh, 2001; Lee et al., 2009; D’Orazio & Di Stefano, 2018)
depends mainly on (a) $\rho_{\star}(t)$, the time-dependent ratio between the
radius of the star and the Einstein radius of the lens $R_{\rm E}(t)$, and (b)
$u(t)$, the time-dependent projected separation between the lens and the
source, expressed in units of the Einstein radius. The ratio $\rho_{\star}$ is
given by

$\rho_{\star}(t)={\frac{R_{\star}}{R_{\rm E}(t)}}\,,$ (1)

where $R_{\star}$ is the radius of the source, and the Einstein radius is
given by

$R_{\rm E}(t)=\sqrt{\frac{4GM}{c^{2}}d(t)\sin(i)\,}\,,$ (2)

where $G$ is the gravitational constant, $M$ is the mass of the compact
object, $c$ is the speed of light, $d(t)$ is the time-dependent separation
between the binary components, and $i$ is the inclination of the binary orbit
to the plane of the sky (note that $i\approx 90^{\circ}$ for self-lensing
binary systems).

The projected separation in units of $R_{\rm E}$ is determined as a function
of time (Rahvar et al., 2011) by

$u(t)={\sqrt{u_{0}^{2}+\left({{\frac{t-t_{0}}{{t_{\rm E}}}}}\right)^{2}}}\,,$
(3)

where $t_{0}$ is the time of closest projected approach and

$u_{\rm 0}=\frac{d(t_{0})}{R_{\rm E}(t_{0})}\sin{\psi}\,$ (4)

is the impact parameter in units of $R_{\rm E}$ (Wiktorowicz et al., 2021) at
$t=t_{0}$, and $\psi$ is the inclination of the source relative to the
observer-lens line-of-sight (so the orbital inclination is
$i=90^{\circ}-\psi$). The Einstein time $t_{\rm E}$ represents the time it
takes for the lens to move across the length of its Einstein radius $R_{\rm
E}$.

Finally, the time-dependent magnification $A(t)$ is expressed in terms of
$u(t)$ and $\rho_{\star}$ by the equation (Hamolli et al., 2015)

$A=\frac{\displaystyle\int_{0}^{2\pi}\int_{0}^{\rho_{\star}}\left(\frac{\xi^{2}+2}{\xi\,\sqrt{\xi^{2}+4}}\right)S\,rdrd\theta}{2\pi\displaystyle\int_{0}^{\rho_{\star}}S\,rdr}\,,$
(5)

where polar coordinates $(r,\theta)$ are set on the projected surface of the
star, and

$\xi^{2}=u^{2}+r^{2}-2ur\cos{\theta}\,.$ (6)

The function $S(r)$ represents the brightness profile of the source, for which
we adopted the linear limb-darkening relation of Hamolli et al. (2015) and Han
(2016), viz.

$S(r)=1-\Gamma\left[1-\frac{3}{2}\sqrt{1-\left(\frac{r}{\rho_{\star}}\right)^{2}}\,\right]\,,$
(7)

where $\Gamma$ is the linear limb-darkening coefficient.

###  3.2 EMCEE Modeling

To fit our model to the data set, we use EMCEE (Foreman-Mackey et al., 2013),
the Python implementation of the Affine Invariant Markov chain Monte Carlo
Ensemble sampler (Goodman & Weare, 2010). The EMCEE code accepts an input
parameter set $\\{x_{i}\\}$ and maximizes the log-likelihood function

$\ln
P(y|\\{x_{i}\\},\sigma_{n})=-\frac{1}{2}\sum_{n}\left[\frac{(y_{n}-m_{n})^{2}}{\sigma_{n}^{~{}2}}+\ln\sigma_{n}^{~{}2}\right],$
(8)

where $y_{n}$ are the light-curve data points, $\sigma_{n}$ are the summations
of the photometric error at each point and additional jitter in the data that
the model doesn’t account for ($f_{LC}$), and $m_{n}$ are the model’s points
given set $\\{x_{i}\\}$. This is accomplished by the Python module
scipy.optimize (Van Rossum & Drake, 1995) which returns the most likely
$\\{x_{i}\\}$ values, but not their uncertainties. These are obtained from the
posterior probability distribution which is derived from
$P(y|\\{x_{i}\\},\sigma_{n})$ and a predetermined prior function $P(m_{n})$,
viz.

$P(m_{n}|\\{x_{i}\\},y_{n},\sigma_{n})\propto
P(m_{n})\,P(y|\\{x_{i}\\},\sigma_{n})\,.$ (9)

The prior allows one to encode any prior information concerning the free
parameters, and thus to preclude nonphysical results.

Figure 3: Our joint fits to the phase folded light curves of KIC 12254688 from
TESS (top left) and Kepler (bottom left) with updated radial velocity data
(right panels) from Masuda et al. (2020). The orange lines represent the 95%
confidence intervals, and the dark gray lines represent the $50^{\rm th}$
percentile estimates.

In order to compare our results with those obtained by KMM18, we also fit our
model to the phase-folded self-lensing pulses originally observed by the
Kepler Space Telescope. As fitting to self-lensing pulses on their own can
prove difficult given the degeneracies between parameters such as white dwarf
mass and stellar radius, we decided to do a joint fit of the self-lensing data
and radial velocity data of the system in similar fashion to KMM18. The radial
velocity (RV) data is as reported in Masuda et al. (2020), which was collected
from the Tillinghast Reflector Echelle Spectrograph (TRES), while the RV model
we use is calculated from the projected 2D velocity components of the star as
found by our post-Newtonian orbital motion code. This data set is an updated
version of the one found in KMM18 and covers nearly a full orbit of the
system. It is also worth noting that these are relative radial velocities,
where the motion of the binary system’s barycenter is compared relative to
another in the sample (whose RV is set to zero). As such, we included a radial
velocity offset term $RV_{0}$, to shift the RV curve to account for this
relative motion.

As the initial conditions of our orbital motion code set the binary components
at apoapsis, an additional parameter we dubbed the time offset ($t_{i}$) is
included in our modeling to shift their orbital position in time to match the
starting point of the radial velocity and light curve data. This joint fit
method was done for both the Kepler and TESS light curves and was accomplished
by maximizing the summation of the likelihood functions of the self-lensing
and radial velocity models. As it is a combination of two models, we added an
additional jitter term, specific to the radial velocity model ($f_{RV}$).

###  3.3 Priors

In the self-lensing models of KIC 12254688, the prior function allows us to
constrain EMCEE by the known properties of the free parameters adopted in this
study: (a) ranges for the mass and the radius of the primary star based on
spectroscopic results (KIC DR25, Mathur et al., 2017), (b) a plausible range
for the mass of the white dwarf companion (Figures 1, 3, 6 in Kepler et al.,
2007), (c) a range for the orbital inclination (Figure 8 in Sorabella et al.,
2022a), and (d) the range for the limb-darkening coefficient, which is set to
(0, 1) by its definition. The priors for the RV model are based closely on
those used by KMM18 while the parameters not found in the KMM18 (or were
derived parameters) were given, within in reason, wide priors.

###  3.4 Best-Fit Models

Our model fits to the TESS and Kepler pulses and their corresponding RV fits
are show in Figure 3, while the best-fit parameters are listed in Table 1. The
corresponding corner plots of the posterior distributions are shown in
Appendix A. Following standard practices, parameter values and their 1$\sigma$
errors were obtained from the 16, 50, and 84 percentiles of the posterior
distributions. The thickness of each orange curve in Figure 3 represents the
$95\%$ confidence interval, as this was determined from $95\%$ of the walker
chains chosen randomly from the set of all walkers. The $50^{th}$ percentile
estimates are shown as the dark gray line.

While our model possesses some differences with that of the KMM18 model, our
best-fit parameter values from the Kepler data set are, in general, very
consistent with the original results of KMM18 (last column in Table 1). The
parameters that showed a significant difference ($>1\sigma$) compared with the
KMM18 results were that of the radial velocity offset and the light curve
jitter. The difference in the RV offset is undoubtedly a consequence of using
the updated data from Masuda et al. (2020), which reported
$RV_{0}=7.49_{-0.3}^{+0.2}$, very much in line with our estimation. As for the
light curve jitter, this could easily stem from a difference in methodology in
phase folding and binning the Kepler data. In both cases, the estimated jitter
in the light curve is very small (about a $0.01\%$ relative flux).

Table 1: Best-fit self-lensing+RV models for KIC 12254688 utilizing the listed
priors and the final 4000 steps of 100 independent walker chains in each data
set. The returned values are obtained from the 16, 50 and 84 percentiles of
their posterior distributions. The best-fit values from Table 4 of Kawahara et
al. (2018) (KMM18) are also listed in the last column for comparison.

| TESS Data | Adopted | Kepler Data |  | Kepler Data |
---|---|---|---|---|---|---
Parameter Name (Unit) | Our Best Fit | Prior | Our Best Fit |  | KMM18 Fit |
White Dwarf Mass ($M_{\odot}$) | $0.60^{+0.04}_{-0.05}$ | (0.1, 1.44) | $0.61^{+0.03}_{-0.04}$ |  | $0.62^{+0.09}_{-0.06}$ |
Primary Star Mass ($M_{\odot}$) | $1.52^{+0.19}_{-0.21}$ | (1.2, 1.8) | $1.57^{+0.15}_{-0.17}$ |  | $1.5^{+0.1}_{-0.1}$ |
Stellar Radius ($R_{\odot}$) | $2.57^{+0.32}_{-0.25}$ | (1.5, 3.5) | $2.16^{+0.08}_{-0.08}$ |  | $2.2^{+0.2}_{-0.1}$ |
Orbit Inclination (deg) | $89.78^{+0.13}_{-0.12}$ | (85, 90) | $89.94^{+0.04}_{-0.04}$ |  | $\left[89.95_{-0.04}^{+0.05}\right]$ | (a)
Limb Darkening Coefficient | $0.50^{+0.33}_{-0.33}$ | (0, 1) | $0.39^{+0.11}_{-0.09}$ |  | $\left[0.50_{-0.35}^{+0.34}\right]$ | (b)
Eccentricity | $0.19^{+0.01}_{-0.01}$ | (0, 1) | $0.18^{+0.01}_{-0.01}$ |  | $0.20^{+0.04}_{-0.04}$ |
Argument of Periapsis (deg) | $129.27^{+1.26}_{-1.34}$ | (100, 150) | $129.29^{+1.29}_{-1.39}$ |  | $129^{+2.00}_{-2.00}$ |
Radial Velocity Offset, RV0 (km s-1) | $-7.45^{+0.05}_{-0.05}$ | ($-15,\,15$) | $-7.45^{+0.05}_{-0.05}$ |  | $-7.7^{+0.4}_{-0.5}$ |
Time Offset, $T_{i}$ (days) | $175.54^{+1.69}_{-1.60}$ | (100, 300) | $175.97^{+1.78}_{-1.63}$ |  | $\cdots$ |
Light Curve Jitter, $\log(f_{LC})$ | $-4.93^{+0.74}_{-0.72}$ | ($-6,\,-3$) | $-4.36^{+0.08}_{-0.08}$ |  | $-3.96^{+0.04}_{-0.04}$ |
RV Jitter, $\log(f_{RV})$ | $-2.63^{+0.31}_{-0.25}$ | ($-3,\,3$) | $-2.63^{+0.30}_{-0.25}$ |  | $-1.92^{+0.70}_{-0.78}$ |

${(a)}$ Inclination $i$ calculated from equation (A1). ${(b)}$ Linear
coefficient $u_{1}$ calculated from the $q$-coefficients in Table 4 of KMM18
($q_{1}=0.7^{+0.2}_{-0.3}$ , $q_{2}=0.3\pm 0.2$) using
$u_{1}=2q_{2}\sqrt{q_{1}}$ (Kipping, 2013) and error propagation.

The results for the TESS fit were largely the same as the Kepler fit, though
with greater uncertainties for most parameters. However, there were a couple
of parameter estimates that showed significant difference compared to the
KMM18 results (light curve jitter and radial velocity offset) and two more
that showed significant difference with both KMM18 and that of our own Kepler
fit (orbital inclination and stellar radius). The difference with the radial
velocity offset certainly stems from the same reason as mentioned for our
Kepler fit. The light curve jitter term actually came out to be smaller than
even the Kepler fit, but this is likely a result of the already much larger
photometric error present in the TESS light curve.

The difference between the Kepler and TESS orbital inclination estimates is
significant though both place the system as nearly edge-on as expected. As for
stellar radius, the estimated value from the TESS fit is much higher than both
the KMM18 result and our own Kepler estimate, even to the point that the error
bars did not overlap with either result. It is also worth noting that the
uncertainty for $\Gamma$ was much higher than in the Kepler fit as well. Both
the uncertainty in $\Gamma$ and the higher value of $R_{\star}$ are certainly
consequences of the significantly less defined nature of the TESS pulse
compared with its Kepler counterpart. Without a definitive shape, it is nearly
impossible to determine the limb-darkening coefficient which ends up driving
the stellar radius to higher values.

While it is clear visually in this TESS data that there is a relative flux
increase of similar magnitude and timescale to that of the Kepler pulses, the
higher intrinsic variability of the TESS light curves makes parameter
estimation more uncertain. While not diverting too far off from the estimates
made by both our and KMM18’s fit to the Kepler data, the larger size of the
error bars and less Gaussian shape of the posterior distributions (see
Appendix A) would inexorably lead to the conclusion that the estimates from
the Kepler data are more reliable. Furthermore, fitting to the RV data without
joint fitting with the TESS data led to markedly similar results as to what we
saw from our TESS joint fit. This would mean much of the parameter
determination with the TESS model stems from the quality of the RV data and
not the light curve. This coupled with the fact that the parameters that are
only determined by the self-lensing model (stellar radius and limb-darkening
coefficient) show significant differences with either Kepler model, is a sign
that the TESS pulse does not have the required signal-to-noise to accurately
estimate the binary parameters. Regardless, the presence of TESS pulses at
exactly the predicted times, with similar timescale and height to that of the
original Kepler events, provides confidence that these are indeed the
anticipated self-lensing pulses, making these signals the first self-lensing
events reported from TESS observations.

##  4 Variability, Synthetic Light Curves, and Discussion

Figure 4: Synthetic light curves of a self-lensing binary system containing a
$1.1M_{\odot}$ star with a radius of $1.5R_{\odot}$ and a $0.35M_{\odot}$ WD
(left panels) or a $7M_{\odot}$ BH (right panels), with orbital periods of
$270$ and $9$ days, respectively. The top panels simulate the light curves as
viewed by the $475$-s ($\sim 8$ min) cadence TESScut exposures, and the bottom
panels are the same light curves as viewed by the 30-min cadence of Kepler.

It is important to note the challenges in detecting self-lensing pulses, such
as those reported here, in TESS data, which is far noisier than Kepler at the
same apparent magnitude. Exposure times for TESS have been shortened in the
extended mission from 2020 through the present. This is a decision which may
merit re-examination for future cycles to restore sensitivity to small-
amplitude pulses. KIC 12254688 was a known system, so we knew exactly when
each pulse would occur; even then, further analysis had to be done to extract
the pulses from the light curves. Had this been an unknown system, it would
have been difficult to discern these pulses from intrinsic variability in the
detrended light curves. This highlights an important limitation for lensing in
TESS data as compared to Kepler data. In the Kepler observations where KMM18
detected lensing events, the pulses stuck out visually in the light curves.
This may not be the case in TESS observations since the other known self-
lensing WD systems exhibit only slightly higher pulses than KIC 12254688.
Thus, despite the projections for many more self-lensing systems containing
WDs rather than NSs or BHs (e.g., Wiktorowicz et al., 2021), their smaller
magnifications make WD lenses harder to detect in TESS data sets.

White dwarf lensing events also have the disadvantage of requiring longer
orbital periods (larger separations) to be discernible, so TESS is likely to
see just one pulse for most targets. Even if such a pulse were observed
clearly in a sector, it would still be difficult to confirm the detection as
an actual lensing event. To underscore the difficulties with WD lenses in
(sub)10-min cadence TESS data, we generated synthetic light curves for a
variety of cases (Figure 4). In all four panels, the source is a
$1.1M_{\odot}$ star with radius $1.5R_{\odot}$ and the same TESS magnitude as
KIC 12254688 ($T=12.67$), in an eccentric orbit ($e=0.5$ and
$\omega=270^{\circ}$). In the left panels, the compact object is a
$0.35M_{\odot}$ WD with an orbital period of $270$ days. The $\sim$ 0.2%
magnification due to self-lensing is about the average of all currently known
WD lensing magnifications. In the right panels, the compact object is a
$7M_{\odot}$ BH with an orbital period of only 9 days. The top panels show the
signals recorded with a $475$-s cadence, the same cadence as for TESS Sector
41; the bottom panels show the same signals recorded with the SNR seen in the
30-min cadence of Kepler.

While the pulse sticks out in the TESS WD light curve, it is hardly as
prominent as the Kepler WD pulse. In the absence of another detection, it
would be difficult to argue for the origin of the TESS pulse. Since the search
for self-lensing pulses inevitably requires searching through thousands of
light curves, it may be easy for researchers to overlook such signals without
doing additional data processing. Additional processing requires more time,
time which could be spent looking for obvious signals in other data sets.
Furthermore, the data gaps in the middle of each TESS light curve last for
days, further reducing the observing duty-cycle for each target.

In the BH panels of Figure 4, there are three self-lensing pulses, each with a
much higher magnification ($\sim$ 0.4%), due to the higher mass of the compact
object. Given the much shorter orbital period, additional physical effects,
such as Doppler boosting and ellipsoidal variations, appear in the light
curves. Including these effects in the models further constrains the free
parameters. For instance, the new effects are highly dependent on the
eccentricity of the orbit, a parameter that the pure lensing model cannot
model on its own with great certainty. Thus, for BH lenses the orbital
parameters can be tightly constrained even without follow-up spectroscopic
measurements, in contrast to our present analysis of the KMM18 WD lens. BH
lenses are expected to be extremely rare, but if detected their light curves
will provide unambiguous self-lensing signals despite the increased noise in
TESS photometry relative to Kepler.

##  5 Summary and Conclusions

We have reported observations of the first two self-lensing pulses by the TESS
Space Telescope. These pulses occurred in the known self-lensing binary KIC
12254688, where a WD gravitationally lenses its stellar companion once per
orbit. We detrended the light curves from TESS Sectors 41 and 56 (Figure 1),
binned the data (Figure 2), and folded them to provide a better sampled pulse.
We joint fitted our self-lensing and radial velocity models to both the folded
TESS pulse and the original folded Kepler pulse and corresponding radial
velocity data from Masuda et al. (2020) (Figure 3). We obtained best-fit
parameters (Table 1) that were, overall, consistent with those determined by
the original study of Kawahara et al. (2018) of the Kepler data. However, we
found inconsistencies in some of the returned TESS parameter estimates
specific to the self-lensing model that are likely caused by the more chaotic,
less defined nature of the folded TESS pulse. These parameters are still
within reason, but are significantly different than what we and KMM18 found
for the Kepler pulses.

It has been suggested by Masuda & Hotokezaka (2019) and Wiktorowicz et al.
(2021) that we are going to discover many more Galactic self-lensing systems
containing WDs than NSs or BHs in TESS data and this could very well be the
case. On the other hand, we have highlighted the difficulty in parameter
estimation in such lower magnification systems in the lower-cadence extended
TESS mission as well as the difficulties that we will face in even
distinguishing such WD self-lensing pulses from intrinsic variability (§ 4).
Thus, it may be wise to first search the higher-cadence TESS data obtained
during the primary mission. We also note the possibility of detecting a BH
lensing event by mere chance. In such a case, it would be much easier to
identify its repeated short-period pulses in the pipeline data, especially in
the lower-cadence TESS data of the extended mission.

We would like to thank Kento Masuda and the members of the TRES team for
providing the radial velocity measurements of KIC 12254688, and the members of
the Lowell Center for Space Science and Technology (LoCSST) for their support
throughout this investigation. Detailed reviewer comments and suggestions are
also greatly appreciated. This work was supported by a NASA FINESST grant, NSF
grant AST-2109004, and the Massachusetts Space Grant Consortium. The data
presented in this paper were obtained from the Mikulski Archive for Space
Telescopes (MAST). The particular sets of observations that we analyzed can be
downloaded from link10.17909/hm87-b313 (catalog 10.17909/hm87-b313).

## References

  * Brasseur et al. (2019) Brasseur, C. E., Phillip, C., Fleming, S. W., et al. 2019, Astrocut: Tools for creating cutouts of TESS images, ASCL record ascl:1905.007.  http://ascl.net/1905.007
  * Cardoso et al. (2018) Cardoso, J. V. d. M., Hedges, C., Gully-Santiago, M., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python, Astrophysics Source Code Library.  http://ascl.net/1812.013
  * D’Orazio & Di Stefano (2018) D’Orazio, D.-J., & Di Stefano, R. 2018, MNRAS, 474, 2975, doi: 10.1093/mnras/stx2936
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65, doi: 10.2140/camcos.2010.5.65
  * Hamolli et al. (2015) Hamolli, L., Hafizi, M., De Paolis, F., & Nucita, A. A. 2015, Advances in Astronomy, 2015, 402303, doi: 10.1155/2015/402303
  * Han (2016) Han, C. 2016, ApJ, 820, 53, doi: 10.3847/0004-637x/820/1/53
  * Hogg & Foreman-Mackey (2018) Hogg, D. W., & Foreman-Mackey, D. 2018, ApJS, 236, 11, doi: 10.3847/1538-4365/aab76e
  * Kawahara et al. (2018) Kawahara, H., Masuda, K., MacLeod, M., et al. 2018, AJ, 155, 144 (KMM18), doi: 10.3847/1538-3881/aaaaaf
  * Kepler et al. (2007) Kepler, S. O., Kleinman, S. J., Nitta, A., et al. 2007, MNRAS, 375, 1315, doi: 10.1111/j.1365-2966.2006.11388.x
  * Kipping (2013) Kipping, D. M. 2013, MNRAS, 435, 2152, doi: 10.1093/mnras/stt1435
  * Kruse & Agol (2014) Kruse, E., & Agol, E. 2014, Science, 344, 275, doi: 10.1126/science.1251999
  * Lee et al. (2009) Lee, C.-H., Riffeser, A., Seitz, S., & Bender, R. 2009, ApJ, 695, 200, doi: 10.1088/0004-637X/695/1/200
  * Maeder (1973) Maeder, A. 1973, A&A, 26, 215
  * Marsh (2001) Marsh, T. R. 2001, MNRAS, 324, 547, doi: 10.1046/j.1365-8711.2001.04293.x
  * Masuda & Hotokezaka (2019) Masuda, K., & Hotokezaka, K. 2019, ApJ, 883, 169, doi: 10.3847/1538-4357/ab3a4f
  * Masuda et al. (2019) Masuda, K., Kawahara, H., Latham, D. W., et al. 2019, ApJ, 881, L3, doi: 10.3847/2041-8213/ab321b
  * Masuda et al. (2020) Masuda, K., Kawahara, H., Latham, D. W., et al. 2020, in White Dwarfs as Probes of Fundamental Physics: Tracers of Planetary, Stellar and Galactic Evolution, ed. M. A. Barstow, S. J. Kleinman, J. L. Provencal, & L. Ferrario, Vol. 357, 215–219, doi: 10.1017/S1743921320000915
  * Mathur et al. (2017) Mathur, S., Huber, D., Batalha, N. M., et al. 2017, ApJS, 229, 30, doi: 10.3847/1538-4365/229/2/30
  * Paczyński (1986) Paczyński, B. 1986, ApJ, 304, 1
  * Rahvar et al. (2011) Rahvar, S., Mehrabi, A., & Dominik, M. 2011, MNRAS, 410, 912, doi: 10.1111/j.1365-2966.2010.17490.x
  * Ricker et al. (2010) Ricker, G. R., Latham, D. W., Vanderspek, R. K., et al. 2010, in American Astronomical Society Meeting Abstracts, Vol. 215, American Astronomical Society Meeting Abstracts #215, 450.06
  * Sorabella et al. (2022a) Sorabella, N. M., Bhattacharya, S., Laycock, S. G. T., Christodoulou, D. M., & Massarotti, A. 2022a, ApJ, 936, 63, doi: 10.3847/1538-4357/ac82b7
  * Sorabella et al. (2022b) —. 2022b, ApJ, 927, 234, doi: 10.3847/1538-4357/ac4a59
  * Sorabella et al. (2023) Sorabella, N. M., Laycock, S. G. T., Neeley, L. J., Christodoulou, D. M., & Bhattacharya, S. 2023, ApJ, 954, 59, doi: 10.3847/1538-4357/ace9df
  * Van Rossum & Drake (1995) Van Rossum, G., & Drake, Jr, F. L. 1995, Python Reference Manual (CWI, Amsterdam).  https://ir.cwi.nl/pub/5008
  * Wiktorowicz et al. (2021) Wiktorowicz, G., Middleton, M., Khan, N., et al. 2021, MNRAS, 507, 374, doi: 10.1093/mnras/stab2135
  * Witt & Mao (1994) Witt, H. J., & Mao, S. 1994, ApJ, 430, 505, doi: 10.1086/174426

##  Appendix A MCMC Statistics

###  A.1 Corner Plots

Corner plots were generated by the Python module corner.py (Foreman-Mackey,
2016; Hogg & Foreman-Mackey, 2018) from the derived posterior distributions of
the free parameters of our models. Figures 5 and 6 display the corner plots
from the TESS and the Kepler data fits, respectively.

In general, the posterior distributions of the TESS fit were much wider and
less Gaussian than that of their Kepler counterparts. The stellar mass (‘m2’
in columns 2) and the limb-darkening coefficient $\Gamma$ were not constrained
particularly well, which had consequences on driving $R_{\star}$ and the
inclination $\psi$ to higher values than what we see in our Kepler estimates.
This is largely due to the less defined shape of the TESS pulse compared with
its Kepler counterpart.

###  A.2 Orbit Inclination and Impact Parameter

KMM18 does not report an orbital inclination in their fitted parameters.
Instead, they fit to a normalized impact parameter denoted by $b\in(0,1)$.
Expressed in terms of our notation, we find that this is equivalent to

$\frac{u_{0}}{(u_{0})_{\rm max}}=\frac{\cos(i)}{\cos(i_{\rm min})}=0.09\,,$
(A1)

where $i_{\rm min}$ is the minimum required orbital inclination for a self-
lensing event at mid-eclipse as found in Wiktorowicz et al. (2021) given by

$i_{\rm min}=\cos^{-1}\left(\frac{R_{\rm
E}(t_{0})+R_{\star}}{d(t_{0})}\\!\right)\,.$ (A2)

For KMM18, $i_{min}\approx 89.49^{\circ}$ which, by plugging into equation
(A1), allows us to find the orbital inclination of $i=89.95^{\circ}$. We then
use error propagation to find the error bars for this value as reported in
Table 1.

Figure 5: One and two dimensional projections of the posterior probability
distributions of the model parameters for the TESS data fit. The returned
values are obtained from the 16, 50 and 84 percentiles, utilizing the final
4000 steps of 100 independent MCMC walker chains. Figure 6: As in Figure 5,
but for the Kepler data fit.


# An Optimization Approach for a Robust and Flexible Control in Collaborative
Applications

Federico Benzi and Cristian Secchi F. Benzi and C. Secchi are with the
Department of Sciences and Methods of Engineering, University of Modena and
Reggio Emilia, Italy<EMAIL_ADDRESS>
###### Abstract

In Human-Robot Collaboration, the robot operates in a highly dynamic
environment. Thus, it is pivotal to guarantee the robust stability of the
system during the interaction but also a high flexibility of the robot
behavior in order to ensure safety and reactivity to the variable conditions
of the collaborative scenario.
In this paper we propose a control architecture capable of maximizing the
flexibility of the robot while guaranteeing a stable behavior when physically
interacting with the environment. This is achieved by combining an energy tank
based variable admittance architecture with control barrier functions. The
proposed architecture is experimentally validated on a collaborative robot.

##  I Introduction

In order to implement an effective human robot collaboration, it is
fundamental to guarantee that the robot can adapt online its behavior in a
stable and efficient way. A lot of work in this direction has been done using
different techniques such as task redundancy [1, 2], adaptive physical
interaction [3, 4] and coaching [5, 6] to name a few. Nevertheless, at the
best of the authors’ knowledge, a control architecture that allows to address
the problems of robust stability and flexibility both in free motion and
during the interaction is still missing.

The behavior of a robot can be defined by specifying the way it physically
interacts with the environment and the set of (dynamic) constraints it has to
satisfy. Adapting the interaction implies to dynamically vary the parameters
(e.g. stiffness, inertia) that determine the way forces are exchanged with the
environment. Dynamically constraining the motion of the robot allows to
dynamically shape its behavior.

Admittance control [7] is a very popular strategy for controlling the
interaction, and it has been widely used in collaborative scenarios (see,
e.g., [4, 3, 8]). It enforces the robot to reproduce a desired passive
physical behavior (e.g. mass-spring-damper), i.e. the admittance dynamics.
Because of the passivity of the admittance dynamics, the behavior of the robot
is robustly stable, i.e. stable both in free motion and when interacting with
a poorly known environment (see e.g.[9]).

When using admittance control, the interaction dynamics is determined by the
choice of the parameters of the admittance dynamics (e.g. inertia and
stiffness) and a specific set of parameters may not be suitable for all stages
of a given application [4].

This severely limits the flexibility of the controller, as changing online the
dynamic parameters can lead to the loss of passivity [10] and, consequently,
to a possibly unstable behavior in human robot collaboration as shown in [3].

Several strategies for implementing a variable admittance control have been
proposed (see e.g. [3, 11]) but they allow only a limited variation of the
interaction dynamics while significantly complicating the overall control
architecture

In [12], an optimization framework has been proposed in order to implement a
variable admittance controller in a robust and flexible way. Here, energy
tanks [10, 13] have been exploited for disembodying passivity from a specific
physical behavior. The energy stored in the tank has been exploited as a
constraint to be satisfied for passively implementing any desired admittance
dynamics. This allows to achieve the maximal flexibility while satisfying the
passivity constraints, i.e. robust stability. Yet [12] does not differentiate
between the interaction phase and the free motion of the robot, and the
admittance control is active in both situations. In a collaborative scenario,
the robot often switches between free motion and interaction and constantly
requiring the satisfaction of both the passivity constraint and the continuous
implementation of a specific admittance dynamics may affect the performance of
the robot in free motion.

Moreover, robust stability is not enough to guarantee the complete safety and
the flexibility of a robot in a shared workspace: additional time-varying
constraints, such as obstacle and human avoidance, self collision and bounding
the joint position, need to be satisfied for ensuring a safe and flexible
behavior. Control barrier functions (CBFs) [14, 15] have been successfully
exploited in many robotics applications for forcing the robot to remain in a
desired subset of the state space, i.e. to dynamically constrain the behavior
of the robot (see e.g. [16, 17, 18]). Time Varying CBFs [19] allow to consider
time-varying constraints and they have been successfully exploited for
enforcing safety in human robot collaboration [20, 21]. Furthermore, tasks to
be executed can be encoded as CBFs as shown in [22]. The control input for
satisfying the dynamic constraints modeled by the CBFs can be found by solving
a convex optimization problem, which makes the CBFs amenable to real-time
implementation.

Passivity can be encoded using Control Barrier Functions [23] but no
interaction controllers are available yet.

In this paper we propose a novel control architecture that merges energy
tanks, for implementing a flexible and passive interactive behavior, and
control barrier functions, for dynamically constraining the behavior of the
robot. The proposed architecture leverages control barrier functions to encode
a series of constraints to be satisfied and energy tanks for passively
implementing a variable admittance controller. The control input is the
solution of a single convex optimization problem that considers both the
passivity constraint and the constraints due to the control barrier functions.
Using the proposed architecture it will be possible to make a collaborative
robot extremely flexible both during free motion and during the interaction
with possibly unknown environments.

The contributions of this paper are:

  * •

A convex optimization problem merging the passivity constraint and control
barrier functions constraints;

  * •

A control architecture capable of maximizing the flexibility of the robot both
during the interaction and during free motion;

  * •

An experimental validation of the proposed architecture on a collaborative
robot.

This paper is organized as follows: in Sec. II the main problem addressed in
this paper is formulated. In Sec. III energy tanks, their usage in passivity
preservation and the overall tank-based architecture is presented. In Sec. IV
the CBF-based architecture for flexible task execution is described, while in
Sec. V the collaborative constraint-oriented control architecture is shown. In
Sec. VI the proposed architecture is validated onto a robotic setup, and
finally in Sec. VII conclusions are drawn and future work is addressed.

##  II Problem Formulation

Figure 1: The admittance control architecture. The input $F_{c}$ represent the
control force of the low level controller actuating the robot.

Consider a velocity controlled fully actuated $n$-DOFs manipulator. The model
of the robot can be described by:

$\dot{x}=J(q)u$ (1)

where $x\in\mathbb{R}^{m}$ is the pose of the end-effector and
$q\in\mathbb{R}^{n}$ is the vector of joint variables.
$J(q)\in\mathbb{R}^{m\times n}$ is the Jacobian of the robot and
$u\in\mathbb{R}^{n}$ is the joint velocity input.

Thanks to the low-level velocity controller, a desired velocity profile
$\dot{x}_{des}\in\mathbb{R}^{m}$ can be reproduced by simply setting
$u=J^{+}(q)\dot{x}_{des}$, where $J^{+}(q)$ denotes the pseudo-inverse of the
Jacobian $J(q)$, such that $\dot{x}(t)\approx\dot{x}_{des}$.

The interaction with the environment is regulated by the admittance control
architecture reported in Fig. 1. In order to enforce a desired dynamic
behavior for the robot, the interaction force $F_{e}(t)\in\mathbb{R}^{m}$ is
integrated by the admittance dynamics that produces a desired velocity
$\dot{x}_{a}$. By setting $\dot{x}_{des}(t)=\dot{x}_{a}\approx\dot{x}$, the
robot exactly reproduces the admittance dynamics. Thus, if the admittance
dynamics is passive, the controlled robot behaves as a passive system and,
therefore, it is stable both in free motion and in contact with the
environment, i.e. robustly stable.

Admittance dynamics model mechanical systems and they can be represented by
the following Euler-Lagrange model

$M(x_{a})\ddot{x}_{a}+C(x_{a},\dot{x}_{a})\dot{x}_{a}+D(x_{a})\dot{x}_{a}+\dfrac{\partial
P}{\partial x_{a}}=F_{e}$ (2)

where $x_{a}\in\mathbb{R}^{m}$ is the pose and $M(x_{a})=M^{T}(x_{a})>0$ is
the inertia matrix, $C(x_{a},\dot{x}_{a})$ represents the Coriolis effects,
$D(x_{a})\geq 0$ is a damping matrix and
$P:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a potential field active on the
system (e.g. elastic potential).
As well known (see, e.g., [9]), (2) is passive with respect to the pair
$(\dot{x}_{a}(t),F_{e}(t))$, using the following storage function:

$H_{a}(x_{a},\dot{x}_{a})=P(x_{a})+\frac{1}{2}\dot{x}_{a}^{T}M(x_{a})\dot{x}_{a}$
(3)

Unfortunately, when considering a variable admittance, i.e. such that the
physical parameters characterizing the desired dynamics are time-varying,
passivity can be lost. In fact, considering the time-varying version of (2)

$M(x_{a},t)\ddot{x}_{a}+C(x_{a},\dot{x}_{a},t)\dot{x}_{a}+D(x_{a},t)\dot{x}_{a}+\dfrac{\partial
P}{\partial x_{a}}(t)=F_{e}$ (4)

we have that (3) becomes then

$H_{a}(x_{a},\dot{x}_{a},t)=P(x_{a},t)+\frac{1}{2}\dot{x}_{a}^{T}M(x_{a},t)\dot{x}_{a}$
(5)

whence, using (4), we can write

$\dot{H}_{a}(x_{a},\dot{x}_{a},t)=F_{e}^{T}\dot{x}_{a}-\dot{x}^{T}_{a}D(x_{a},t)\dot{x}_{a}+\dfrac{\partial
H_{a}}{\partial t}$ (6)

where the last term is sign indefinite and can introduce energy into the
system and, therefore, lead to a loss of passivity.

Thus, standard admittance control is not suitable for implementing an
interaction that is flexible and passive at the same time.

Besides passivity, various dynamic and time-varying constraints are necessary
for establishing safety and other desired behaviors for the robot. These
constraint cannot be simply embedded in the standard admittance control
architecture.

We aim at developing a control architecture that allows to passively
implementing a time varying admittance dynamics as the one in (4) and to
enforce a set of dynamic and time-varying constraints on the robot.

##  III Energy Tank Architecture

In this section we show how to exploit energy tanks [10],[13] for reproducing
a variable admittance dynamics (4) on the robot (1) following the approach
presented in [12].

Energy tanks are energy storing elements which can be represented by:

$\begin{cases}\dot{x}_{t}=u_{t}\\\
y_{t}=\frac{\partial{T}}{\partial{x_{t}}}=x_{t}(t)\end{cases}$ (7)

where $x_{t}\in\mathbb{R}$ is the state of the tank, while the pair
$(u_{t},y_{t})\in\mathbb{R}\times\mathbb{R}$ represents the power port the
tank can exchange energy through and

$T(x_{t})=\frac{1}{2}x^{2}_{t}$ (8)

is the function representing the energy stored in the tank.

Energy tanks are completely disconnected from a specific physical dynamics.
Thus, the energy stored in the tank can be exploited for implementing any
desired behavior. This can be achieved by interconnecting the power port of
the tank $(u_{t},y_{t})$ with the power port $(F_{e},\dot{x}_{des})$ of the
implemented admittance dynamics as:

$\begin{cases}u_{t}(t)=A^{T}(t)F_{e}(t)\\\
\dot{x}_{des}(t)=A(t)y_{t}(t)\end{cases}$ (9)

where $A(t)\in\mathbb{R}^{n}$ is defined as

$A(t)=\frac{\gamma(t)}{x_{t}(t)}$ (10)

and $\gamma(t)\in\mathbb{R}^{n}$ is the desired value for the output
$\dot{x}_{des}(t)$, i.e. the velocity implementing the desired admittance. In
fact, if we plug (9) into (7) we get

$\begin{cases}\dot{x}_{t}=A^{T}(t)F_{e}(t)\\\
\dot{x}_{des}=A(t)y_{t}(t)=\gamma(t)\end{cases}$ (11)

which implies that any desired port behavior can be obtained by an appropriate
modulation of the energy exchanged with the tank. From (7),(8) and (11) it
follows that

$\dot{T}=u_{t}y_{t}=A^{T}(t)F_{e}y_{t}=\gamma^{T}F_{e}$ (12)

which clearly shows that the energy needed for implementing the desired
behavior is extracted from/injected into the tank.

The modulated tank (11) faces a singularity whenever $x_{t}(t)=0$ because of
the definition of the modulation factor in (11). The singularity happens when
the tank is empty and no behavior can be implemented using the energy stored
in the tank [10]. In order to avoid this condition, it is necessary to
initialize $x_{t}$ such that $T(x_{t}(0))\geq\underline{\varepsilon}>0$ and to
ensure that $T(x_{t}(t))\geq\underline{\varepsilon}$ $\forall t>0$. As shown
in [12], the following result holds:

######  Proposition 1.

If $T(x_{t})\geq\underline{\varepsilon}$ for all $t\geq 0$, then the modulated
tank (11) remains passive independently of the desired output $\gamma(t)$.

Thus, as long as the tank is not depleted, any desired port behavior can be
passively implemented by modulating the energy stored in the tank.
Furthermore, as shown in [24] for the single port case and in [25] for the
multi-port case, any passive dynamics can be reproduced without depleting the
energy tank. Thus, modulated energy tanks can be used as a generalized
admittance control since they allow to reproduce any physical dynamics.

In case we aim at reproducing a non passive admittance dynamics, i.e. the time
varying dynamics in (4), we can exploit Prop. 1 for enforcing a passive
implementation of the desired behavior. In fact, it is possible to guarantee
the passivity of the modulated tank by enforcing the following constraint:

$T(x_{t})\geq\underline{\varepsilon}\quad\forall t\geq 0$ (13)

Thus, as shown in [12], it is possible to find the best passive implementation
of any desired admittance dynamics by solving the following optimization
problem

$\displaystyle\underset{\dot{x}_{des}}{\text{minimize}}$
$\displaystyle||\dot{x}_{des}-\dot{x}_{a}||^{2}$ (14) subject to
$\displaystyle\int_{0}^{t}F_{e}^{T}(\tau)\dot{x}_{des}(\tau)d\tau\geq-T(x_{t}(0))+\underline{\varepsilon}$

The problem (14) allows to find the best passive approximation of the desired
behavior $\dot{x}_{a}$, taking into account the amount of energy stored in the
tank. The solution $\dot{x}_{des}$ is then utilized to correctly tune the
modulation matrix $A(t)$ in (11). In this way, a passive energy balance is
guaranteed even if the admittance parameters change online.
In [12], a discrete time version of the problem was introduced, which turns
the integral constraint in (14) into a linear one. The resulting optimization
problem is then convex, which makes it computationally fast and simple in
formulation, thus suitable for being executed in real time and embedded into a
larger optimization-based framework.

##  IV Control Barrier Functions based architecture

In this section we show how to exploit CBFs for enforcing multiple time-
varying dynamic constraints on a robot described by (1). The presented
framework leverages the techniques exposed in [26], [27], where both kinematic
limits and desired tasks to be executed by the robot can be represented as
time-varying constraints on the control input. Thus, it is possible to model
tasks and physical limitations as constraints to be fulfilled by the
controlled system and CBFs are exploited for formulating an optimization
problem for finding the best input that satisfies all the constraints, e.g.
that can lead to the optimal execution of all the tasks.

We consider tasks whose execution can be expressed as the minimization of a
non negative, possibly time-varying, continuously-differentiable cost function
$C:\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}$. Using the robot model
in (1) and considering as an output variable the time-varying task variable
$\sigma\in\mathbb{R}^{n}$, this can be described by the following optimization
problem:

$\displaystyle\underset{u}{\text{minimize}}$ $\displaystyle C(\sigma,t)$ (15)
subject to $\displaystyle\dot{x}=J(q)u$ $\displaystyle\sigma=k(x,t)$

The optimization problem can be reformulated and made computationally simple
using Control Barrier Functions. Let $\mathcal{C}\subset\mathbb{R}^{n}$ be the
subset where the task is considered executed, i.e. where $C(\sigma,t)=0$.
Consider the a control barrier function
$h:\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}$ defined as
$h(\sigma,t)=-C(\sigma,t)$. By construction $h$ is non negative only in the
region of satisfaction of the task, i.e. when $C(\sigma,t)=0$. Thus, enforcing
the non negativity of $h$ is equivalent to enforcing the execution of the task
$\sigma$. This can be encoded in the following convex optimization problem
[26]:

$\displaystyle\underset{\dot{q}}{\text{minimize}}$
$\displaystyle||\dot{q}||^{2}$ (16) subject to $\displaystyle\frac{\partial
h}{\partial t}+\frac{\partial h}{\partial\sigma}\frac{\partial\sigma}{\partial
x}J(q)\dot{q}+\alpha(h(\sigma,t))\geq 0$

where $\alpha(\cdot)$ is an extended class $\mathcal{K}$ function111An
extended class $\mathcal{K}$ is a function
$\phi:\mathbb{R}\rightarrow\mathbb{R}$ such that $\phi$ is strictly increasing
and $\phi(0)=0$ and where we have considered that the input of (1) is
$u=\dot{q}$. The formulation can easily be extended in order to enable the
simultaneous execution of multiple tasks. Let us consider a set of M different
tasks $T_{1},\dots,T_{M}$ which have to be executed, each encoded respectively
by the cost functions $C_{1},\dots,C_{M}$. The execution of these tasks can
then be accomplished by solving the following convex optimization problem:

$\displaystyle\underset{\dot{q},\delta}{\text{minimize}}$
$\displaystyle||\dot{q}||^{2}+l||\delta||^{2}$ (17) subject to
$\displaystyle\frac{\partial h_{m}}{\partial t}+\frac{\partial
h_{m}}{\partial\sigma}\frac{\partial\sigma}{\partial x}J(q)\dot{q}$
$\displaystyle+\alpha(h_{m}(\sigma,t))\geq-\delta_{m}\quad
m\in\\{1,\dots,M\\}$

in which $h_{m}(\sigma,t)=-C_{m}(\sigma,t)$ and
$\delta=[\delta_{1},\dots,\delta_{M}]^{T}$ is the vector of slack variables
corresponding to each constraint, while $l\geq 0$ is a scaling factor. Each
$\delta_{i}$ indicates how much the constraint corresponding to the $i-$th
task can be relaxed, in order to guarantee the feasibility of the problem even
if conflicting constraints are active at the same time.

##  V Collaborative Constraint-Oriented Control Architecture

In this section we propose a collaborative constraint-oriented control
architecture that merges the flexible energy management architecture
illustrated in Sec. V with the CBFs based architecture for task execution
shown in Sec. IV in order to obtain a control strategy that guarantees a
maximally flexible and robustly stable behavior of the robot. We will keep the
exposition of the algorithm in continuous time, but the overall discussion can
be easily extended to discrete time (see [12] in particular for the constraint
on passivity).

It is possible to enforce a passive behavior of the robot controlled using
control barrier functions by inserting the passivity constraint of (14) in
(17) and by using the output of the optimization problem for modulating the
matrix $A(t)$ in (11). This leads to:

$\displaystyle\underset{\dot{q},\delta}{\text{minimize}}$
$\displaystyle||\dot{q}||^{2}+l||\delta||^{2}$ (18) subject to
$\displaystyle\frac{\partial h_{m}}{\partial t}+\frac{\partial
h_{m}}{\partial\sigma}\frac{\partial\sigma}{\partial x}J(q)\dot{q}$
$\displaystyle+\alpha(h_{m}(\sigma,t))\geq-\delta_{m}\quad
m\in\\{1,\dots,M\\}$
$\displaystyle\int_{0}^{t}F_{e}^{T}(\tau)J(q)\dot{q}(\tau)d\tau\geq-T(x_{t}(0))+\underline{\varepsilon}$

and the tank is modulated by setting $\gamma=J(q)\dot{q}$.

Notice that, since the passivity constraint is linear, the overall
optimization problem maintains its convexity. Furthermore, the passivity
constraint in (18) does not present a dedicated slack variable to relax it.
This is due to the fact that the preservation of passivity is of top priority
as it is connected to robust stability. At the same time, this choice does not
over-constrain the system during free motion, as demonstrated in the following
proposition.

######  Proposition 2.

If $T(x_{t}(t))\geq\underline{\varepsilon}$ during the interaction, then the
passivity constraint in (18) is automatically satisfied also in free motion.

###### Proof.

Assume that free motion starts at time $\overline{t}>0$. We can then subdivide
the integral term of the passivity constraint in (18) in two parts:

$\int_{0}^{\overline{t}}F_{e}^{T}(\tau)J(q)\dot{q}(\tau)d\tau+\int_{\overline{t}}^{t}F_{e}^{T}(\tau)J(q)\dot{q}(\tau)d\tau\geq\\\
\geq-T(x_{t}(0))+\underline{\varepsilon}$ (19)

Since $F=0$ after $\overline{t}$, (19) reduces to

$\int_{0}^{\overline{t}}F_{e}^{T}(\tau)J(q)\dot{q}(\tau)d\tau\geq-T(x_{t}(0))+\underline{\varepsilon}$
(20)

which is equivalent to the constraint active during the interaction phase,
thus concluding the proof. ∎

Figure 2: The final collaborative constrained-control architecture.

Hence, the passivity constraint is automatically satisfied during free motion,
meaning that the optimizer disregards it completely during this phase,
focusing on satisfying the remaining $M$ constraints in the stack.
Following this procedure, we have therefore obtained a control architecture
which is capable not only of accomplishing concurrent constraints in a
flexible way but also to guarantee a robust behavior while interacting with a
poorly known environment, thanks to the preservation of passivity. The
developed framework acts as a sort of ”armor”, protecting the robot against
unstable behaviors during the interaction, while still allowing it to
implement all the encoded tasks at the best of its dexterity. Finally,
everything is encompassed into a single convex optimization problem, for a
fast and efficient resolution.

In the same fashion, we can insert the desired admittance into our own
optimization problem, by properly modifying the previous objective function in
(18) as follows:

$\displaystyle\underset{\dot{q},\delta}{\text{minimize}}$
$\displaystyle||\dot{q}-\dot{q}_{a}||^{2}+l||\delta||^{2}$ (21) subject to
$\displaystyle\frac{\partial h_{m}}{\partial t}+\frac{\partial
h_{m}}{\partial\sigma}\frac{\partial\sigma}{\partial x}J(q)\dot{q}$
$\displaystyle+\alpha(h_{m}(\sigma,t))\geq-\delta_{m}\quad
m\in\\{1,\dots,M\\}$
$\displaystyle\int_{0}^{t}F_{e}^{T}(\tau)J(q)\dot{q}(\tau)d\tau\geq-T(x_{t}(0))+\underline{\varepsilon}$

in which $\dot{q}_{a}=J(q)^{+}\dot{x}_{a}$ is the desired admittance expressed
in the joint space. In this way, we have encompassed the whole variable
admittance control scheme into our task execution architecture, allowing for a
flexible interaction with the environment, while still ensuring the overall
passivity via the energetic constraint, as well as satisfying the other time-
varying constraints in the stack.

The implementation of the desired admittance can be enforced at the price of
relaxing the tasks encoded by the CBFs. A possible strategy for achieving this
goal is to dynamically weight the first term in the function to optimize in
(21). A possible solution is to use the following metric matrix
$W(t)=(1+\kappa\|F_{e}(t)\|^{2})I_{n}$, where $I_{n}$ is the identity matrix
of order $n$ and $\kappa>0$ that can be used for tuning the effect of the
external force. Thus, replacing $\|\dot{q}-\dot{q}_{a}\|^{2}$ in (21) with

$\|\dot{q}-\dot{q}_{a}\|^{2}_{W}(t)=(\dot{q}-\dot{q}_{a})^{T}W(t)(\dot{q}-\dot{q}_{a})$
(22)

allows to penalize the deviation from the desired admittance dynamics the
more, the higher the interaction force is. In free motion (22) reduces to the
standard euclidean norm used in (21).

The final formulation of the optimization problem is as follows:

$\displaystyle\underset{\dot{q},\delta}{\text{minimize}}$
$\displaystyle||\dot{q}-\dot{q}_{adm}||^{2}_{W}+l||\delta||^{2}$ (23) subject
to $\displaystyle\frac{\partial h_{m}}{\partial t}+\frac{\partial
h_{m}}{\partial\sigma}\frac{\partial\sigma}{\partial x}J(q)\dot{q}$
$\displaystyle+\alpha(h_{m}(\sigma,t))\geq-\delta_{m}\quad
m\in\\{1,\dots,M\\}$
$\displaystyle\int_{0}^{t}F_{e}^{T}(\tau)J(q)\dot{q}(\tau)d\tau\geq-T(x_{t}(0))+\underline{\varepsilon}$

The overall control architecture is reported in Fig. 2, in which we defined
$\dot{x}_{opt}=J(q)\dot{q}$ as the optimal velocity in the task space.

The velocity controlled robot is wrapped by the modulated tank that is
exploited for ensuring the passivity of the controlled system. The CBF-based
optimizer receives as an input the desired admittance velocity $\dot{x}_{a}$
and, by solving (23) determines the best value $\dot{x}_{opt}$ for passively
implementing the desired admittance behavior and all the tasks that are
encoded by the CBFs. The matrix $A(t)$ is modulated by setting
$\gamma=\dot{x}_{opt}$.

##  VI Experiments

The framework proposed in this paper has been validated on a Universal Robot
10e manipulator, equipped with an on-board 6-axis force/torque (F/T) sensor.
Both the robot and the sensor run with a sampling time of 2ms.
The robot is employed to accomplish a set of tasks, including obstacle
avoidance, joint control and position control of the end effector, as well as
the satisfaction of the passivity constraint and the implementation of a
variable admittance dynamics. Each task is encoded using a specific CBF,
following the formulations exposed in Sec. IV.
The obstacle avoidance task is encoded by the following CBF:

$h_{safe}=-\xi_{s}(d^{2}-D_{min}^{2})$ (24)

in which $\xi_{s}=10$ and $D_{min}$ is the minimum distance value between the
obstacle and the tip of the end-effector, while $d=d(x,t)$ is defined as

$d=||x-x_{obs}||^{2}$ (25)

in which $x_{obs}$ is the Cartesian position of the closest obstacle.
Following the optimization-based formulation (23), for each CBF the
correspondent gain $\xi_{(.)}$ is tuned and, after simple computations, the
resulting constraint is obtained and inserted into the optimization problem.
For the sake of simplicity, we chose the function $\alpha$ to be the identity
function, such that $\alpha(h_{m}(\sigma,t))=h_{m}(\sigma,t)\,\forall
m\in\\{1,\dots,M\\}$. In the case of (24), for example, the resulting
constraint is

$2(d-D_{min})J(q)\dot{q}\geq-(h_{safe}(x,t))+\delta_{s}$ (26)

in which $\delta_{s}$ is the dedicated slack variable. The joint control task
consists in maintaining the joint variables within an upper and lower limit.
For each $i-$th joint, a dedicated CBF has been implemented as:

$h_{lim_{i}}=\xi_{l}\frac{(q^{+}_{i}-q_{i})(q_{i}-q^{-}_{i})}{(q^{+}_{i}-q^{-}_{i})}\quad
i\in\\{1,\dots,n\\}$ (27)

in which $\xi_{l}=1$, while $q^{+}_{i}$ and $q^{-}_{i}$ represent the real
joint limits of the $i$-th joint, in the joint space. Observe how ensuring
that $h_{lim_{i}}\geq 0$ is equivalent to keeping the value of
$q_{i}\in[q^{+}_{i},q^{-}_{i}]$.
Finally, the position control task is encoded through the following CBF:

$h_{pos}=-\xi_{p}||x-x_{goal}||^{2}$ (28)

in which $\xi_{p}=5$ and $x_{goal}$ is the desired Cartesian position of the
end effector.
Additionally, the desired admittance dynamics presents a time-varying
repulsive potential $P(x,t)$ which is centered in the current position of the
obstacle. The generated force $\frac{\partial P}{\partial x}(t)$ grows as the
robot approaches the obstacle, such that the human can perceive an hazardous
area during the interaction. The formulation of the resulting force is as
follows

$\frac{\partial P}{\partial
x}(t)=\begin{cases}K_{rep}\Bigl{(}\frac{1}{d(p,t)}-\frac{1}{D^{*}}\Bigr{)}\frac{p(t)-p_{obs}(t)}{d^{3}(p,t)}\,\,\,if\,d(p,t)<D^{*}\\\
0\quad otherwise\end{cases}$ (29)

in which $D^{*}>D_{min}$ is the activation distance for the potential. As seen
before, such a behavior is non-passive, and energy needs to be extracted from
the tank for preserving a passive energy balance.
A series of experiments have been conducted, in order to concurrently validate
each component of the architecture.
First, the robot is moved to a desired goal using (28). Along the way, the
operator interacts with it, moving it around in the workspace. As soon as the
interaction stops, the robot switches to free-motion and resumes its motion
towards the goal, unaffected by the previous interaction. Fig. 3 shows the
evolution of (28) over time, together with the intensity of the force acting
on the system.
Secondly, the same position control task is implemented, this time with the
presence of a time-varying obstacle moving in the workspace. As visible in
Fig. 4, as soon as the distance $d$ approaches the minimum value
$D_{min}=0.25m$, the constraint forces the robot to move away from the
obstacle, preserving the overall safety. Once the obstacle is removed, the
robot resumes its motion towards the goal.
Finally, the passivity of the interaction is validated using a variable
admittance control. The admittance model encompasses a constant mass and
damping term $M(x_{a})$ and $D(x_{a})$, as well as the time-varying potential
in (29). In particular, the desired inertia is chosen as a diagonal matrix
with elements equal to $0.75kg$ and $0.25kg$ for translations and rotations
respectively, while the damping is a diagonal matrix with constant elements
equal to $0.05\frac{Ns}{m}$ (translations) and $0.025\frac{Ns}{m}$
(rotations).
During the experiment, the operator tries to guide the robot towards the
obstacle. As the value of $d$ drops below $D^{*}=0.5m$, the operator feels the
repulsive force (29) generated by the potential. Fig. 5 shows the evolution of
the energy in the tank: as soon as the current value of $T(x_{t}(t))$
approaches $\underline{\varepsilon}=0.1$, the constraint (13) ensures that the
tank is never depleted by implementing a passive approximation of the desired
behaviour.
Additional graphs are shown in the accompanying video, in which the three
experiments are portrayed.

Figure 3: Value of the CBF $h_{pos}$ over time and norm of the force measured
by the F/T sensor. Figure 4: Value of the CBF $h_{pos}$ over time and distance
between the end effector of the robot and the obstacle. Figure 5: Evolution of
the tank energy over time.

##  VII Conclusions and Future Work

In this paper we have proposed an optimization-based framework which combines
CBFs and energy tanks for achieving a robust and flexible interaction, while
respecting a set of time-varying constraints, thus guaranteeing the complete
safety of the robotic application in a HRC scenario. We have formulated
passivity as a constraint and inserted it, along other safety-relevant
constraints, into a single convex optimization problem. Future work aims at
establishing priorities among the different constraints in an optimal fashion
and predicting the evolution of the energy in the tank.

## References

  * [1] H. Su, C. Yang, G. Ferrigno, and E. De Momi, “Improved human–robot collaborative control of redundant robot for teleoperated minimally invasive surgery,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 1447–1453, 2019.
  * [2] A. M. Zanchettin, N. M. Ceriani, P. Rocco, H. Ding, and B. Matthias, “Safety in human-robot collaborative manufacturing environments: Metrics and control,” _IEEE Transactions on Automation Science and Engineering_ , vol. 13, no. 2, pp. 882–893, 2016.
  * [3] F. Ferraguti, C. T. Landi, L. Sabattini, M. Bonfè, C. Fantuzzi, and C. Secchi, “A variable admittance control strategy for stable physical human–robot interaction,” _The International Journal of Robotics Research_ , vol. 38, no. 6, pp. 747–765, 2019. [Online]. Available: https://doi.org/10.1177/0278364919840415
  * [4] F. Dimeas and N. Aspragathos, “Online stability in human-robot cooperation with admittance control,” _IEEE Transactions on Haptics_ , vol. 9, no. 2, pp. 267–278, 2016.
  * [5] C. T. Landi, F. Ferraguti, C. Secchi, and C. Fantuzzi, “A passivity-based strategy for manual corrections in human-robot coaching,” _Electronics_ , vol. 8, no. 3, pp. 1–16, 2019.
  * [6] M. Karlsson, A. Robertsson, and R. Johansson, “Autonomous interpretation of demonstrations for modification of dynamical movement primitives,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , 2017, pp. 316–321.
  * [7] L. Villani and J. De Schutter, “Force control,” in _Springer Handbook of Robotics_ , B. Siciliano and O. Khatib, Eds. Cham: Springer International Publishing, 2016.
  * [8] Y. Aydin, O. Tokatli, V. Patoglu, and C. Basdogan, “Stable physical human-robot interaction using fractional order admittance control,” _IEEE Transactions on Haptics_ , vol. 11, no. 3, pp. 464–475, 2018.
  * [9] C. Secchi, S. Stramigioli, and C. Fantuzzi, _Control of interactive robotic interfaces: A port-Hamiltonian approach_. Springer Science & Business Media, 2007, vol. 29.
  * [10] F. Ferraguti, N. Preda, A. Manurung, M. Bonfè, O. Lambercy, R. Gassert, R. Muradore, P. Fiorini, and C. Secchi, “An energy tank-based interactive control architecture for autonomous and teleoperated robotic surgery,” _IEEE Transactions on Robotics_ , vol. 31, no. 5, pp. 1073–1088, 2015.
  * [11] S. Grafakos, F. Dimeas, and N. Aspragathos, “Variable admittance control in phri using emg-based arm muscles co-activation,” in _2016 IEEE International Conference on Systems, Man, and Cybernetics (SMC)_ , 2016, pp. 001 900–001 905.
  * [12] C. Secchi and F. Ferraguti, “Energy optimization for a robust and flexible interaction control,” in _ICRA 2019_ , 05 2019, pp. 1919–1925.
  * [13] M. Franken, S. Stramigioli, S. Misra, C. Secchi, and A. Macchelli, “Bilateral telemanipulation with time delays: A two-layer approach combining passivity and transparency,” _IEEE Transactions on Robotics_ , vol. 27, no. 4, pp. 741–756, 2011.
  * [14] 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, 2017.
  * [15] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in _2019 18th European Control Conference (ECC)_ , 2019, pp. 3420–3431.
  * [16] G. Notomista, S. F. Ruf, and M. Egerstedt, “Persistification of robotic tasks using control barrier functions,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 2, pp. 758–763, 2018.
  * [17] D. Zhang, G. Yang, and R. P. Khurshid, “Haptic teleoperation of uavs through control barrier functions,” _IEEE Transactions on Haptics_ , vol. 13, no. 1, pp. 109–115, 2020.
  * [18] J. Fu, Y. Lv, G. Wen, X. Yu, and T. Huang, “Velocity and input constrained coordination of second-order multi-agent systems with relative output information,” _IEEE Transactions on Network Science and Engineering_ , vol. 7, no. 3, pp. 1925–1938, 2020.
  * [19] G. Notomista and M. Egerstedt, “Persistification of robotic tasks,” _IEEE Transactions on Control Systems Technology_ , pp. 1–12, 2020.
  * [20] F. Ferraguti, M. Bertuletti, C. T. Landi, M. Bonfè, C. Fantuzzi, and C. Secchi, “A control barrier function approach for maximizing performance while fulfilling to iso/ts 15066 regulations,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 5921–5928, 2020.
  * [21] F. Ferraguti, C. Talignani Landi, S. Costi, M. Bonfè, S. Farsoni, C. Secchi, and C. Fantuzzi, “Safety barrier functions and multi-camera tracking for human-robot shared environment,” _Robotics and Autonomous Systems_ , vol. 124, p. 103388, 2020. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0921889019306426
  * [22] G. Notomista, S. Mayya, M. Selvaggio, M. Santos, and C. Secchi, “A set-theoretic approach to multi-task execution and prioritization,” 2020.
  * [23] G. Notomista, X. Cai, J. Yamauchi, and M. Egerstedt, “Passivity-based decentralized control of multi-robot systems with delays using control barrier functions,” in _2019 International Symposium on Multi-Robot and Multi-Agent Systems (MRS)_ , 2019, pp. 231–237.
  * [24] P. R. Giordano, A. Franchi, C. Secchi, and H. H. Bülthoff, “A passivity-based decentralized strategy for generalized connectivity maintenance,” _The International Journal of Robotics Research_ , vol. 32, no. 3, pp. 299–323, 2013\.
  * [25] G. Riggio, C. Fantuzzi, and C. Secchi, “On the use of energy tanks for multi-robot interconnection,” in _IROS 2018_ , 10 2018, pp. 3738–3743.
  * [26] G. Notomista and M. Egerstedt, “Constraint-driven coordinated control of multi-robot systems,” _CoRR_ , vol. abs/1811.02465, 2018. [Online]. Available: http://arxiv.org/abs/1811.02465
  * [27] G. Notomista, S. Mayya, S. Hutchinson, and M. Egerstedt, “An optimal task allocation strategy for heterogeneous multi-robot systems,” _CoRR_ , vol. abs/1903.08641, 2019.