# Spectral properties of soft X-ray transient MAXI J0637$-$430 using AstroSat

Neal Titus Thomas *  — Department of Physics and Electronics, CHRIST (Deemed
to be University)
Bangalore Central Campus, Bengaluru - 560029, India  Shivappa B. Gudennavar *
— Department of Physics and Electronics, CHRIST (Deemed to be University)
Bangalore Central Campus, Bengaluru - 560029, India  Ranjeev Misra  Inter-
University Centre of Astronomy and Astrophysics
Post Bag 4, Ganeshkhind, Pune - 411007, India  Bubbly S. G. *  — Department of
Physics and Electronics, CHRIST (Deemed to be University)
Bangalore Central Campus, Bengaluru - 560029, India

(Received : 22-10-2021; Accepted : 2-12-2021)

###### Abstract

Soft X-ray transients are systems that are detected when they go into an
outburst, wherein their X-ray luminosity increases several orders of
magnitude. These outbursts are markers of the poorly understood change in the
spectral state of these systems from low/hard state to high/soft state. We
report the spectral properties of one such soft X-ray transient: MAXI
J0637$-$430, with data from the SXT and LAXPC instruments on-board AstroSat
mission. The source was observed for a total of $\sim$ 60 ks over two
observations on 8th and 21st November, 2019 soon after its discovery. Flux
resolved spectral analysis of the source indicates the presence of a multi-
colour blackbody component arising from the accretion disk and a thermal
Comptonization component. The stable low temperature ($\sim$ 0.55 $keV$) of
the blackbody component, points to a cool accretion disk with an inner disk
radius of the order of a few hundred $km$. In addition, we report the presence
of a relativistically broadened Gaussian line at 6.4 $keV$. The disk dominated
flux and photon power law index of $\gtrapprox 2$ and a constant inner disk
radius indicate the source to be in the soft state. From the study we conclude
that MAXI J0637$-$430 is a strong black hole X-ray binary candidate.

Accretion — X-ray binaries — Transients — Black Hole physics

††journal: The Astrophysical Journal

: 3-2-2022

##  1 Introduction

Soft X-ray transients are a subclass of low mass X-ray binaries (LMXBs) that
appear as extremely faint sources ($L=10^{30}-10^{33}$
$erg\leavevmode\nobreak\ s^{-1}$) during most of their lifetime. They are
characterized by a non-steady transfer of mass onto the compact object and
they occasionally undergo sporadic outbursts, which occur at intervals of 1 -
60 years (Chen et al., 1997; Tetarenko et al., 2016). This causes their X-ray
luminosity to increase by a factor of upto $10^{7}$ (Paradijs & McClintock,
1995) which then decays back to quiescence with an e - folding timescale of
$\sim 30$ days (Chen et al., 1997). During the short outbursts, they emit
enough X-rays, which makes them the brightest X-ray sources
($L=10^{37}-10^{38}$ $erg\leavevmode\nobreak\ s^{-1}$) in the sky. The
occurrence of these outbursts is attributed to instabilities in the accretion
disk that are both viscous and thermal in nature (Meyer & Meyer-Hofmeister,
1981; Cannizzo et al., 1995; King & Ritter, 1998; Lasota, 2001). Soft X-ray
transients, especially the ones harbouring a black hole, usually go
undetected; and are discovered only when they undergo an outburst.
MAXI J0637$-$430 is one such source, which was first detected by the MAXI/GSC
nova alert system during seven scan transits from $2^{nd}$ $-$ $3^{rd}$
November, 2019 in the $2-4$ $keV$ and $4-10$ $keV$ bands. The scans revealed
the source to be located at RA (J2000) = 06 h 38 m 54 s, Dec (J2000) = $-$42 h
45 m 57 s in the soft band ($2-4$ $keV$) and at RA (J2000) = 06 h 37 m 43 s,
Dec (J2000) = $-$43 h 03 m 15 s in the hard band ($4-10$ $keV$) with 90%
confidence level (Negoro et al., 2019). Since its first detection, MAXI
J0637$-$430 underwent a considerable increase in its flux from $59$ $\pm$ $6$
mCrab to $\sim 200$ mCrab in the $2-4$ $keV$ band and from $32$ $\pm$ $6$
mCrab to $\sim 50$ mCrab in the $4-10$ $keV$ band (Negoro et al., 2019). A
Target of Opportunity (ToO) observation performed on $3^{rd}$ November, 2019
by the Swift mission detected MAXI J0637$-$430 at RA (J2000) = 06 h 36 m 23.59
s, Dec (J2000) = $-$42 h 52 m 04.1 s, which was consistent with MAXI’s hard
band localization. Its X-ray spectrum modelled using absorbed disk blackbody +
power law with inner disk temperature ($kT_{in}$) of 0.9 $\pm$ $0.1$ $keV$ and
power-law index ($\Gamma$) of $2.3$ $\pm$ $0.8$ indicated that the source
underwent an outburst and transitioned from hard to soft spectral state. An
optical counterpart with a brightness of u = $14.87$ $\pm$ $0.02$ (Vega) was
detected by the UVOT instrument on-board Swift at RA (J2000) = 06 h 36 m 23.23
s Dec (J2000) = $-$42 h 52 m 04.25 s. Since there are no known stars at this
position, it was deemed that this optical source underwent a significant
brightening, as is common for the optical counterpart of black hole low mass
X-ray binaries (BH-LMXBs) during outbursts (Kennea et al., 2019). Follow up
observation by the NuStar mission found the source to be at a flux of $\sim$
95 mCrab. Preliminary analysis of spectrum of the source in the energy range
$3-79$ $keV$ with a thermal disk blackbody component, power law and a
reflection component yielded a $kT_{in}$ of 0.628 $\pm$ 0.004 $keV$ and a
$\Gamma$ of 2.40 $\pm$ 0.04 (90% confidence errors) (Tomsick et al., 2019).
MAXI J0637$-$430 was also observed in radio band with the ATCA with flux
densities of $66\pm 15$ $\mu$Jy at 5.5 GHz and $60\pm 10$ $\mu$Jy at 9 GHz
(Russell et al., 2019). However, as the nature of the radio jet emission could
not be deciphered, the source could not be properly classified using the
radio/X-ray correlation in X-ray binaries. The source was also observed in
infrared band with the simultaneous imaging camera SIRIUS attached to 1.4 m
telescope InfraRed Survey Facility (IRSF), where the estimated magnitudes in
the J, H and K bands were $17.40\pm 0.01$, $17.69\pm 0.02$ and $17.96\pm
0.05$, respectively (Murata et al., 2019). Spectral analysis of the X-ray data
from Swift observations with an absorbed disk blackbody model showed the
source to have a $kT_{in}$ of 0.675 $\pm$ 0.003 $keV$ (Knigge et al., 2019),
which is consistent with the value obtained from NuStar observations. Since
its discovery, MAXI J0637$-$430 was observed by the NICER mission continuously
with a cadence of 1-2 days, which observed it undergo a spectral state
transition. After $\sim$ 23 days since its discovery, it was reported that
source transitioned into the hard state (Remillard et al., 2020). Subsequent
observations by Swift showed that MAXI J0637$-$430 could possibly be
approaching its quiescence level (Tomsick & Lazar, 2020). A multi-wavelength
study of the source was carried out by Tetarenko et al. (2021) using data from
Swift-XRT and UVOT, Gemini/GMOS, ATCA and AAVSO. This study made use of an
irradiated accretion disk model - $\mathtt{(diskir)}$ (Gierliński et al.,
2009) to derive the time-series evolution of its spectral parameters over the
entire outburst cycle. Analysis of NICER data by (Jana et al., 2021) revealed
the source to comprise of an ultra-soft thermal component ($kT_{in}\lesssim$
0.6 $keV$) and a power law tail. The study also showed that its spectra do not
need a thermal component corresponding to the emission from neutron star
surface, thus suggesting that the compact object in the MAXI J0637$-$430 is
most likely a black hole. Mass of the black hole inferred from this study is 5
$-$ 12 $M_{\odot}$ for a source distance of $d<10\leavevmode\nobreak\ kpc$ and
the distance to the source is found to have a lower limit of 6.5 $kpc$. Baby
et al. (2021) found that the 0.5 $-$ 25 $keV$ spectra of the source could be
modelled with a multi-colour disc emission $\mathtt{(diskbb)}$ convolved with
a thermal Comptonisation component $\mathtt{(thcomp)}$. Spectral fitting with
the $\mathtt{kerbb}$ model in conjunction with the soft-hard transition
luminosity, favour a black hole with mass between 3 $-$ 19 $M_{\odot}$ and
retrograde spin at a distance $<$ 15 $kpc$. Broadband spectral study on NuSTAR
data of the source showed that a two-component model, comprising of a
combination of multi-color disk blackbody and thermal Comptonization component
is adequate to fit the spectra only upto 10 $keV$. When higher energies are
considered, scenarios involving a plunging region and reprocessing of
returning disk radiation are equally possible (Lazar et al., 2021).
Encouraged by the observation campaigns carried out by various satellite and
ground based telescopes, ToO observations of MAXI J06347$-$430 were performed
using the Soft X-ray Telescope (SXT) and Large Area X-ray Proportional Counter
(LAXPC) instruments on-board AstroSat in the 0.3 $-$ 80 $keV$ energy range on
8th and 21st November, 2019. Here, we report the results of spectral and
temporal studies carried out on MAXI J0637$-$430 data from the SXT and LAXPC
instruments. The details of AstroSat observations and the data reduction
procedures are described in Section 2. In Section 3, lightcurve and hardness
intensity diagram (HID) are presented. In Section 4, we present the results of
spectral analysis. The findings and summary of the results are discussed in
Section 5.

##  2 Observations and Data reduction

ToO observations of MAXI J06347$-$430 (Thomas et al., 2019) in the 0.3 $-$ 80
$keV$ energy range were carried out using SXT and LAXPC on-board AstroSat for
a total of $\sim$ 60 ks on 8th (hereafter, Observation 1), 15th and 21st
November, 2019 (hereafter, Observation 2). We did not include the 15th
November data in our study as it contains 9-pointing safety observations for
the UVIT instrument on-board AstroSat, each with different pointing and
offset. Due to this, the spectra from the individual pointing could not be
combined as the effective area of the instrument changes with the offset.
Moreover, since the LAXPC pointings were also different, flux measurement
using SXT+LAXPC data could not be made. AstroSat observations used for our
study, marked on the 2 $-$ 20 $keV$ MAXI lightcurve in Figure 1 shows that
Observation 1 was carried out shortly after the outburst peak, whereas
Observation 2 was performed midway during the outburst decay. The Photon
Counting mode (PC) was employed for observation with SXT, whereas for LAXPC,
the observation was carried out in the Event Analysis Mode (EA). A log of
observations used for this study is given in Table 1.
SXT is a focusing telescope equipped with a Charged Coupled Device (CCD)
camera that performs X-ray imaging in the 0.3 $-$ 8.0 $keV$ energy range with
a spectral resolution of $\sim$ 150 $eV$ at 6 $keV$ (Singh et al., 2016). SXT
data of MAXI J06347$-$430 was processed using the standard SXT pipeline -
AS1SXTLevel2-1.4b111https://www.tifr.res.in/ãstrosat_sxt/sxtpipeline.html.
This yielded Level 2 event files for individual orbits of the observation,
which were then merged into one master event file using the SXT Event Merger
Tool222https://www.tifr.res.in/ãstrosat_sxt/dataanalysis.html. The merged
event file was then used to extract source images with the help of XSELECT
V2.4k. The source was selected between the region of 8′ and 5′ (inner radius)
for Observations 1 and 2, respectively, and 15′ (outer radius) to reduce the
effect of pile-up of the CCD. The response matrix333sxt_pc_mat_g0to12.rmf and
background444SkyBkg_comb_EL3p5_Cl_Rd16p0_v01.pha files provided by the SXT
Payload Operations Centre (POC) were used for the analysis. Off-axis Auxiliary
Response File (ARF) was created with the sxt_ARFModule2. SXT data in the range
0.5 $-$ 5.0 $keV$ for Region 1 and 3; and 0.5 $-$ 4.8 $keV$ for Region 2 were
used (Figure 4) as the data quality above and below these energy ranges was
poor.
LAXPC is a cluster of three co-aligned proportional counters (LAXPC-10,
LAXPC-20, LAXPC-30) that operates in the 3 $-$ 80 $keV$ energy range with an
absolute temporal resolution of 10 $\mu s$ (Yadav et al., 2016; Agrawal, 2017;
Antia et al., 2017). Data from LAXPC was processed with the LAXPCSOFT (Format
A)555http://astrosat-ssc.iucaa.in/?q=laxpcData to obtain event files, Good
Time Interval (GTI) files, lightcurves, source and background energy spectra,
Response Matrix Files (RMF) and power density spectra. LAXPC-20 data alone was
used for our study as it was reported by the POC that LAXPC-10 underwent an
abnormal change in its gain on $28^{th}$ March, 2018 and LAXPC-30 was not
operational during this time. Moreover, as the energy spectrum above 20 $keV$
was background dominated, spectral studies using LAXPC were restricted to 4
$-$ 20 $keV$ energy range.

##  3 Lightcurve and Hardness Intensity Diagram (HID)

Net lightcurves of the source were obtained in the 0.7 $-$ 7.0 $keV$ range
from the SXT instrument; and 4.0 $-$ 5.0 $keV$ and 5.0 $-$ 30.0 $keV$ ranges
from the LAXPC instrument. These lightcurves were binned to $\sim$ 50 s. It is
seen that during the beginning of the outburst decay i.e. in Observation 1,
the source intensity remains fairly constant at $\sim$ 30 counts/s and $\sim$
28 counts/s in the 4.0 $-$ 5.0 $keV$ (Panel 2 in Figure 2) and 5.0 $-$ 30
$keV$ (Panel 3 in Figure 2) ranges respectively. This then changes as the
count rates in both energy ranges increase by a small factor towards the end
of the observation. This jump in the count rate is reflected in the hardness-
time diagram too (Panel 4 in Figure 2). In comparison, the LAXPC net flux
along with the hardness ratio is seen to decrease monotonically through the
latter part of the outburst decay i.e. in Observation 2 (Panels 2, 3 and 4 in
Figure 3). Using LAXPC-20 data of both the observations, a combined, 50 s
binned HID was generated with hardness defined as the ratio of counts in the
5.0 $-$ 30.0 $keV$ range to 4.0 $-$ 5.0 $keV$ range and intensity defined as
the sum of counts in the 4.0 $-$ 30 $keV$ range. From the pattern traced by
the HID, it is not possible to determine if the source showed characteristics
of the q-diagram exhibited by BH-LMXBs (Remillard & McClintock, 2006) or the Z
or Atoll pattern exhibited by NS-LMXBs (Hasinger & van der Klis, 1989).
However, it is seen that the HID showed variability in hardness from $\sim$
0.9 to $\sim$ 1.9 (Figure 4). Further, in order to investigate the hardness-
intensity relation in energy range $<$ 4 $keV$, data from the SXT instrument,
corresponding to the three regions in the LAXPC-HID was used to generate an
HID (Figure 5). The hardness of this SXT-HID was defined as the ratio of
counts in the 1.0 $-$ 7.0 $keV$ range to 0.3 $-$ 1.0 $keV$ range and the
intensity was defined as the sum of counts in the 0.3 $-$ 7.0 $keV$ range, for
time corresponding to all the three regions from the LAXPC-HID. Regions 1 and
3 in this SXT-HID show variations in average flux and hardness, whereas Region
2 has too less data points to make a definitive conclusion.

Table 1: Observation log Obs. ID | Date | MJD | Exposure (ks)
---|---|---|---
|  |  | SXT | LAXPC
9000003290 | 08-11-2019 | 58795 | 8.7 | 6.5
9000003328 | 21-11-2019 | 58808 | 19.1 | 24.8
Figure 1: 1 day binned MAXI lightcurve of MAXI J0637$-$430 in 2 $-$ 20 $keV$
band Figure 2: SXT lightcurve of the source in 0.7 $-$ 7.0 $keV$ energy range
(Panel 1), LAXPC-20 net lightcurve in the energy ranges 4.0 $-$ 5.0 $keV$
(Panel 2) and 5.0 $-$ 30 $keV$ (Panel 3) and their hardness ratio (Panel 4)
from Observation 1 Figure 3: SXT lightcurve of the source in 0.7 $-$ 7.0 $keV$
energy range (Panel 1), LAXPC-20 net lightcurve in the energy ranges 4.0 $-$
5.0 $keV$ (Panel 2) and 5.0 $-$ 30 $keV$ (Panel 3) and their hardness ratio
(Panel 4) from Observation 2 Figure 4: HID of the source from LAXPC-20 data,
which is divided into 3 regions for flux resolved spectral studies. Figure 5:
HID of the source from SXT data showing Regions 1, 2 and 3 as resolved in the
LAXPC-HID.

##  4 Spectral Analysis

The simultaneous broadband X-ray spectral coverage of AstroSat with the SXT
and LAXPC instruments was used to perform flux resolved spectroscopy. This was
done by dividing the LAXPC-HID into 3 regions: Region 1, 2 and 3. The
divisions were made such that each region denotes an isolated cluster of
points in the LAXPC-HID (Figure 4). Simultaneous GTIs for both the SXT and
LAXPC instruments were generated, using which simultaneous spectra for both
instruments were generated for all three regions of the LAXPC-HID. The spectra
were then fit with the spectral modelling tool $\mathtt{XSPEC}$
$\mathtt{version:12.10.1o}$ (Arnaud, 1996) in the energy ranges 0.5 $-$ 5.0
$keV$ from the SXT instrument for Regions 1 and 3; and 0.5 $-$ 4.8 $keV$ for
Region 2; whereas form the LAXPC instrument, the energy range 4.0 $-$ 20 $keV$
was chosen for all three regions. These energy ranges were chosen as the
spectra below 0.5 and above 4.8 $keV$ (from the SXT instrument) showed very
high residuals and that above 20.0 $keV$ (from the LAXPC was dominated by the
background. A multi-colour blackbody model - $\mathtt{diskbb}$ (Mitsuda et
al., 1984), was used along with a convolution Comptonization model -
$\mathtt{simpl}$ (Steiner et al., 2009), in order to explain the emission from
the accretion disk and thermally Comptonized corona, respectively. The energy
range for the spectral fits were extended using $\mathtt{energies}$
$\mathtt{0.01}$ $\mathtt{100}$ $\mathtt{500}$ $\mathtt{log}$ to supply an
energy-binning array. To account for absorption in the interstellar medium, we
used the Tuebingen-Boulder Inter-Stellar Medium absorption model -
$\mathtt{tbabs}$, with the solar abundance table given by Wilms et al. (2000).
Further, a multiplicative constant factor was included to address
uncertainties caused due to cross calibration of the SXT and LAXPC
instruments. As prescribed by the POC, a systematic error of 3% was added to
all spectral fit (Bhattacharya, 2017). In addition to this, gain fit was
performed for the SXT data to account for the non-linear change in the
detector gain. Slope of the gain was frozen to 1 leaving the offset to vary.
There were positive residuals around 6.4 $keV$ indicating possible presence of
a disk reflection feature. A Gaussian component with its line energy frozen at
6.4 $keV$ was later added to account for this. The addition of the Gaussian
component yielded a small change in $\Delta\chi^{2}$ from 1.03 and 0.92 to 1.0
and 0.76 in Regions 1 and 2 respectively; whereas for Region 3 it remained
constant. The model combination -
$\mathtt{constant*tbabs(simpl*diskbb+gaussian)}$ yielded good fits for all the
three regions, with $\Delta\chi^{2}$ of $\sim$ 0.9 (Figure 6). The best fit
spectral parameters of this fit are given in Table 2 and the corresponding
spectra are give in Figure 6. The Norm of $\mathtt{diskbb}$ remains fairly
constant in all the three regions. Hence, in order to understand the observed
variation in the HID, spectral fit was repeated with the Norm of
$\mathtt{diskbb}$ fixed at 2075. This value was obtained by fitting a constant
through the Norms of all the three regions. Further, the unabsorbed total and
disk flux were calculated in the 0.5 $-$ 20 $keV$ range using the
$\mathtt{cflux}$ model. The best fit parameters of this spectral fit are given
in Table 3.
In addition to this, we also carried out temporal analysis with data from the
LAXPC instrument in the 4 $-$ 30 $keV$ range. However, it did not yield
substantial results.

Table 2: Best fit model spectral parameters Model | Parameter | Region 1 | Region 2 | Region 3
---|---|---|---|---
$\mathtt{tbabs}$ | $N_{H}(10^{20}cm^{-2})$ | 3.300 $\pm$ 0.006 | 2.120 $\begin{subarray}{c}+0.045\\\ -0.021\end{subarray}$ | 1.250 $\pm$ 0.005
$\mathtt{simpl}$ | $\Gamma$ | 2.00 $\begin{subarray}{c}+0.14\\\ -0.16\end{subarray}$ | 1.95 $\begin{subarray}{c}+0.17\\\ -0.20\end{subarray}$ | 2.46 $\pm\leavevmode\nobreak\ 0.03$
| $FracSctr$ | 0.016 $\begin{subarray}{c}+0.004\\\ -0.003\end{subarray}$ | 0.031 $\begin{subarray}{c}+0.010\\\ -0.009\end{subarray}$ | 0.100 $\begin{subarray}{c}+0.009\\\ -0.008\end{subarray}$
$\mathtt{diskbb}$ |  $kT_{in}$ (keV) | 0.610 $\pm 0.005$ | 0.650 $\pm$ 0.030 | 0.520 $\pm$ 0.007
| $Norm$ | 2103 $\begin{subarray}{c}+165\\\ -152\end{subarray}$ | 1455 $\begin{subarray}{c}+714\\\ -465\end{subarray}$ | 2083 $\begin{subarray}{c}+104\\\ -98\end{subarray}$
$\mathtt{Gaussian}$ |  $Line$ (keV) | 6.4 (f) | 6.4 (f) | 6.4 (f)
|  $Width$ (keV) | 1.06 $\begin{subarray}{c}+0.35\\\ -0.40\end{subarray}$ | 1.49 $\begin{subarray}{c}+0.42\\\ -0.47\end{subarray}$ | 0.20 $\begin{subarray}{c}+0.56\\\ -0.37\end{subarray}$
|  $Norm$ ($10^{-3}$) | 1.39 $\begin{subarray}{c}+0.60\\\ -0.50\end{subarray}$ | 3.24 $\begin{subarray}{c}+1.60\\\ -1.40\end{subarray}$ | 0.34$\pm\leavevmode\nobreak\ 0.30$
Reduced $\chi^{2}$ /dof |  | 1.03/363 | 0.76/109 | 1.33/460
Table 3: Best fit model spectral parameters with fixed $\mathtt{diskbb}$ Norm Model | Parameter | Region 1 | Region 2 | Region 3
---|---|---|---|---
$\mathtt{tbabs}$ | $N_{H}(10^{20}cm^{-2})$ | $3.31\pm 0.006$ | $3.87\begin{subarray}{c}+0.04\\\ -0.03\end{subarray}$ | $1.12\pm 0.005$
$\mathtt{simpl}$ | $\Gamma$ | $2.00\begin{subarray}{c}+0.14\\\ -0.16\end{subarray}$ | $2.00\begin{subarray}{c}+0.15\\\ -0.17\end{subarray}$ | $2.49\pm 0.03$
| $FracSctr$ | $0.016\begin{subarray}{c}+0.004\\\ -0.003\end{subarray}$ | $0.028\pm 0.007$ | $0.100\pm 0.008$
$\mathtt{diskbb}$ | $kT_{in}\leavevmode\nobreak\ \textit{(keV)}$ | $0.610\pm 0.003$ | $0.620\pm 0.005$ | $0.520\pm 0.005$
| $Norm$ | $2075(f)$ | $2075(f)$ | $2075(f)$
$\mathtt{Gaussian}$ | $Line\leavevmode\nobreak\ \textit{(keV)}$ | $6.4(f)$ | $6.4(f)$ | $6.4(f)$
| $Width\leavevmode\nobreak\ \textit{(keV)}$ | $1.06\begin{subarray}{c}+0.35\\\ -0.40\end{subarray}$ | $1.32\begin{subarray}{c}+0.34\\\ -0.42\end{subarray}$ | $0.37\begin{subarray}{c}+0.56\\\ -0.36\end{subarray}$
| $Norm\leavevmode\nobreak\ (10^{-3})$ | $1.39\begin{subarray}{c}+0.60\\\ -0.50\end{subarray}$ | $3.14\begin{subarray}{c}+1.50\\\ -1.30\end{subarray}$ | $0.40\begin{subarray}{c}+0.30\\\ -0.20\end{subarray}$
Unabsorbed disk flux | $(10^{-9}erg\leavevmode\nobreak\ cm^{-2}\leavevmode\nobreak\ s^{-1})$ | $5.31\begin{subarray}{c}+1.05\\\ -0.95\end{subarray}$ | $4.77\begin{subarray}{c}+1.20\\\ -0.84\end{subarray}$ | $2.58\begin{subarray}{c}+1.05\\\ -0.95\end{subarray}$
Unabsorbed total flux | $(10^{-9}erg\leavevmode\nobreak\ cm^{-2}\leavevmode\nobreak\ s^{-1})$ | $5.53\begin{subarray}{c}+1.05\\\ -0.95\end{subarray}$ | $5.15\begin{subarray}{c}+1.18\\\ -0.85\end{subarray}$ | $3.07\begin{subarray}{c}+1.04\\\ -0.96\end{subarray}$
Unabsorbed disk flux/total flux |  | $0.96$ | $0.92$ | $0.84$
$L/L_{Edd}$ |  | $0.025$ | $0.023$ | $0.141$
Reduced $\chi^{2}/dof$ |  | $1.02/364$ | $0.77/110$ | $1.34/461$

(a)

(b)

(c)

Figure 6: SXT+LAXPC unfolded spectra of Regions 1, 2 and 3 fit with the model
combination $\mathtt{constant*tbabs(simpl*diskbb+gaussian)}$. The residuals
($\chi=(data-model)/error$) are plotted in the bottom panels.

##  5 Results and Discussion

In this work, we analysed the SXT and LAXPC data from AstroSat observations
(8th and 21st November, 2019) of MAXI J0637$-$430 in the 0.5 $-$ 20 $keV$
energy range. The analysis revealed three distinct clusters in the LAXPC-HID
of the source (Figure 4). However, these clusters do not form a clear pattern
to give insights regarding the exact nature of the source. Monitoring the
source through its entire outburst cycle with NICER and MAXI has revealed its
HID to exhibit signatures of the various states of a BH-LMXB in outburst (Jana
et al., 2021; Baby et al., 2021). Flux resolved spectral analysis showed that
the spectra can be characterized by a multi-colour blackbody component arising
from the accretion disk along with a thermal Comptonization component. Our
choice of model - $\mathtt{diskbb}$, to characterise the soft component is in
agreement with erstwhile studies carried out on the source. However, different
models have been used to characterise its hard component - $\mathtt{powerlaw}$
(Tetarenko et al., 2021), $\mathtt{nthcomp}$ (Jana et al., 2021; Lazar et al.,
2021) and $\mathtt{thcomp}$ (Baby et al., 2021). The accretion disk
temperatures, 0.61, 0.65 and 0.52 $keV$ in Regions 1, 2 and 3 respectively,
point to a cool disk. This is in agreement with the studies carried out by
Tetarenko et al. (2021), Jana et al. (2021) and Baby et al. (2021), where the
disk temperature is seen to decay from $\sim$ 0.6 to $\sim$ 0.1 $keV$ during
the course of the outburst. The slightly lower value of disk temperature in
Region 3 correlates with the region being harder in the HID (Figure 4). The
$\mathtt{diskbb}$ Norms of 2103, 1455 and 2083 in Regions 1, 2 and 3, indicate
constant accretion disk radius throughout both the observations which is
consistent with that exhibited by many BH-LMXBs in the soft state (Done et
al., 2007). Similar results have been found by Baby et al. (2021) who used the
convolution model $\mathtt{thcomp}$ along with $\mathtt{diskbb}$ to
characterise the AstroSat spectra of the source. Assuming a source distance of
10 $kpc$, inclination angle of 70∘ and a colour hardening factor of 1.7
(Shimura & Takahara, 1995), we calculated the inner disk radius to be $\sim
98$, $\sim 81$ and $\sim 97$ $km$ in Regions 1, 2 and 3 respectively. For a
black hole of $20$ $M_{\odot}$, keeping the disk normalization constant at
2075, we estimated the inner disk radius to be $\sim 6$ $R_{g}$. This is the
distance at which the Innermost Stable Circular Orbit (ISCO) is located for a
non-rotating black hole. The increase in scatter fraction (obtained using
$\mathtt{simpl}$ model) from 0.016 $\begin{subarray}{c}+0.004\\\
-0.003\end{subarray}$ in Region 1 to 0.031 $\begin{subarray}{c}+0.01\\\
-0.009\end{subarray}$ in Region 2 reflects the increased LAXPC count rate in
the HID. In addition, the presence of a Gaussian line at 6.4 $keV$ points to a
reflection feature from the accretion disk. The width of this line ($\sim$ 1
$keV$ in Regions 1 and 2) shows that it is broadened due to relativistic
effects around the vicinity of the central compact object. It is to be noted
that the width and the Norm of the Gaussian component is significantly smaller
in Region 3. The unabsorbed total flux and disk flux of the source in the 0.5
$-$ 20 $keV$ range imply Eddington fraction of 0.025, 0.023 and 0.0141 for
Regions 1, 2 and 3 respectively (Table 3). The ratio of unabsorbed disk flux
to the total flux ($\sim$ 0.9) in all three Regions (Table 3), suggests that
the total flux is dominated by emission from the accretion disk. This
consistent disk dominated flux combined with the photon power law index of
$\gtrapprox$ 2, across all the three regions of the HID, shows the source to
be in soft state.

##  6 Conclusions

We carried out flux resolved spectral studies on two observations ($\sim$ 60
ks) of the soft X-ray transient source MAXI J0637 $-$ 430 in the 0.5 $-$ 20
$keV$ energy range using the SXT and LAXPC instruments on-board AstroSat.
Spectral analysis of shows the source to have a cool accretion disk having
temperature $\sim$ 0.55 $keV$ with a reflection feature at 6.4 $keV$. The
value of photon index of $\gtrapprox 2$ and the ratio of unabsorbed disk flux
to the total disk of $\sim 0.9$ points that MAXI J0637 $-$ 430 was in the soft
spectral state when observed by AstroSat. Also, it is seen that the value of
the disk normalization is consistent with being constant (Table 2) and points
to an accretion disk with an inner disk radius of 11.1 $R_{g}$. This is
observed in several BH-LMXBs in the soft state. We conclude from our study
that MAXI J0637$-$430 is a strong black hole X-ray binary candidate. Further
observations and in-depth studies of the source during its future outbursts
are essential to confirm its nature and unravel other physical parameters.

## Acknowledgements

We thank the SXT POC at TIFR, Mumbai and for the LAXPC POC team for their
support, timely release of data and providing the necessary software tools.
This work has made use of software provided by HEASARC. The authors
acknowledge the financial support of ISRO under AstroSat Archival Data
Utilization Program (No. DS-2B-13013(2)/9/2019-Sec.2 dated April 29, 2019).
This publication uses data from the AstroSat mission of the Indian Space
Research Organisation (ISRO), archived at the Indian Space Science Data Centre
(ISSDC). One of the authors (SBG) thanks the Inter-University Centre for
Astronomy and Astrophysics (IUCAA), Pune for the Visiting Associateship.

## References

  * Agrawal (2017) Agrawal, P. 2017, Journal of Astrophysics and Astronomy, 38, doi: 10.1007/s12036-017-9449-6
  * Antia et al. (2017) Antia, H. M., Yadav, J. S., Agrawal, P. C., et al. 2017, The Astrophysical Journal Supplement Series, 231, 10, doi: 10.3847/1538-4365/aa7a0e
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Baby et al. (2021) Baby, B. E., Bhuvana, G. R., Radhika, D., et al. 2021, MNRAS, 508, 2447, doi: 10.1093/mnras/stab2719
  * Bhattacharya (2017) Bhattacharya, D. 2017, Journal of Astrophysics and Astronomy, 38, 51, doi: 10.1007/s12036-017-9461-x
  * Cannizzo et al. (1995) Cannizzo, J. K., Chen, W., & Livio, M. 1995, ApJ, 454, 880, doi: 10.1086/176541
  * Chen et al. (1997) Chen, W., Shrader, C. R., & Livio, M. 1997, ApJ, 491, 312, doi: 10.1086/304921
  * Done et al. (2007) Done, C., Gierliński, M., & Kubota, A. 2007, A&A Rev., 15, 1, doi: 10.1007/s00159-007-0006-1
  * Gierliński et al. (2009) Gierliński, M., Done, C., & Page, K. 2009, MNRAS, 392, 1106, doi: 10.1111/j.1365-2966.2008.14166.x
  * Hasinger & van der Klis (1989) Hasinger, G., & van der Klis, M. 1989, A&A, 225, 79
  * Jana et al. (2021) Jana, A., Jaisawal, G. K., Naik, S., et al. 2021, Monthly Notices of the Royal Astronomical Society, 504, 4793, doi: 10.1093/mnras/stab1231
  * Kennea et al. (2019) Kennea, J. A., Bahramian, A., Evans, P. A., et al. 2019, The Astronomer’s Telegram, 13257, 1
  * King & Ritter (1998) King, A. R., & Ritter, H. 1998, MNRAS, 293, L42, doi: 10.1046/j.1365-8711.1998.01295.x
  * Knigge et al. (2019) Knigge, C., Rao, A., Gandhi, P., & Altamirano, D. 2019, The Astronomer’s Telegram, 13295, 1
  * Lasota (2001) Lasota, J.-P. 2001, New Astron. Rev., 45, 449, doi: 10.1016/S1387-6473(01)00112-9
  * Lazar et al. (2021) Lazar, H., Tomsick, J. A., Pike, S. N., et al. 2021, The Astrophysical Journal, 921, 155, doi: 10.3847/1538-4357/ac1bab
  * Meyer & Meyer-Hofmeister (1981) Meyer, F., & Meyer-Hofmeister, E. 1981, A&A, 104, L10
  * Mitsuda et al. (1984) Mitsuda, K., Inoue, H., Koyama, K., et al. 1984, PASJ, 36, 741
  * Murata et al. (2019) Murata, K. L., Adachi, R., Kawai, N., et al. 2019, The Astronomer’s Telegram, 13292, 1
  * Negoro et al. (2019) Negoro, H., Miike, K., Nakajima, M., et al. 2019, The Astronomer’s Telegram, 13256, 1
  * Paradijs & McClintock (1995) Paradijs, J., & McClintock, J. 1995, Optical and Ultraviolet Observations of X-Ray Binaries, Cambridge University Press
  * Remillard et al. (2020) Remillard, R., Pasham, D., Gendreau, K., et al. 2020, The Astronomer’s Telegram, 13427, 1
  * Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49, doi: 10.1146/annurev.astro.44.051905.092532
  * Russell et al. (2019) Russell, T. D., Miller-Jones, J. C. A., Sivakoff, G. R., & Tetarenko, A. J. 2019, The Astronomer’s Telegram, 13275, 1
  * Shimura & Takahara (1995) Shimura, T., & Takahara, F. 1995, ApJ, 445, 780, doi: 10.1086/175740
  * Singh et al. (2016) Singh, K. P., Stewart, G. C., Chandra, S., et al. 2016, in Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, T. Takahashi, & M. Bautz, Vol. 9905, International Society for Optics and Photonics (SPIE), 389 – 398, doi: 10.1117/12.2235309
  * Steiner et al. (2009) Steiner, J. F., Narayan, R., McClintock, J. E., & Ebisawa, K. 2009, PASP, 121, 1279, doi: 10.1086/648535
  * Tetarenko et al. (2021) Tetarenko, B. E., Shaw, A. W., Manrow, E. R., et al. 2021, MNRAS, 501, 3406, doi: 10.1093/mnras/staa3861
  * Tetarenko et al. (2016) Tetarenko, B. E., Sivakoff, G. R., Heinke, C. O., & Gladstone, J. C. 2016, ApJS, 222, 15, doi: 10.3847/0067-0049/222/2/15
  * Thomas et al. (2019) Thomas, N. T., Gudennavar, S. B., Misra, R., & Bubbly, S. G. 2019, The Astronomer’s Telegram, 13296, 1
  * Tomsick & Lazar (2020) Tomsick, J. A., & Lazar, H. 2020, The Astronomer’s Telegram, 13800, 1
  * Tomsick et al. (2019) Tomsick, J. A., Garcia, J., Fabian, A., et al. 2019, The Astronomer’s Telegram, 13270, 1
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914, doi: 10.1086/317016
  * Yadav et al. (2016) Yadav, J. S., Agrawal, P. C., Antia, H. M., et al. 2016, in Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, T. Takahashi, & M. Bautz, Vol. 9905, International Society for Optics and Photonics (SPIE), 374 – 388, doi: 10.1117/12.2231857


# Learning Model Predictive Control with Error Dynamics Regression for
Autonomous Racing

Haoru Xue1, Edward L. Zhu2, and Francesco Borrelli2 1Haoru Xue is with the
Robotics Institute, Carnegie Mellon University, Pittsburgh, PA, 15213, USA
<EMAIL_ADDRESS>L. Zhu and Francesco Borrelli are with the
Department of Mechanical Engineering, University of California Berkeley,
Berkeley, CA 94720, USA {edward.zhu<EMAIL_ADDRESS>
###### Abstract

This work presents a novel Learning Model Predictive Control (LMPC) strategy
for autonomous racing at the handling limit that can iteratively explore and
learn unknown dynamics in high-speed operational domains. We start from
existing LMPC formulations and modify the system dynamics learning method. In
particular, our approach uses a nominal, global, nonlinear, physics-based
model with a local, linear, data-driven learning of the error dynamics. We
conduct experiments in simulation, 1/10th scale hardware, and deployed the
proposed LMPC on a full-scale autonomous race car used in the Indy Autonomous
Challenge (IAC) with closed loop experiments at the Putnam Park Road Course in
Indiana, USA. The results show that the proposed control policy exhibits
improved robustness to parameter tuning and data scarcity. Incremental and
safety-aware exploration toward the limit of handling and iterative learning
of the vehicle dynamics in high-speed domains is observed both in simulations
and experiments.

##  I Introduction

Recent progress in autonomous racing research has empowered full-size
autonomous race cars at or near the top speed of a professional human racing
series. Algorithms that were previously experimented with on small-scale
vehicles [1] are deployed and further researched on professional race car
platforms. Specifically, the teams in the Indy Autonomous Challenge (IAC) [2],
starting in 2021, have been demonstrating multi-agent competition above 240
$\mathrm{km}\text{/}\mathrm{h}$ [3].

In the field of control, which is tasked with trajectory tracking at the limit
of vehicle dynamics handling, various Model Predictive Control (MPC)
techniques remain the state-of-the-art on full-size race car systems. These
variations of MPC algorithms solve a finite-horizon optimal control problem
(FHOCP) by computing optimal actuation commands for a vehicle modeled with
kinematic or dynamic models, subject to certain track boundary and state
constraints [4, 5, 6, 7]. However, for high-speed racing, the performance of
classic, nonadaptive MPC methods is ultimately bounded by the fidelity of the
vehicle dynamics model. Inaccurate vehicle models will not only induce
suboptimal behavior such as weaving, hunting, and steady-state tracking error,
but also become dangerous at the handling limit, especially when the model
overpredicts this limit. This introduces a chicken-and-egg problem for
autonomous racing research at the handling limit: to safely achieve higher
speed levels requires more vehicle dynamics data; but to acquire high-speed
data requires the controller to achieve higher speed levels in the first
place. Therefore, the motivation of this work is to break the circular problem
mentioned earlier by performing incremental and safety-aware exploration
toward the limit of handling and iterative learning of the vehicle dynamics in
high-speed domains.

Figure 1: Overlay of the IAC car’s trajectory running LMPC with error dynamics
regression. Thw data is collected at the Putnam Park Road Course in Indiana,
USA. The zoomed-in image of the hairpin turn shows iterative improvement
towards optimal cornering.

In this work, we present a modification to the LMPC approach presented in [8]
where we improve significantly upon the model learning strategy. In
particular, we propose to learn the error dynamics between a given nominal
model and the state evolution of the true system collected during run-time.
This is in contrast with the prior approach which attempts to learns a model
of the true system directly. We show through simulation and hardware
experiments on a 1/10th scale vehicle that our approach exhibits improved
robustness to parameter tuning and reduced modeling errors in scenarios when
data is scarce for model learning. We finally demonstrate the effectiveness of
our approach on a full-size IAC autonomous race car, the results of which are
shown in Figure 1. An open source C++ implementation of our work has been made
available at https://github.com/MPC-Berkeley/Racing-LMPC-ROS2.

_Related Work:_ Learning unknown or changing system dynamics has been studied
within some classic control frameworks such as adaptive control [9, 10].
System identification is usually part of the controller design, which reasons
about some modeling error or noise in the system and adapts the system model
to account for them. In autonomous racing, previous works have proposed
dynamics learning strategies using Gaussian Processes [11, 12, 13] or
implicitly with Gaussian noise assumption [14]. The method is further extended
to model multi-vehicle interactions and solve for the optimal overtaking
maneuver [15, 16].

Another branch of works exploit the repetitive nature of autonomous racing in
a single-vehicle scenario, and improves controller performance without any
prior assumption about the distribution of the modeling error or noise. Data-
driven methods, such as iterative learning control, can optimize the tracking
of a repetitive trajectory given data from previous iterations by solving for
control variable corrections [17]. More recently, the learning model
predictive control (LMPC) framework has been proposed with a direct regression
of local affine dynamics [18, 19, 8]. This combination of LMPC and regressive
dynamics learning can iteratively learn the unknown system dynamics and reduce
lap-time without an explicit prior reference trajectory or a nominal dynamics
model. The approach has also been demonstrated on scaled vehicle hardware
platforms in [18].

##  II Problem Formulation

In this work, we model the vehicle using the planar dynamic bicycle model [20,
21] with state and input vectors:

$\displaystyle x$
$\displaystyle=[v_{x},v_{y},\omega_{z},e_{\psi},s,e_{y}]^{\top}\in\mathbb{R}^{6}$
$\displaystyle u$ $\displaystyle=[a,\delta]^{\top}\in\mathbb{R}^{2}$

where $v_{x}$, $v_{y}$, and $\omega_{z}$ are the vehicle’s longitudinal
velocity, lateral velocity, and yaw rate at the center of gravity. $s$,
$e_{y}$, and $e_{\psi}$ are the pose of the vehicle expressed in the Frenet
frame w.r.t. a parametric path $\tau:[0,L]\mapsto\mathbb{R}^{2}$, which is
assumed to be continuously twice differentiable and of total length $L$. In
particular, when given an arclength $s\in[0,L]$, $\tau$ returns the global
$x$-$y$ position of the path. In the case where $\tau$ forms a closed circuit,
i.e. $\tau(0)=\tau(L)$, $\tau^{\prime}(0)=\tau^{\prime}(L)$,
$\tau^{\prime\prime}(0)=\tau^{\prime\prime}(L)$, where $\tau^{\prime}$ and
$\tau^{\prime\prime}$ are the element-wise first and second derivatives of
$\tau$, we define the circuit as $\bar{\tau}(s)=\tau(\text{mod}(s,L))$ for
$s\geq 0$, where mod is the modulo operator. Given a path $\tau$ and a vehicle
with global position $p=(x,y)$ and heading $\psi$, the Frenet frame pose is
defined as the vehicle’s path progress $s(p)=\arg\min_{s}\|\tau(s)-p\|_{2}$,
its lateral deviation from the path $e_{y}(p)=\min_{s}\|\tau(s)-p\|_{2}$, and
its heading deviation from the path tangent
$e_{\psi}(p,\psi)=\psi-\arctan(\tau^{\prime}(s(p)))$. The inputs to the model
are the longitudinal acceleration $a$ and the front steering angle $\delta$.
We write the discretized nonlinear dynamics as:

$\displaystyle x_{k+1}=f(x_{k},u_{k})$ (1)

The vehicle’s state and input are subject to the constraints
$\mathcal{X}=\\{x\ |\ -W/2\leq e_{y}\leq W/2\\}$ and $\mathcal{U}=\\{u\ |\
u_{l}\leq u\leq u_{u}\\}$, which describe the boundary constraints for a track
of constant width $W$ and the limits on achievable acceleration and steering
angle.

The objective of this work is to design a control policy which solves the
infinite-horizon, minimium lap time, optimal control problem:

$\displaystyle T=\min_{u_{0},u_{1},\dots}\ $ $\displaystyle\
\sum_{k=0}^{\infty}\mathds{1}_{\mathcal{F}}(x_{k})$ (2a) subject to
$\displaystyle\ x_{0}=\bar{x},$ (2b) $\displaystyle\ x_{k+1}=f(x_{k},u_{k}),$
(2c) $\displaystyle\ x_{k}\in\mathcal{X},\ u_{k}\in\mathcal{U},$ (2d)

where $\bar{x}\in\mathcal{S}=\\{x\ |\ s=0\\}$ is a state at the start of the
track and $\mathcal{F}=\\{x\ |\ s\geq L\\}$ is the set of states beyond the
finish line of the circuit. If for some time step $k$, $x_{k}\in\mathcal{F}$
for the first time, then the vehicle has finished the lap and has transitioned
onto the next lap. $\mathds{1}_{\mathcal{F}}$ is the indicator function for
the set $\mathcal{F}$, which is defined as

$\displaystyle\mathds{1}_{\mathcal{F}}(x)=\begin{cases}0\quad\text{if}\
x\in\mathcal{F}\\\ 1\quad\text{o.w.}\end{cases}$

Due to the infinite horizon cost in (2a), it should be clear that the optimal
control problem posed in (2) is not amenable to online implementation.
Moreover, it is likely that there is mismatch between the vehicle model in
(2c) and the true vehicle dynamics, the effects of which would be especially
apparent at high speeds where the vehicle is operating at the limit of
handling. To address these challenges, previous work on LMPC [19, 22] proposed
a finite-horizon convex approximation of (2), which uses state and input
trajectories from prior laps to synthesize a convex terminal cost-to-go
function and target set. The same data is additionally used to identify affine
time-varying (ATV) models which approximate the true dynamics of the vehicle.
In this work, we carry over the same approach for terminal cost-to-go function
and target set synthesis, but propose a modification to the system
identification procedure which results in an LMPC control policy that achieves
similar performance to [19, 22] w.r.t. minimum lap time, while exhibiting a
significant improvement in robustness to tuning parameters.

##  III Local Learning Model Predictive Control

In this section, we briefly describe the procedure for synthesizing LMPC
policies over iterations (or laps) of the minimum time control task. For
further details, we refer the reader to [19, 22].

The main idea of LMPC is to use the data from iterations $1$ to $j-1$ to
synthesize a finite-horizon local convex approximation of (2) for iteration
$j$, which is then solved in a receding horizon manner. Let us first denote
the available dataset at iteration $j$ as
$\mathcal{D}^{j}=(\mathbf{x}^{j},\mathbf{u}^{j})\cup\mathcal{D}^{j-1}$, where
$\mathbf{x}^{j}=\\{x_{0}^{j},x_{1}^{j},\dots,x_{T^{j}}^{j}\\}$ and
$\mathbf{u}^{j}=\\{u_{0}^{j},u_{1}^{j},\dots,u_{T^{j}}^{j}\\}$ are the closed-
loop state and input sequence for iteration $j$ and $T^{j}$ denotes the time
step at which the lap was completed for iteration $j$, i.e. the first time
step where $\mathds{1}_{\mathcal{F}}(x_{k})=0$. Note that on the first
iteration, the dataset $\mathcal{D}^{0}$ may be initialized using human-driven
laps or closed-loop trajectories from a simple low-speed center-line tracking
controller. In addition, we assume that the trajectories stored in
$\mathcal{D}^{j}$ are feasible w.r.t. the constraints (2d) and reach the
finish line in a finite amount of time.

###  III-A Local Convex Target Set

As the closed-loop trajectories stored in $\mathcal{D}^{j}$ are from the
actual vehicle, it is straightforward to see that if we can control the
vehicle to any state in $\mathcal{D}^{j}$, then there exists a known feasible
control sequence which can be used to complete the lap from that state. LMPC
uses this fact to construct a safe set as a discrete collection of the states
in $\mathcal{D}^{j}$, which is then used as the terminal set for a FHOCP [22].
In this work, we construct convex terminal sets which are local about a given
state $x$. These terminal sets sacrifice the theoretical guarantee of
persistent feasibility for the FHOCP in favor of reduced computational
complexity for fast online control. In particular, for some given state $x$,
we first take the $K$-nearest neighbors from the each of the previous $P$ laps
by evaluating the weighted Euclidean distance over the relevant states:

$\displaystyle(x_{k}^{i}-x)^{\top}D(x_{k}^{i}-x),\ \forall$ $\displaystyle
i\in\\{j,j-1,\dots,j-(P-1)\\},$ $\displaystyle k\in\\{0,1,\dots,T^{i}\\}$ (3)

for some given $P$ and $D\succeq 0$. Let
$\mathbf{X}^{j}(x;\mathcal{D}^{j})\in\mathbb{R}^{n\times KP}$ denote the
matrix formed by the $KP$ states closest to $x$, we then define the target set
as the convex hull of these states:

$\displaystyle\mathcal{X}_{N}^{j}(x;\mathcal{D}^{j})=\\{\bar{x}\in\mathbb{R}^{n}\
|\ \exists\lambda\in\mathbb{R}^{KP}$ $\displaystyle,\ 0\leq\lambda\leq 1,\
\mathbf{1}^{\top}\lambda=1,$
$\displaystyle\mathbf{X}^{j}(x;\mathcal{D}^{j})\lambda=\bar{x}\\}$

###  III-B Local Terminal Cost-to-go Function

By assumption, the closed-loop state trajectories stored in $\mathcal{D}^{j}$
reach the finish line in a finite amount of time. We can therefore compute the
cost-to-go for any state $x_{k}^{i}$ in $\mathcal{D}^{j}$ as $T^{i}-k$, which
represents the time remaining to finish the $i$-th lap for
$i\in\\{0,\dots,j\\}$. Let $J_{N}^{j}(x;\mathcal{D}^{j})\in\mathbb{R}^{KP}$
denote the vector of cost-to-go values calculated for the states in
$\mathbf{X}^{j}(x;\mathcal{D}^{j})$. We can therefore define the minimum cost-
to-go function for any $\bar{x}\in\mathcal{X}_{N}^{j}(x;\mathcal{D}^{j})$ as:

$\displaystyle Q_{N}^{j}(\bar{x},x;\mathcal{D}^{j})$
$\displaystyle=\min_{\lambda\in\mathbb{R}^{KP}}\ $ $\displaystyle\
J_{N}^{j}(x;\mathcal{D}^{j})^{\top}\lambda$ subject to $\displaystyle\
0\leq\lambda\leq 1,\ \mathbf{1}^{\top}\lambda=1,\
\mathbf{X}^{j}(x;\mathcal{D}^{j})\lambda=\bar{x}$

This function can be thought of as a local convex approximation to the optimal
cost-to-go (from solving (2)) at iteration $j$, which when incorporated as the
terminal cost-to-go function, allows an FHOCP to reason about infinite-horizon
behavior of the vehicle, to the extent as it is captured in the dataset.

###  III-C Local LMPC

Using the terminal target set and cost-to-go function synthesized from the
dataset $\mathcal{D}^{j-1}$, we may now construct the FHOCP at time step $k$
of iteration $j$ for our local LMPC policy as follows:

$\displaystyle J_{k}^{j}(x_{k},u$
${}_{k-1},\bar{\mathbf{z}}_{k};\mathcal{D}^{j-1})=$
$\displaystyle\min_{\mathbf{x},\mathbf{u},\lambda}\ $ $\displaystyle\
\sum_{t=0}^{N-1}\mathds{1}_{\mathcal{F}}(x_{t})+c_{u}\|u_{t}\|_{2}^{2}+c_{\Delta
u}\|u_{t}-u_{t-1}\|_{2}^{2}$
$\displaystyle\qquad+J_{N}^{j-1}(\bar{x}_{k+N};\mathcal{D}^{j-1})^{\top}\lambda$
(4a) subject to $\displaystyle\ x_{0}=x_{k},\ u_{-1}=u_{k-1}$ (4b)
$\displaystyle\
x_{t+1}=A(\bar{z}_{k+t};\mathcal{D}^{j-1})x_{t}+B(\bar{z}_{k+t};\mathcal{D}^{j-1})u_{t}$
$\displaystyle\quad+C(\bar{z}_{k+t};\mathcal{D}^{j-1}),\
t\in\\{0,\dots,N-1\\},$ (4c) $\displaystyle\ x_{t}\in\mathcal{X},\
u_{t}\in\mathcal{U},\ t\in\\{0,\dots,N\\}$ (4d) $\displaystyle\
\mathbf{X}^{j-1}(\bar{x}_{k+N};\mathcal{D}^{j-1})\lambda=x_{N},$ (4e)
$\displaystyle\ 0\leq\lambda\leq 1,\ \mathbf{1}^{\top}\lambda=1$ (4f)

where $\bar{z}_{k}=(\bar{x}_{k},\bar{u}_{k})$ and
$\bar{\mathbf{z}}_{k}=\\{\bar{z}_{k},\dots,\bar{z}_{k+N}\\}$ is a state and
input sequence which is used to form the local convex approximation. In
practice, this is chosen as the solution to (4) from the previous time-step.
Construction of the $A$, $B$, and $C$ matrices in the ATV dynamics (4c) is the
main contribution of this work and will be discussed in the next section. Note
that (4) is a convex program which can be solved efficiently using existing
solvers.

Finally, let $\mathbf{u}^{\star}=\\{u_{0}^{\star},\dots,u_{N-1}^{\star}\\}$ be
the optimal control sequence which solves (4) at time step $k$. The input
$u_{k}=u_{0}^{\star}$ is applied to the vehicle and the FHOCP (4) is repeated
at time step $k+1$ for the measured state $x_{k+1}$ until the vehicle reaches
the finish line, at which point $\mathcal{D}^{j}$ will be constructed using
the closed-loop trajectory for lap $j$.

##  IV Learned Vehicle Dynamics Model

We now introduce the main contribution of our work, which is a learning
strategy on the dynamics states $v_{x},v_{y},\omega_{z}$ based on regressive
error. The learning process take as inputs the nominal dynamics model and a
set of neighboring states, and outputs a learned local affine approximation of
the true vehicle dynamics.

We assume that the true dynamics of the system can be written as follows

$\displaystyle x^{+}=f(x,u)+e(x,u),$ (5)

where $f$ are the nominal dynamics from (1) and $e$ denotes some unknown
modeling error which we would like to identify using data. For a given state
and input $x$ and $u$, let us denote the prediction of the nominal model as
$\hat{x}^{+}=f(x,u)$. We may now linearize the error dynamics about a
reference state and input $\bar{z}$ as follows:

$\displaystyle x^{+}-\hat{x}^{+}=A^{e}x+B^{e}u+C^{e},$ (6)

where $A^{e}$ and $B^{e}$ are the Jacobians of $e$ w.r.t. $x$ and $u$
evaluated at $\bar{z}$ and
$C^{e}=e(\bar{x},\bar{u})-A^{e}\bar{x}-B^{e}\bar{u}$. It can be clearly seen
from (6) that the linearization of the unknown error term is related to the
system state and input $x$ and $u$, the actual state evolution $x^{+}$, and
the prediction of the nominal model $\hat{x}^{+}$, which are all quantities
that can be easily obtained from the dataset $\mathcal{D}^{j-1}$ at iteration
$j$. As such, we query $\mathcal{D}^{j-1}$ to find the $M$ nearest neighbors
$\mathbf{z}=\\{z_{1},\dots,z_{M}\\}$ of $\bar{z}$ and their corresponding
state evolutions $\mathbf{x}^{+}=\\{x_{1}^{+},\dots,x_{M}^{+}\\}$ (using a
technique similar to (III-A)). We then compute the prediction of nominal model
at each of the state and input pairs in $\mathbf{z}$ to obtain
$\hat{\mathbf{x}}^{+}=\\{\hat{x}_{1}^{+},\dots,\hat{x}_{M}^{+}\\}$. As
mentioned before, we are interested in learning only the error dynamics
associated with the velocity components of the state $x$ as the kinematic
components of the state are well-understood. Therefore, using the datasets
$\mathbf{z}$, $\mathbf{x}^{+}$, and $\hat{\mathbf{x}}^{+}$, we define the
residuals for each velocity state as follows:

$\displaystyle{A^{e}}[11:13]x_{m}[1:3]+{B^{e}}[11]a_{m}+{C^{e}}[1]$
$\displaystyle=v_{x,m}^{+}-\hat{v}_{x,m}^{+},$
$\displaystyle{A^{e}}[21:23]x_{m}[1:3]+{B^{e}}[22]\delta_{m}+{C^{e}}[2]$
$\displaystyle=v_{y,m}^{+}-\hat{v}_{y,m}^{+},$
$\displaystyle{A^{e}}[31:33]x_{m}[1:3]+{B^{e}}[32]\delta_{m}+{C^{e}}[3]$
$\displaystyle=\omega_{z,m}^{+}-\hat{\omega}_{z,m}^{+},$ (7)

where $m\in\\{1,\dots,M\\}$ and we use the notation $A[11:13]$ to index the
first three elements of the first row of matrix $A$. Note that we have
injected additional structure into the individual residuals by making explicit
the dependence of each velocity state on a certain subset of the input
components. From these three expressions, we define the following regressor
vectors:

$\displaystyle\Gamma_{v_{x}}=\begin{bmatrix}{\scriptstyle{A^{e}}[11:13]}\\\
{\scriptstyle{B^{e}}[11]}\\\
{\scriptstyle{C^{e}}[1]}\end{bmatrix},\Gamma_{v_{y}}=\begin{bmatrix}{\scriptstyle{A^{e}}[21:23]}\\\
{\scriptstyle{B^{e}}[22]}\\\
{\scriptstyle{C^{e}}[2]}\end{bmatrix},\Gamma_{\omega_{z}}=\begin{bmatrix}{\scriptstyle{A^{e}}[31:33]}\\\
{\scriptstyle{B^{e}}[32]}\\\ {\scriptstyle{C^{e}}[3]}\end{bmatrix}.$

We solve the regression problem in a manner similar to [8], where we weigh the
importance of data point $m$ using the Epanechnikov kernel function [23] with
bandwidth $h$:

$\displaystyle K(u)$
$\displaystyle=\begin{cases}\frac{3}{4}(1-u^{2}/h^{2}),&|u|<1\\\
0,&\text{otherwise}\end{cases}.$

such that larger penalties are assigned to the residuals for data points that
are close to the linearization point $\bar{x}$. For
$l=\\{v_{x},v_{y},\omega_{z}\\}$, the weighted least squares problem is
defined as follows:

$\displaystyle\Gamma_{l}^{\star}=\operatorname*{argmin}_{\Gamma_{l}}\sum_{m=1}^{M}K\left(||\bar{x}-x_{m}||^{2}_{Q}\right)y_{m}^{l}(\Gamma_{l}),$

where $Q$ denotes the relative scaling between the variables, and
$y_{m}^{l}(\Gamma_{l})$ are the $l^{2}$-norms of the residuals in (7).

Finally, we can unpack the optimal regressor vectors $\Gamma_{l}^{\star}$ to
construct the matrices $A^{e}$, $B^{e}$, and $C^{e}$, and add them to the
nominal ATV dynamics to obtain the complete ATV model for (4c):

$\displaystyle A(\bar{z};\mathcal{D}^{j-1})$
$\displaystyle=A^{f}+\begin{bmatrix}\begin{matrix}\Gamma_{v_{x}}^{\star}[1:3]\\\
\Gamma_{v_{y}}^{\star}[1:3]\\\ \Gamma_{\omega_{z}}^{\star}[1:3]\\\
\mathbf{0}_{3\times 3}\end{matrix}&\mathbf{0}_{6\times 3}\end{bmatrix},$
$\displaystyle B(\bar{z};\mathcal{D}^{j-1})$
$\displaystyle=B^{f}+\begin{bmatrix}\begin{matrix}\Gamma_{v_{x}}^{\star}[4]&0\\\
0&\Gamma_{v_{y}}^{\star}[4]\\\ 0&\Gamma_{\omega_{z}}^{\star}[4]\end{matrix}\\\
\mathbf{0}_{3\times 2}\end{bmatrix},$ $\displaystyle
C(\bar{z};\mathcal{D}^{j-1})$
$\displaystyle=C^{f}+\begin{bmatrix}\Gamma_{v_{x}}^{\star}[5]&\Gamma_{v_{y}}^{\star}[5]&\Gamma_{\omega_{z}}^{\star}[5]&\mathbf{0}_{1\times
3}\end{bmatrix}^{\top},$

where $A^{f}$ and $B^{f}$ are the Jacobians of $f$ w.r.t. $x$ and $u$
evaluated at $\bar{z}$ and
$C^{f}=f(\bar{x},\bar{u})-A^{f}\bar{x}-B^{f}\bar{u}$. We note that the key
difference of this procedure compared to [8] is that the regressive learning
is performed on the error between the prediction of the nominal model and the
actual state evolution data instead of on the entire model in (5) all at once.
This strategy not only provides a safety fallback when the regressive learning
fails at run-time, but also demonstrates, in comparison to [8], improved
robustness to parameter tuning and reduced modeling errors in scenarios when
data is scarce.

##  V Experiments & Results

We now present a set of experiments in simulation and hardware to demonstrate
the effectiveness of combining LMPC and error dynamics regression in
autonomous racing. The key metrics throughout these experiments are

  * •

20th iteration lap time (ILT-20): the average lap time after running the LMPC
for 20 iterations. The shorter the lap time, the better the performance of the
control policy.

  * •

Iteration to fail (ITF): the average number of iterations before an experiment
fails due to loss of control or track boundary violation. We set the maximum
number of iterations in an experiment to 20 laps around the track.

###  V-A Robustness Study on Learning Parameters

This experiment aims to study how LMPC parameter tunings affect the safety of
the learning process. The motivation is that there are not many trial-and-
error opportunities when tuning parameters on safety-critical hardware such as
a race car, and a robust control setup should be able to handle a broad range
of parameter settings without destabilizing the system. We set up simulation
experiments to compare our approach against the LMPC implementation in [8] in
terms of their robustness to the change in the control-rate cost (CRC)
($c_{\Delta u}$ in (4a)). Intuitively the CRC is analogous to the “learning
rate” for LMPC. A higher CRC will penalize drastic changes in the control
input and keep the vehicle closer to the safe set over the MPC horizon. It is
therefore a key parameter when tuning the behavior of the LMPC policy.

In simulation, we ran both LMPC implementations with 5 different CRCs across
two orders of magnitude and recorded the lap time and average velocity at
every iteration. The simulation model is formulated to be as close as possible
to the 1/10th scale vehicle used in the following hardware experiment. For
fairness of comparison, we purposely introduced a mismatch in the tire-road
friction coefficient $\mu$ in our nominal dynamics model. Whereas a value of
0.9 is used in the simulation model, a value of 1.2 is chosen for the LMPC
nominal dynamics model in (1), which is typically only seen on full scale race
cars with racing slick tires. This would have induced dangerous over-
confidence in the model about the surface friction, as shown in the following
hardware experiments, if no dynamics learning is used.

Table I presents the key results of this experiment. We swept CRC from 1.0 to
0.01 and compared the performance of the two LMPC implementations. With a high
CRC penalty (experiments 1 and 2), both [8] and ours fail to reach the optimal
lap time in 20 iterations, with [8] achieving a more significant lap time
reduction. However, once we decrease the CRC to 0.1, [8] starts to fail at
around lap 8, whereas our method exhibits greater robustness to different
settings of the CRC and can converge to the optimal lap time within 20
iterations once the CRC drops below 0.1. Figure 3 compares the per-iteration
performance of both methods, showing that [8] fails with more than half of the
CRC settings, while our approach was able to keep the system safe for all
settings.

To understand how the accuracy of the learned model compares between the two
approaches, we inspect the average norm of the difference between the learned
dynamics matrix $A(\bar{z},\mathcal{D}^{j-1})$ and the actual local dynamics
matrix of the simulation model, which is computed by taking the Jacobian of
the true system model at the same points $\bar{z}$. Figure 2 shows the
evolution of the mismatch in the dynamics matrix over laps of a simulation
experiment with a CRC value of 1. We highlight that although the modeling
error of full dynamics regression in [8] eventually reaches the same level as
our error dynamics regression, our learning approach better approximates the
actual local dynamics consistently and always only shows a small model
mismatch after learning. This is especially important in the initial
iterations when there is relatively scarce data to accurately fit an ATV
model. We believe that this large initial model mismatch with [8] may have
contributed to the failure cases with lower settings of CRC, as the LMPC is
unintentionally reaching the handling limit due to exploiting a poorly fitted
model before sufficient data can be collected to perform an accurate
regression of the full dynamics.

Figure 2: Comparison of average learned model error per iteration for a
simulation experiment.

Overall, our simulation experiments suggest that LMPC with error dynamics
regression exhibits greater stability and robustness to key parameter settings
and works well even when scarce dynamics data is available. These properties
could potentially make the system more suitable for hardware deployment in
safety-critical applications.

Experiment # | Control Rate Cost | ILT-20 | ITF
---|---|---|---
[8] | ours | [8] | ours
1 | 1.0 | 8.5 | 11.8 | 20+ | 20+
2 | 0.5 | 6.6 | 9.5 | 20+ | 20+
3 | 0.1 | - | 6 | 9 | 20+
4 | 0.05 | - | 4.9 | 8 | 20+
5 | 0.01 | - | 4.8 | 4 | 20+
TABLE I: Results of the robustness study. An ITF of 20+ means that the
algorithm did not fail in the 20-iteration experiment. Figure 3: Per-iteration
results of the robustness study. The left column shows [8] and the right
column shows our lap time and average speed with different control-rate cost
(CRC) tunings. A red cross is placed where failure occurs due to unsafe
handling of the vehicle by the LMPC.

###  V-B Hardware Experiment Comparing Learning Performance

In this section, we present results from physical experiments with the
Berkeley Autonomous Race Car (BARC) platform. These robotic vehicles are based
on 1/10 scale, off-the-shelf RC cars that have been heavily modified to
support autonomous vehicle architecture. Their specifications are listed in
Table II. LMPC computation is done off-board on a laptop computer and
localization is provided by a motion capture system.

Specification | BARC | IAC Car
---|---|---
Mass | 2.2 $\mathrm{kg}$ | 800 $\mathrm{kg}$
Wheelbase | 0.26 $\mathrm{m}$ | 4 $\mathrm{m}$
Top Speed |  $\sim$10 $\mathrm{km}\text{/}\mathrm{h}$ |  $\sim$320 $\mathrm{km}\text{/}\mathrm{h}$
Top Lateral Acceleration |  $\sim$10-15 $\mathrm{m}\text{/}{\mathrm{s}}^{2}$ |  $\sim$15-20 $\mathrm{m}\text{/}{\mathrm{s}}^{2}$
Top Longitudinal Acceleration |  $\sim$10-15 $\mathrm{m}\text{/}{\mathrm{s}}^{2}$ |  $\sim$10-15 $\mathrm{m}\text{/}{\mathrm{s}}^{2}$
Top Braking |  $\sim$10-15 $\mathrm{m}\text{/}{\mathrm{s}}^{2}$ |  $\sim$20-30 $\mathrm{m}\text{/}{\mathrm{s}}^{2}$
Max Engine Power |  $\sim$100 $\mathrm{W}$ |  $\sim$340 $\mathrm{kW}$
TABLE II: Vehicle specifications of BARC and IAC platforms used in this work.

In total, we conducted four sets of experiments. In the first two sets, the
LMPC proposed in [8] is run on the BARC platform without any regressive
learning of the vehicle dynamics. In the first case, we purposefully under-
estimate the vehicle handling limit by reducing the tire-road friction
coefficient to $\mu=0.7$ in the nominal model. In the second case, we over-
estimate by setting the coefficient to $\mu=1.2$. Note that the actual tire-
road friction coefficient was measured to be 0.9. The last two sets of
experiments are done with the full dynamics learning from [8] and our error
dynamics learning strategy, respectively. Each set of experiments consists of
10 runs under constant LMPC tuning with $\text{CRC}=0.1$, each lasting at most
20 iterations. We record the lap times and record any failure cases before the
20-iteration threshold.

Table III compares the metrics in these four sets of experiments. When the
LMPC is run with $\mu=0.7$, the car is no longer able to reach the optimal lap
time, presumably because the mismatched dynamics lead the controller to under-
estimate the handling limit. With $\mu=1.2$, although the LMPC is still able
to achieve the same lap time as our work, it fails during several trials after
17-19 iterations due to track constraint violation, exhibiting borderline
unsafe behaviors. We next compare the two regressive learning strategies.
Under this low CRC setting, LMPC with full dynamics regression [8] attempts to
improve the lap time quickly, but eventually fails after no more than 5
iterations in all 10 trials. On the other hand, our method converges to a best
lap time safely in all experiments, while exhibiting faster lap time reduction
when compared to the cases of no dynamics learning, as shown in Figure 4. We
note that although Figure 4 seems to suggest that simply over-estimating the
vehicle handling limit could also lead to optimal lap times, it comes at the
cost of a lower ITF value and is potentially more dangerous. Figure 5
visualizes the vehicle’s trajectory over iterations with our regressive
learning strategy, which shows convergence towards an optimal race line.

LMPC Configuration | Avg. ILT-20 | Avg. ITF
---|---|---
No regression ($\mu=0.7$) | 6.292 | 20+
No regression ($\mu=1.2$) | 5.42 | 18.9
Full Dynamics Regression ([8]) | - | 4.2
Error Dynamics Regression (Ours) | 5.51 | 20+
TABLE III: Key metrics of the BARC hardware experiment. A hyphen in the ILT-20
column means that the algorithm fails before the 20th iteration. An ITF of 20+
means that the algorithm did not fail over the 20 laps. Figure 4: Per-
iteration results of the BARC hardware experiment. The four smaller plots on
the left show the lap-time reduction of the 10 trials. The combined plot on
the right compares the average lap time achieved. A red cross is placed where
failure occurs due to unsafe handling of the vehicle by the LMPC. Figure 5:
Visualization of the 1st, 5th and 20th iteration’s trajectory of LMPC with
error dynamics regression on the BARC platform.

###  V-C Hardware Demonstration on a Full-Size Race Car

In this demonstration, we deployed our LMPC and error dynamics regression
strategy on a full-size IndyLights race car which has been used in the Indy
Autonomous Challenge [3]. The key specifications of the vehicle are listed in
Table II. The race track used for this demonstration is the Putnam Park Road
Course in Indiana, USA. We construct the initial dataset $\mathcal{D}^{0}$ for
the vehicle by running a tracking MPC on the center line of the track at 10
$\mathrm{m}\text{/}\mathrm{s}$ for 3 laps. We note that this step could also
be done with a more rudimentary controller. We then ran LMPC with error
dynamics regression for 10 iterations. Figure 1 visualizes the trajectory
driven by the vehicle, with a zoomed-in view to highlight how LMPC optimizes
the shape of the trajectory for a hairpin turn. It can be clearly seen in
Figure 1 that the vehicle starts at lower speeds closer to the center line and
incrementally increases its speed through the corners and approaches the track
boundary constraint.

##  VI Conclusion

This work proposes a new error dynamics learning approach under the LMPC
framework to iteratively and safely explore the boundary of vehicle handling
limit in autonomos racing. We based the learning on a nominal dynamics model,
which provides explanability and safety fallback for the learning process. We
then derived the regression method on top to learn a local approximation of
the nonlinear model mismatch. Through simulation and hardware experiments, we
showed that our method exhibits greater robustness against parameter tuning
and data scarcity compared with previous LMPC works, and offers new potential
solutions to address the challenge of vehicle dynamics data acquisition in
high-speed autonomous racing.

## References

  * [1] M. O’Kelly, H. Zheng, D. Karthik, and R. Mangharam, “F1TENTH: An Open-source Evaluation Environment for Continuous Control and Reinforcement Learning,” in _Proceedings of the NeurIPS 2019 Competition and Demonstration Track_. PMLR, Aug. 2020, pp. 77–89.
  * [2] “Indy Autonomous Challenge,” Feb. 2023. [Online]. Available: https://www.indyautonomouschallenge.com
  * [3] A. Wischnewski, M. Geisslinger, J. Betz, T. Betz, F. Fent, A. Heilmeier, L. Hermansdorfer, T. Herrmann, S. Huch, P. Karle, F. Nobis, L. Ögretmen, M. Rowold, F. Sauerbeck, T. Stahl, R. Trauth, M. Lienkamp, and B. Lohmann, “Indy autonomous challenge – autonomous race cars at the handling limits,” 2022.
  * [4] A. Wischnewski, T. Herrmann, F. Werner, and B. Lohmann, “A Tube-MPC Approach to Autonomous Multi-Vehicle Racing on High-Speed Ovals,” _IEEE Transactions on Intelligent Vehicles_ , vol. 8, no. 1, pp. 368–378, Jan. 2023.
  * [5] J. Kabzan, M. I. Valls, V. J. F. Reijgwart, H. F. C. Hendrikx, C. Ehmke, M. Prajapat, A. Bühler, N. Gosala, M. Gupta, R. Sivanesan, A. Dhall, E. Chisari, N. Karnchanachari, S. Brits, M. Dangel, I. Sa, R. Dubé, A. Gawel, M. Pfeiffer, A. Liniger, J. Lygeros, and R. Siegwart, “AMZ Driverless: The full autonomous racing system,” _Journal of Field Robotics_ , vol. 37, no. 7, pp. 1267–1294, 2020.
  * [6] J. Betz, T. Betz, F. Fent, M. Geisslinger, A. Heilmeier, L. Hermansdorfer, T. Herrmann, S. Huch, P. Karle, M. Lienkamp, B. Lohmann, F. Nobis, L. Ögretmen, M. Rowold, F. Sauerbeck, T. Stahl, R. Trauth, F. Werner, and A. Wischnewski, “TUM autonomous motorsport: An autonomous racing software for the Indy Autonomous Challenge,” _Journal of Field Robotics_ , vol. 40, no. 4, pp. 783–809, 2023.
  * [7] J. Funke, P. Theodosis, R. Hindiyeh, G. Stanek, K. Kritatakirana, C. Gerdes, D. Langer, M. Hernandez, B. Müller-Bessler, and B. Huhnke, “Up to the limits: Autonomous Audi TTS,” in _2012 IEEE Intelligent Vehicles Symposium_ , Jun. 2012, pp. 541–547.
  * [8] U. Rosolia and F. Borrelli, “Learning How to Autonomously Race a Car: A Predictive Control Approach,” _IEEE Transactions on Control Systems Technology_ , vol. 28, no. 6, pp. 2713–2719, Nov. 2020.
  * [9] M. Bujarbaruah, X. Zhang, U. Rosolia, and F. Borrelli, “Adaptive mpc for iterative tasks,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 6322–6327.
  * [10] A. Sasfi, M. N. Zeilinger, and J. Köhler, “Robust adaptive mpc using control contraction metrics,” _Automatica_ , vol. 155, p. 111169, 2023.
  * [11] J. Kabzan, L. Hewing, A. Liniger, and M. N. Zeilinger, “Learning-based model predictive control for autonomous racing,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 3363–3370, 2019.
  * [12] L. Hewing, A. Liniger, and M. N. Zeilinger, “Cautious NMPC with Gaussian Process Dynamics for Autonomous Miniature Race Cars,” in _2018 European Control Conference (ECC)_ , Jun. 2018, pp. 1341–1348.
  * [13] J. Ning and M. Behl, “Vehicle dynamics modeling for autonomous racing using gaussian processes,” 2023.
  * [14] D. Kalaria, Q. Lin, and J. M. Dolan, “Delay-aware robust control for safe autonomous driving,” in _2022 IEEE Intelligent Vehicles Symposium (IV)_ , 2022, pp. 1565–1571.
  * [15] E. L. Zhu, F. L. Busch, J. Johnson, and F. Borrelli, “A gaussian process model for opponent prediction in autonomous racing,” 2023.
  * [16] T. Brüdigam, A. Capone, S. Hirche, D. Wollherr, and M. Leibold, “Gaussian process-based stochastic model predictive control for overtaking in autonomous racing,” 2021.
  * [17] D. Bristow, M. Tharayil, and A. Alleyne, “A survey of iterative learning control,” _IEEE Control Systems Magazine_ , vol. 26, no. 3, pp. 96–114, 2006.
  * [18] U. Rosolia and F. Borrelli, “Learning Model Predictive Control for Iterative Tasks: A Computationally Efficient Approach for Linear System,” _IFAC-PapersOnLine_ , vol. 50, no. 1, pp. 3142–3147, Jul. 2017.
  * [19] M. Brunner, U. Rosolia, J. Gonzales, and F. Borrelli, “Repetitive learning model predictive control: An autonomous racing example,” in _2017 IEEE 56th annual conference on decision and control (CDC)_. IEEE, 2017, pp. 2545–2550.
  * [20] J. Kong, M. Pfeiffer, G. Schildbach, and F. Borrelli, “Kinematic and dynamic vehicle models for autonomous driving control design,” _IEEE Intelligent Vehicles Symposium, Proceedings_ , vol. 2015-August, pp. 1094–1099, 2015.
  * [21] F. Christ, A. Wischnewski, A. Heilmeier, and B. Lohmann, “Time-optimal trajectory planning for a race car considering variable tyre-road friction coefficients,” _Vehicle System Dynamics_ , vol. 59, no. 4, pp. 588–612, Apr. 2021.
  * [22] U. Rosolia and F. Borrelli, “Learning model predictive control for iterative tasks. a data-driven control framework,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 7, pp. 1883–1896, 2017.
  * [23] V. A. Epanechnikov, “Non-Parametric Estimation of a Multivariate Probability Density,” _Theory of Probability & Its Applications_, vol. 14, no. 1, pp. 153–158, Jan. 1969.