# On the status of the star Schulte 12 in the association Cyg OB2

V.G. Klochkova<EMAIL_ADDRESS>E.S. Islentieva  V.E. Panchuk
Special Astrophysical Observatory, Nizhnij Arkhyz, 369167 Russia

###### Abstract

High-resolution spectra of the LBV candidate Schulte 12 in the Cyg OB2
association were obtained at the 6-meter BTA telescope with the NES echelle
spectrograph on the arbitrary dates in 2001–2022. Variability of the emission
profile of H$\alpha$ and absorptions of He I, Si II with time was found. Based
on the radial velocity measurements at 10 observation dates, radial velocity
variability with an amplitude of $\Delta$Vr$\approx$8 km/s relative to the
average value of the heliocentric velocity Vr=$-15.6\pm 2.6$ km/s was
revealed. This indicates the presence of a companion in the system. Based on
the reliable intensity measurements of a sample of DIBs, color excess
E(B-V)=1.74$\pm 0.03^{m}$ was determined. This results in the interstellar
extinction value Av$\approx 5.6^{m}$ that is only about half of the total
extinction. Taking current Schulte 12 parameters, including Gaia EDR3
parallax, we estimated its absolute magnitude as Mv$\approx-9.2^{m}$ and
luminosity log$(L/L_{\odot})\approx 5.5$, which does not exceed the
Humphreys–Davidson limit.
Keywords: LBV candidates, Schulte 12, optical spectroscopy

##  1 Introduction

Hot high-luminosity star Schulte 12 (numbered according to Schulte (1958)), a
descendant of a far evolved massive star, is a member of the Cyg OB2
association. Being quite close (according to the Gaia EDR3 data, its parallax
$\pi=0.5895\pm 0.0518$ mas), this star has been studied quite well – the
SIMBAD database contains over 390 publications on it. Cyg OB2 is known as the
richest group in the Galaxy in terms of the number of massive OB stars Kiminki
et al. (2007); Massey and Thompson (1991). Insignificant distance to the
association ($\approx$1.6 Kpc) makes it possible to study in detail evolved
stars located there, which is important for the studies of the peculiarities
of the evolution of massive stars. Schulte 12 (=MT 304 in the list Massey and
Thompson (1991)) has long been considered as the brightest star in the Galaxy
(see paper Humphreys (1978) and references therein). Considering high
luminosity of the star in combination with its photometric and spectral
instability, authors Massey and Thompson (1991) classified it as a forming LBV
star. Clark et al. Clark et al. (2012), using data of spectral and photometric
monitoring of Schulte 12, estimated its parameters and concluded that the
luminosity of this star exceeds the Humphreys–Davidson limit. During the
following decade, a lot of studies were devoted to solution of this mysterious
phenomenon. The main method of the studies was long-term photometric and
spectral monitoring. The results and conclusions of different authors are
sometimes contradictory. For instance, Naze Naze et al. (2019), having
analyzed the data of long-term observations of Schulte 12 in the X-ray and in
the optical ranges, found possible values of the period of variability and
calculated preliminary parameters of the binary system. However, based on the
Gaia DR2 parallax, these authors concluded that in the Hertzsprung-Russell
diagram Schulte 12 is located far from the Humphreys–Davidson limit and, more
likely, is not a hypergiant, but a normal supergiant.

Chentsov et al. (2013) accomplished spectroscopy of Schulte 12 in 2001–2011
and found the signs of the radial velocity (hereafter – Vr) variability of
this star. In the next decade, the search for Vr variability continued (see
papers Mahy et al. (2022); Oskinova et al. (2017) and references therein).
However, being not so distant, Schulte 12 is weak in the optical region due to
extremely high extinction Maryeva et al. (2016), which complicates its high-
resolution spectral monitoring. The available publications still are based on
a restricted samples of spectra Klochkova and Chentsov (2004); Clark et al.
(2012); Naze et al. (2019); Chentsov et al. (2013). A full-fledged spectral
monitoring of Schulte 12, which is necessary to determine the cause of Vr
variability and to find the parameters of variability, has not yet been
carried out. This stimulated long-term observations with the 6-m BTA
telescope, and in this paper we present the results obtained on the basis of
Schulte 12 spectroscopy in 2001–2022 with high spectral resolution in the
250–300 nm wide wavelength range.

The main aim of our work is to search for the variability of spectral features
and the behavior of the radial velocity pattern with time. The methods of
observation and processing of spectral material are briefly described in
Section 2, its results are presented in Section 3. Their discussion and
conclusions follow in Section 4.

##  2 ECHELLE SPECTROSCOPY ON BTA

The spectra of Schulte 12 were obtained with the NES echelle spectrograph
Panchuk et al. (2017) stationary located in the Nasmyth focus of the 6-m BTA
telescope. The dates of observations of the star in 2001–2022 are listed in
Table 1. In the recent years, the spectrograph has been equipped with a CCD
array with the number of elements 4608$\times$2048; element size is
0.0135$\times$0.0135 mm; read noise is 1.8 e-. The wavelength range recorded
per exposure is $\Delta\lambda$=470–778 nm. In 2001–2011, we used a CCD with
2048$\times$2048 elements.

To reduce the flux loss at the entrance slit, the NES spectrograph is equipped
with a star image slicer. Using a slicer, each spectral order is repeated
three times. The spectral resolution of the NES is
R=$\lambda/\Delta\lambda\geq$60000\. In the spectra of Schulte 12, the signal-
to-noise ratio along the order varies from 10 to 40. In this work, we also
used the spectrum of Schulte 12, obtained in 2001 with the PFES echelle
spectrograph in the primary focus of the BTA, which provides a moderate
resolution. A detailed description of this spectrograph is presented in
Panchuk et al. (1997).

Table 1: Results of the measurements of the heliocentric velocity Vr in the spectra of Schulte 12 for different types of spectral details. Numbers in the parentheses indicate the number of spectral features used in the averaging.  Date | $\Delta\lambda$ | Vr, km/s
---|---|---
JD | nm | absorptions | absorptions “main” | DIBs
1 | 2 | 3 | 4 | 5
12.06.2001 | 454–793 | $-10.1\pm 2.0$ | $-13.5\pm 1.4$ | $-9.0\pm 0.5$
2452073.48 |  | (28) | (11) | (32)
12.04.2003 | 528–676 | $-24.0\pm 0.6$ | $-22.4\pm 0.6$ | $-10.7\pm 0.4$
2452742.47 |  | (24) | (12) | (18)
08.12.2006 | 447–594 | $-22.7\pm 1.1$ | $-22.1\pm 1.3$ | $-9.5\pm 1.2$
2454078.27 |  | (17) | (10) | (19)
26.09.2010 | 522–669 | $-9.4\pm 0.8$ | $-8.4\pm 1.4$ | $-9.4\pm 0.6$
2455466.32 |  | (30) | (6) | (30)
19.11.2010 | 522–669 | $-10.5\pm 0.7$ | $-11.0\pm 0.8$ | $-10.0\pm 0.3$
2455520.28 |  | (27) | (9) | (16)
11.08.2017 | 470–778 | $-2.6\pm 0.6$ | $-1.84\pm 0.8$ | $-10.0\pm 0.5$
2457976.51 |  | (47) | (9) | (45)
09.12.2019 | 470–778 | $-21.3\pm 0.6$ | $-19.7\pm 1.2$ | $-10.4\pm 0.4$
2458826.22 |  | (38) | (11) | (45)
07.08.2020 | 470–778 | $-20.0\pm 1.2$ | $-23.5\pm 1.4$ | $-10.7\pm 0.2$
2459068.52 |  | (34) | (10) | (38)
26.10.2021 | 470–778 | $-10.1\pm 2.4$ | $-8.1\pm 1.9$ | $-10.8\pm 0.4$
2459514.24 |  | (26) | (6) | (41)
11.06.2022 | 470–778 | $-25.0\pm 0.8$ | $-25.4\pm 1.1$ | $-10.2\pm 0.2$
2459742.44 |  | (47) | (11) | (15)

Extraction of one-dimensional data from two-dimensional echelle spectra was
performed using the ECHELLE context of the MIDAS package, modified considering
the geometry of our echelle frame. All details of the procedure are described
by Yushkin and Klochkova (2005). The traces of cosmic particles have been
removed by the standard method, by means of median averaging of a pair of
spectra obtained sequentially. A Th–Ar lamp was used to calibrate the
wavelengths. All subsequent steps in the processing of one-dimensional spectra
were performed using the modern version of the DECH 20t package developed by
G. Galazutdinov. The systematic error of the measurements of the heliocentric
radial velocity Vr from a set of telluric features and interstellar lines of
the Na I doublet does not exceed 0.25 km/s along one line; the measurement
error of Vr for broad absorptions does not exceed 0.5 km/s. For the average
values of the velocity in Table 1, the errors are 0.6–2.4 km/s, depending on
the type of the measured lines. To identify details in the spectra of Schulte
12, we used the lists of lines from Klochkova and Chentsov (2004). In
addition, the studies of the spectroscopy of related objects from Chentsov et
al. (2003); Chentsov and Sarkisyan (2007) were used.

##  3 RESULTS

###  3.1 The profiles of H$\alpha$ and other lines

The optical spectrum of Schulte 12 is not particularly impressive compared to
the spectra of other hot evolved high–luminosity stars. For comparison, we use
the results of spectroscopy of MWC 314, which is close to Schulte 12 over a
number of parameters. MWC 314 has an extremely high luminosity and a
temperature close to Schulte 12, has no significant photometric variability,
and is therefore considered an LBV candidate. But unlike Schulte 12, the
optical spectrum of MWC 314 contains a rich variety of features expected in
the spectrum of an LBV star: powerful H I and He I emissions, permitted (Si
II, Fe II) and forbidden emissions of [N II], [Ca II], [Fe II] and of other
metals, often with a two-peak profile Miroshnichenko et al. (1998); Frasca et
al. (2016). The spectrum of Schulte 12 contains mainly absorptions of ions of
light elements, with the exception of the emission of H$\alpha$ with an
unusual and variable profile, which is clearly seen in Fig. 1. The type of the
H$\alpha$ profile and its character resembles the picture of the variability
of this emission in the spectra of the blue hypergiant HD 183143 (see Fig. 5
in Chentsov et al. (2003)). However, in the spectra of HD 183143, the lower
part of the pedestal of the H$\alpha$ profile is much narrower than in the
spectra of Schulte 12.

Figure 1: H$\alpha$ profiles in 2010–2021: September 26, 2010 (blue line),
November 8, 2017 (green), September 12, 2019 (red), and October 26, 2021
(black). Here and below in the profiles in subsequent figures, the telluric
details are not removed.

In addition to the H$\alpha$ emission profile, the spectra of Schulte 12
contain weak emissions at the wavelengths $\lambda\approx$7495 and 7513 Å.
Both details are identified as Fe II emissions; they are also present in the
spectra of the related hypergiant HD 183143 Chentsov et al. (2003). The
intensity of the second Fe II emission $\lambda$ 7513.17 Å in the spectrum of
Schulte 12 is higher, its equivalent width for the different dates of our
observations is Wλ=0.04–0.08 Å. Based on our six spectra with the Fe II
emission of $\lambda$ 7513 Å, including the spectrum with a moderate
resolution from June 12, 2001, we obtain the average value of the velocity
corresponding to the position of this emission: Vr(7513)=$-19.0\pm 3.1$ km/s.
The average value without considering the spectrum dated June 12, 2001,
Vr(7513)=$-21.7\pm 1.3$ km/s, differs little from the previous one, but has a
significantly lower error. The region of formation of this emission is
probably the gaseous circumstellar medium, therefore, the value Vsys=$-21.7\pm
1.3$ km/s can be taken as the systemic velocity for Schulte 12. This velocity
value does not coincide with the average Vsys$\approx-10.3$ km/s value for the
Cyg OB2 members presented in Kiminki et al. (2007), but on the velocity
histogram for the association it falls into the normal range of Vsys. Note
that the purely emission profile of the Mg II 10952 Å line is centered near
Vr$\approx-25$ km/s in the spectra of Schulte 12 shown in the fragment of Fig.
5 in Naze et al. (2019). These authors used the average position of two
emissions, Fe II 7513 Å and Mg II 10952 Å, to estimate the parameters of the
suspected binary system with a period of about 108 day.

Figure 1 shows significant variations in the shape and intensity only in the
upper part of the H$\alpha$ profile. The lower part of the profile with the
wings wider than 400 km/s looks like a symmetrical and stable pedestal. Figure
2 illustrates the H$\alpha$ profile variations in the Schulte 12 spectra for
different observation dates in more detail. Position of the vertical dashed
line in all fragments of this Figure corresponds to the velocity
Vr(7513)$\approx-21.7$ km/s, taken as the systemic one.

Figure 2: Variability of the H$\alpha$ profile in the spectra of Schulte 12.
Position of the vertical line corresponds to the velocity Vr(7513)=$-21.7$
km/s, we taken as the systemic one

###  3.2 The pattern of radial velocity

Results of our measurements of the heliocentric radial velocity Vr using lines
of various types in the spectra of Schulte 12 are summarized in Table 1. The
short-wavelength part of the spectra is unsuitable for measurements due to
extreme extinction; therefore, we measured mainly the details of the spectra
in the wavelength region $\lambda\geq$5000 Å. In addition, He I and Si II
lines were not included in the calculation of average values because of the
significant asymmetry and variability of their profiles. These features of the
He I lines are clearly seen in the example of the He I 5876 Å line in Fig. 3,
which shows the profiles for several observation dates. Position of the
absorption core of the line varies in the range from +30 to $-$38 km/s. Figure
4 shows a fragment of two spectra containing Si II absorptions $\lambda$ 6347
and 6371 Å for two observations dates. Here one can clearly see variability of
the shape of the profile of Si II lines and of the position of their cores, as
well as the elongation of the short-wavelength wing of their profiles due to
the effect of the wind. At the same time, as expected, the interstellar bands
in this Figure are stationary: neither their shape nor position change.

Table 2: List of absorptions “main” in the spectra of Schulte 12 used for the values of Vr in the column 4 of Table 1 $\lambda$, Å | Element |  $\lambda$, Å | Element
---|---|---|---
4793.65 | NII(20) | 5133.12 | CII(16)
4803.29 | NII(20) | 5139.17 | CII(16)
4895.11 | NII(1) | 5142.34 | SII(1)
4994.36 | NII(24.64) | 5143.49 | CII(16)
5001.14 | NII(19) | 5145.16 | CII(16)
5001.48 | NII(19) | 5151.09 | CII(16)
5002.70 | NII(4) | 5512.70 | OI(25)
5010.62 | NII(4) | 5526.25 | SII(11)
5025.66 | NII(19) | 5535.35 | CII(10)
5040.72 | NII(19) | 5573.47 | FeIII(68)
5041.03 | SiII(5) | 5639.97 | SII(14)
5045.10 | NII(4) | 5640.33 | SII(11)
5047.29 | SII(15) | 5645.67 | SII(6)
5055.96 | SiII(5) | 5666.63 | NII(3)
5056.31 | SiII(5) | 5676.02 | NII(3)
5063.46 | FeIII(5) | 5679.56 | NII(3)
5073.90 | FeIII(5) | 5686.21 | NII(3)
5086.72 | FeIII(5) | 5696.60 | AlIII(2)
5093.56 | FeII | 5710.77 | NII(3)
5097.27 | FeII | 5722.73 | AlIII(2)
5100.74 | FeII(35) | 5730.65 | NII(3)
5103.34 | SII(7) | 5739.73 | SiIII(4)
5121.82 | CII(12) | 5747.30 | NII(9)
5125.20 | CII | 6482.05 | NII(8)
5127.35 | FeIII(5) | 6578.05 | CII(2)
|  | 6582.88 | CII(2)

Column 3 in Table 1 contains the average value Vr(abs) for each date for all
identified symmetric absorptions in the corresponding spectrum. Column 4 shows
the average values of Vr for the sample of the most reliable absorptions
without visible peculiare features of the profiles (N II, C II, Si III, Al
III), which we will call “main” in this text (see Table 2). The average values
in columns 3 and 4 of the Table 1 vary significantly from date to date. The
average values from column 3 change with the amplitude $\Delta$Vr(abs)$\approx
7.8$ km/s relative to the average for all dates Vr(abs)=$-15.5\pm$2.6 km/s.
Based on the measurements of the positions of the “main” lines listed in
column 4, the date-averaged value Vr(main)= $-15.6\pm 2.6$ km/s differs little
from the average velocity over all symmetric absorptions Vr(abs), but has a
slightly larger amplitude of velocity variability $\Delta$Vr$\approx 8.1$
km/s. Thus, the radial velocity variability in the spectrum of Schulte 12 with
an amplitude of about 8 km/s is beyond doubt.

Figure 3: Variability of the profile of the HeI 5876 Å line in the spectra of
Schulte 12 for the dates April 13, 2003 (black thin line), December 8, 2006
(blue), November 19, 2010 (lilac), August 11, 2017 (violet), December 9, 2019
(red), October 26, 2021 (black thick line), November 6, 2022 (dark green).

The authors of a recent paper Mahy et al. (2022) studied the variability of
several LBV and LBV candidates, including Schulte 12. Using high-resolution
spectra obtained with the HERMES spectrograph and a cross-correlation method
based on the selected lines, they obtained a high value of binarity rate
($\approx 62\%$), among these objects. For Schulte 12, using several N II
lines, He I 5876 and Si II 6347 Å lines, they obtained a total variation of
velocity $\Delta$Vr(max)=24.2 km/s over 3658 days, which is also close to our
data presented in Table 1.

Figure 4: Fragments of the spectra of Schulte 12 containing absorptions of Si
II and DIBs obtained in 2019 (red line) and 2021 (black line).

###  3.3 Circumstellar and interstellar details

Schulte 12 has long been known as one of the most reddened stars. Back in
1954, based on photoelectric observations, Morgan et al. (1954) noted that the
value of the color index of Schulte 12, E(B-V)=$+3.01^{m}$, is the largest
among hot stars and results in the interstellar extinction value from 9 to
10m. The authors of publications Clark et al. (2012); Maryeva et al. (2016),
modeling the spectrophotometric data for Schulte 12, found an equally high
total absorption Av$\approx 10^{m}$.

Schulte 12 is also known as a star with an optical spectrum containing many
interstellar bands (DIBs), reliably identified by Klochkova and Chentsov
(2004). We emphasize that DIBs dominate among other absorptions in the
spectrum of the star. Two such bands are contained in the fragment in Fig. 4,
which also includes two absorptions of SiII 6347 and 6374 Å. In column 5 of
the Table 1 we present information on the mean radial velocities Vr(DIBs) for
each of the spectra we have. Velocity value averaged over all observation
dates, Vr(DIBs)=$-10.0\pm 0.2$ km/s, is in good agreement with the analogous
value measured in the spectra of 12 other members of the Cyg OB2 association
Chentsov et al. (2013).

Multicomponent profile of the Na I 5889 Å line in the spectra of Schulte 12 in
2019 and 2021 is shown in Fig. 5. Short vertical line in this Figure indicates
the position of K I 7699 Å interstellar line. Its position averaged over our
five spectra, Vr(KI)=$-9.5\pm 0.2$ km/s, is in good agreement with the average
velocity Vr(DIBs)=$-10.0\pm 0.2$ km/s. The shortest-wavelength absorption “1”
of the Na I D–lines profile, whose position (Vr=$-32.0\pm 0.4$ km/s) does not
change from the date to the date of observations, probably is formed in the
circumstellar envelope. The circumstellar origin of this feature is also
confirmed by the small width of its profile, about 10 km/s. Absorption with
close position, Vr$\approx-30$ km/s, was also recorded in the spectra of other
B stars in Cyg OB2 Chentsov et al. (2013). Taking into account the accepted
value of the systemic velocity Vsys=$-21.7\pm 1.3$ km/s, we can estimate the
envelope expansion velocity: Vexp$\approx$10 km/s. The next equally narrow
component “2” on the NaD line profiles, Vr$\approx-16$ km/s, may also have a
circumstellar origin. However, this still remains at the level of a
conjecture. Component “3” arises in the interstellar medium, where many DIBs
and the K I line are formed.

Figure 5: Multicomponent profile of the NaI 5889 Å line in the Schulte 12
spectra obtained in 2019 (red line) and 2021 (black line). Position of the
short vertical line corresponds to the velocity Vr=$-9.5$ km/s for the
interstellar line K I 7699 Å.

Earlier, Chentsov et al. (2013) used the BTA+NES spectra to study the
absorption for 13 hot stars in Cyg OB2, including Schulte 12. Using the
interstellar components of the Na I D-lines and the DIB 5797 interstellar
band, these authors confirmed the extremely high reddening of Schulte 12 and
assumed the contribution of the circumstellar medium to it, which is well
illustrated in Fig. 16 in the publication Chentsov et al. (2013). Maryeva et
al. (2016), having studied the extinction behavior in the Cyg OB2 association,
found an increase in the reddening in the direction of Schulte 12. These
authors also concluded that a significant proportion of the color excess is
due to the contribution of the circumstellar medium.

Table 3: Equivalent widths Wλ(DIBs) averaged over Schulte 12 spectra. In column 3, the values of color excess according to the calibrations by Luna et al (2008) are printed in boldface, while those according to the calibrations by Kos and Zwitter (2013) in regular font. $\lambda$, | Wλ, | E(B-V)
---|---|---
Å | mÅ | mag
5780.48 | 829 |  1.8; |  $>$1.4
5797.06 | 269 |  1.7; |  $>$1.4
6195.98 | 109 |  1.8; |  $>$1.4
6269.75 | 199 |  |  $>$1.4
6379.32 | 159 |  1.7; |  1.4
6445.20 | 58 |  |  $>$1.4
6613.62 | 366 |  1.7; |  $>$1.4
6660.71 | 56 |  |  1.4

With high-quality spectral material, we estimated the color excess of Schulte
12 corresponding to the equivalent widths of DIBs for which calibrations Wλ
depending on E(B-V) have been published. In Table 3, we present for selected
DIBs and corresponding color excesses using the calibrations from Luna et al.
(2008); Kos and Zwitter (2013). Due to the too high interstellar extinction
toward Schulte 12, we unfortunately could not obtain quantitative estimates
using wellknown calibrations by Kos and Zwitter (2013), since the equivalent
widths of DIBs in the spectrum of Schulte 12 are far beyond the calibration
limit in this work. According to the calibrations by Luna et al. (2008), the
average color excess calibrations [21], the average color excess
E(B-V)=1.74$\pm 0.03^{m}$. Applying the standard ratio of extinction to color
excess, R=3.1, we obtain the interstellar extinction for Schulte 12:
Av(DIBs)$\approx 5.6^{m}$. It follows that about half of the extreme
extinction of Schulte 12 is due to the existence of a powerful circumstellar
envelope.

Taking the current value of the distance to Schulte 12 kpc from the catalog
Bailer-Jones (2021) based on the Gaia EDR3 data, effective temperature 13.7 kK
Clark et al. (2012) and the total extinction from Maryeva et al. (2016), we
estimated the absolute magnitude of this star as Mv$\approx-9.2\pm 0.15^{m}$
and its luminosity log$(L/L_{\odot})\approx$5.5. Note that the parallax
inaccuracy caused by the presence of a powerful circumstellar envelope, which
creates the effect of a pseudophotosphere, increasing the apparent angular
diameter of the star and distorting its parallax, can also introduce
uncertainty in the value of the absolute magnitude of Schulte 12. This effect
is known for high-luminosity stars with envelopes. For Schulte 12, it was
noted in papers Naze et al. (2019); Chentsov (2021).

##  4 DISCUSSION OF RESULTS AND CONCLUSIONS

It is known that, along to the extremely high luminosity, LBV stars are
distinguished by their significant photometric and spectral variability. The
main criterion that ensures reliable identification of an LBV star among other
massive evolved stars is the fixation of a cardinal change in its brightness
and color. An excellent illustration of the specific behavior of the
brightness of LBV stars can be found in the review by van Genderen (2001),
where Fig. 11 shows the light curve for $\eta$ Car over the period of its
observations between 1600 and 2000.

As it follows from the ASAS-SN data Kochanek et al. (2017), during the last
1200 nights of observations, Schulte 12 did not experience significant changes
in brightness relative to the mean value V=11.6m. Earlier photometric data for
the last decade of the 20th century, with variability of the same level, are
presented by Clark et al. (2012). Strict periodicity has not yet been revealed
in the long-term series of photometric data for Schulte 12. The observed
irregular brightness instability of Schulte 12 at the level of 0.1m is similar
to the microvariability in LBV stars Lamers et al. (1998), and can be
explained as the manifestation of wind instability. The absence of LBV
eruptions may indicate that Schulte 12 is in a “dormant” period. Such an
inactive period in the life of an LBV can be long; for example, a giant
eruption of the LBV candidate P Cyg took place around 1600. The presence of a
powerful circumstellar envelope around Schulte 12 indicates that this star
also has already experienced episodes of a large loss of matter.

Using measurements of DIBs equivalent widths in our Schulte 12 spectra and
calibration by Luna et al. (2008), we estimated the average color excess:
E(B-V)=1.74$\pm 0.03^{m}$. Applying the standard extinction to color excess
ratio, R=3.1, we obtained the value of the interstellar extinction for Schulte
12: Av(DIBs)$\approx 5.6^{m}$. Thus, we conclude that about half of the
extreme extinction of Schulte 12 is due to the existence of a powerful
circumstellar envelope

The most important physical parameter of an LBV candidate is the value of the
absolute luminosity. The problem of the extremely high luminosity of Schulte
12 has been considered by many authors. First of all, we note the already
mentioned work Clark et al. (2012), the authors of which simulated photometric
(from UV to radio range) and spectral data for Schulte 12 and obtained extreme
values of its parameters: absolute magnitude Mv$\approx-9.85^{m}$, luminosity
log$(L/L_{\odot})$=6.28 and mass $\mathcal{M/M}_{\odot}\approx$110\. As a
result, Clark et al. (2012) came to the conclusion that it is impossible,
within the existing concepts, to establish the evolutionary status for such a
massive and cold hypergiant.

In the summary of physical parameters for a sample of LBV and LBV candidates,
Mahy et al. (2022) give two luminosities for Schulte 12: the first value
log$(L/L_{\odot})$=6.1–6.3 is based on the Gaia EDR3 parallax and the second
one is based on the Gaia 2 data according to Smith et al. (2019), which is
significantly lower: log$(L/L_{\odot})$=5.7–5.9. Our luminosity estimate based
on the distance from the catalog Bailer-Jones (2021) based on Gaia EDR3, the
modern values of the effective temperature 13.7 kK Clark et al. (2012) and
total extinction Av$\approx 10^{m}$ Maryeva et al. (2016) is even lower:
log$(L/L_{\odot})$=5.5. Thus, the set of Schulte 12 parameters leads to its
normal position in the Hertzsprung–Russell diagram. To illustrate this
important conclusion, one can use, for example, the diagram for massive
evolved stars in the review by Kraus and Torres (2021).

One of our new results is evaluation of the expansion rate of the
circumstellar gaseous medium. The shortest-wavelength absorption of the Na I
D-lines profile, whose position Vr=$-32.0\pm 0.4$ km/s does not change between
the dates of observations, is probably formed in the circumstellar envelope.
Considering the accepted value of the systemic velocity Vsys=$-21.7\pm$1.3
km/s, we can estimate the value of the envelope expansion velocity:
Vexp$\approx$10 km/s. It can be assumed that the next, also stable component,
Vr$\approx-16.0$ km/s, is also formed in the circumstellar medium of Schulte
12.

Thanks to numerous observations, we recorded the variability of the H$\alpha$
emission profile, the reason for which still remains unclear. Similar changes
in the H$\alpha$ profile were recorded by Clark et al. (2012) in the spectra
of Schulte 12 in 1998, 2000, and 2008. To explain the variability of the
profile shape, these authors suggested that the complex H$\alpha$ profile may
be due to the presence of an additional short-wavelength nebular emission.
However, they stressed that such a conclusion now seems premature, since the
spectrum of Schulte 12 lacks the expected lines, primarily the forbidden [N
II], [S II] emissions characteristic for nebulae spectra.

The observed variability of the He I and Si II absorption profiles is due to
the influence of the wind. It should be noted here that this variability of
the He I 5876, Si II 6347 Å profiles makes an additional contribution to the
Schulte 12 velocity variability found by the authors Mahy et al. (2022).

We consider the reliable detection of radial velocity variability to be an
important result of our study: according to the spectra for 10 observation
dates, the amplitude of variability is $\Delta$Vr$\approx$8 km/s relative to
the average velocity Vr=$-15.6\pm$2.6 km/s, which indicates the presence of a
companion in the Schulte 12 system. The proof of the presence of companion is
a fundamental point in solving the problem of extremely high luminosity of
Schulte 12 at such a low effective temperature. This non-standard combination
of parameters leads to the deviation of the position of the star in the H–R
diagram from the theoretical isochrone for the members of the Cyg OB2
association (see Clark et al. (2012), Fig. 14). Recall, thar Maryeva et al.
(2016) discovered the companion of Schulte 12 by speckle interferometry at the
6-m telescope and estimated the orbital period of the binary system as 100–200
yr. Analyzing the X–ray spectrum of Schulte 12, Oskinova et al. (2017) also
came to the conclusion about the presence of an O-companion and indicated its
presence as the cause of the collision of winds in the binary system and the
appearance of X–ray. The authors of a current review by Weis and Bomans (2020)
made a more general conclusion that the binarity of a massive star may be a
valid indication that the star is observed on its way to the LBV stage.

The search for companions is an urgent task in the understanding of the nature
of stars with extreme masses and luminosities, as it is well illustrated by
recent speckle observations of stars in the nucleus of the R 136 cluster
Kalari et al. (2022). These observations with the Gemini telescope made it
possible, in particular, to resolve the components of R 136a1, the most
massive star known. The results by Kalari et al. (2022) made the estimates of
the masses of stars in the nucleus of R 136 by a factor of 1.5–2 lower.

To confirm our results and to estimate the parameters of the binary system, a
long-term continuation of the Schulte 12 spectral monitoring with high
spectral resolution and high signal–to–noise ratio is required, which will
provide the necessary series of homogeneous kinematic data. New spectral
observations will serve as the basis for the modeling of the profiles of
characteristic features in the spectrum of Schulte 12 and studying the
behavior of their parameters over time.

## Financing

We acknowledge financial support of the Russian Science Foundation (grant No.
22-12-00069). Observations by 6-m telescope of the Special Astrophysical
Observatory of RAS were supported by the Ministry of Science and Higher
Education of Russian Federation.

## Acknowledgements

This study made use of the astronomical data bases SIMBAD, VALD, SAO/NASA ADS,
and Gaia DR3.

## References

  * Schulte (1958) D.H. Schulte, Astrophys. J. 128 41 (1958).
  * Massey and Thompson (1991) P. Massey and A.B. Thompson, Astron. J. 101 1408 (1991).
  * Kiminki et al. (2007) D. C. Kiminki, H.A. Kobulnicky, K. Kinemuchi, et al., Astrophys. J. 664 1102 (2007).
  * Humphreys (1978) R.M. Humphreys, Astrophys. J. Suppl. 38 309 (1978).
  * Clark et al. (2012) J.S. Clark, F. Najarro, I. Negueruela, et al., Astron. Astrophys. 541 A145 (2012).
  * Naze et al. (2019) Y. Naze, G. Rauw, S. Czesla, L. Mahy, and F. Campos, Astron. Astrophys. 627 A99 (2019).
  * Chentsov et al. (2013) E.L. Chentsov, V.G. Klochkova, V.E. Panchuk, M.V. Yushkin, and D.S. Nasonov, Astron. Rep. 57 527 (2013).
  * Mahy et al. (2022) L. Mahy, C. Lanthermann, D. Hitsemekers, et al., Astron. Astrophys. 657 A4 (2022).
  * Oskinova et al. (2017) L.M. Oskinova, D.P. Huenemoerder, W.-R. Hamann, et al., Astrophys. J. 845 39 (2017).
  * Maryeva et al. (2016) O.V. Maryeva, E.L. Chentsov, V.P. Goranskij, et al., Mon. Not. R. Astron. Soc. 458 491 (2016).
  * Klochkova and Chentsov (2004) V.G. Klochkova and E.L. Chentsov, Astron. Rep. 48 1005 (2004).
  * Panchuk et al. (2017) V.E. Panchuk, V.G. Klochkova, and M.V. Yushkin, Astron. Rep. 61 820 (2017).
  * Panchuk et al. (1997) V.E. Panchuk, I.D. Najdenov, V.G. Klochkova, et al., Bull. Spec. Astrophys. Observ. 44 127 (1997).
  * Yushkin and Klochkova (2005) M.V. Yushkin and V.G. Klochkova, Preprint Spets. Astrofiz. Observ., No. 206 (2005).
  * Chentsov and Sarkisyan (2007) E.L. Chentsov and A.N. Sarkisyan, Astrophys. Bull. 62 257 (2007).
  * Chentsov et al. (2003) E.L. Chentsov, S.V. Ermakov, V.G. Klochkova, et al., Astron. Astrophys. 397 1035 (2003).
  * Miroshnichenko et al. (1998) A.S. Miroshnichenko, Y. Fremat, L. Houziaux, et al., Astron. Astrophys. Suppl. Ser. 131 469 (1998).
  * Frasca et al. (2016) A. Frasca, A.S. Miroshnichenko, C. Rossi, et al., Astron. Astrophys. 585, A60 (2016).
  * Morgan et al. (1954) W.W. Morgan, H.L. Johnson, and N.G. Roman, Publ. Astron. Soc. Pacif. 66 85 (1954).
  * Kos and Zwitter (2013) J. Kos and T. Zwitter, Astrophys. J. 774 72 (2013).
  * Luna et al. (2008) R. Luna, N.L.J. Cox, M.A. Satorre, et al. Astron. Astrophys. 480 133 (2008).
  * Bailer-Jones (2021) A.L. Bailer-Jones, J. Rybizki, M. Fouesneau, M. Demleitner, and R. Andrae, Astron. J. 161 147 (2021).
  * Chentsov (2021) E.L. Chentsov, Astrophys. Bull. 76 286 (2021).
  * van Genderen (2001) A.M. van Genderen, Astron. Astrophys. 366 508 (2001).
  * Kochanek et al. (2017) B.S. Kochanek, B.J. Shappee, K.Z. Stanek, et al., Publ. Astron. Soc. Pacif. 129 (980), (2017).
  * Lamers et al. (1998) H.J.G.L.M. Lamers, M.V. Bastiaanse, C. Aerts, and H.W.W. Spoon, Astron. Astrophys. 335 605 (1998).
  * Smith et al. (2019) N. Smith, M. Aghakhanloo, J.W. Murphy, et al. Mon. Not. R. Astron. Soc. 488 1760 (2019).
  * Kraus and Torres (2021) M. Kraus and A.F. Torres, Im: Pulsations along Stellar Evolution, Workshop Ser. Assoc. Argentina Astron. 12 234 (2021).
  * Weis and Bomans (2020) K. Weis and D.J. Bomans, Galaxies 8 20 (2020).
  * Kalari et al. (2022) V.M. Kalari, E.P. Horch, R. Salinas, et al., Astrophys. J. 935 162 (2022).


# Optimizing BioTac Simulation for Realistic Tactile Perception

Wadhah Zai El Amri  L3S Research Center, Leibniz Universität Hannover
Hannover, Germany
0000-0002-0238-4437 Nicolás Navarro-Guerrero  L3S Research Center, Leibniz
Universität Hannover
Hannover, Germany
0000-0003-1164-5579

###### Abstract

Tactile sensing presents a promising opportunity for enhancing the interaction
capabilities of today’s robots. BioTac is a commonly used tactile sensor that
enables robots to perceive and respond to physical tactile stimuli. However,
the sensor’s non-linearity poses challenges in simulating its behavior. In
this paper, we first investigate a BioTac simulation that uses temperature,
force, and contact point positions to predict the sensor outputs. We show that
training with BioTac temperature readings does not yield accurate sensor
output predictions during deployment. Consequently, we tested three
alternative models, i.e., an XGBoost regressor, a neural network, and a
transformer encoder. We train these models without temperature readings and
provide a detailed investigation of the window size of the input vectors. We
demonstrate that we achieve statistically significant improvements over the
baseline network. Furthermore, our results reveal that the XGBoost regressor
and transformer outperform traditional feed-forward neural networks in this
task. We make all our code and results available online on
https://github.com/wzaielamri/Optimizing_BioTac_Simulation.

###### Index Terms:

BioTac, XGBoost, Transformer, Tactile Perception

##  I Introduction

Tactile sensing sensors offer robots valuable information that can be used to
enhance and complement knowledge coming from other modalities such as vision
or audio, especially in situations where this knowledge is entirely or
partially not available [1, 2].

In situations where the sensor is unavailable or experiment repetitions are
costly, the value of a reliable, real-time simulation becomes evident. Such a
simulation can effectively estimate sensor outputs for various touch
scenarios. Such simulation would offer a good alternative to gathering data in
different setups and environments [1]. However, simulations only approximate
real-world data, and models trained on this data alone can typically not be
directly used in real robots. This problem is also known as the reality gap,
and it can be mitigated with techniques such as sim2real [3] or continual
learning [4]. Nevertheless, the quality of the simulated data should be as
high as possible to minimize the reality gap and increase the success of other
modules down the pipeline.

Various studies have presented simulations for a range of tactile sensors,
including TACTO [5] for vision-based sensors like DIGIT [6] and OmniTact [7],
touch simulation for fabric-based tactile sensors [8], a simulation for the
iCub skin [9], and several simulations for the BioTac sensor [10, 11, 12].

However, in this paper, we focus on the BioTac [13, 14], one of the most
widely used tactile sensors [15]. It consists of a rigid core covered by an
elastomeric skin filled with an incompressible conductive fluid. This sensor
proved to be helpful in different tasks, such as grasp stability [16], object
identification [17], localizing artificial tumors [18, 19], etc.

Some simulations of the BioTac sensor exist. Ruppel et al. [10] collected a
real-world dataset with the BioTac 2P sensor [20] mounted on a Shadow Robot
Hand [21]. The dataset captures the tactile sensor readings while touching
with different forces, positions, and orientations, an indenter with a
spherical tip of radius equal to $2\;mm$, attached to a calibrated force-
torque sensor (ATi nano17e [22]). They propose a neural network model that can
estimate electrodes, vibration, and pressure signals using the force and
position vector of the contact point alongside temperature values. However,
their approach requires using the temperature values in the input vector,
which are not available in rigid body simulation environments like Gazebo
[23].

Narang et al. [11] developed a finite element model (FEM) that can mimic the
deformation behavior of the elastomeric skin and the liquid gel inside. From
the point cloud of the FEM nodal field data, a neural network based on the
PointNet++ architecture [24] is used to estimate the electrode values of the
BioTac 2P [20]. To validate the FEM model and train the neural network, they
collected data on different trajectories using nine different indenters
interacting with three different tactile sensors. The FE simulation takes 7
minutes to simulate one single trajectory, which makes it a big drawback.
Later, they improved the approach using the NVIDIA Isaac Gym simulator [25],
achieving a 75 times improvement [26]. However, a single trajectory simulation
of $6\;mm$ takes $5.57\;s$ to simulate.

Furthermore, the authors improved the electrodes’ estimation, by using a self-
supervised latent learning network that maps between the FEM Mesh and
electrode signals. Despite these improvements, this solution is still not fast
enough to be used in real-time systems. Moreover, this solution estimates only
the electrodes’ values and not the pressure or vibration signals recorded by
the BioTac sensor, which are valuable modalities when used in grasping and
manipulation tasks [27].

Zaoata-Impata and Gil [12] used vision to estimate the electrode values of the
BioTac SP sensor [28]. A modified Semi-Regression Generative Adversarial
Network (SR-GAN) [29] was used to synthesize electrode values out of point
cloud representation of grasped deformable daily objects and their grasping
data by using labeled and unlabeled input data. The authors collected 4000
samples of grasping positions, electrodes, and pressure values of two BioTac
SP sensors mounted on a Shadow Robot Hand [21] and 3D point clouds recorded
with an Intel RealSense D415 depth camera. Since collecting data is not
trivial, a generator was used to generate synthetic BioTac data for the point
clouds and the grasping positions in order to augment the dataset and improve
the learning of the discriminator, whose task was to estimate the correct
BioTac outputs. The unlabeled samples boosted the performance of the network.
In addition, this work sheds light on the possibility of using vision, in this
case through point clouds, to estimate signal data for the BioTac sensor.
However, this implementation only used 4000 samples.

Since vibration and pressure readings add valuable information when used for
manipulation tasks [27], we must simulate these alongside electrode values.
The solution of Narang et al. [26] only estimates the electrode signals and is
comparatively slow. In comparison, Zaoata-Impata and Gil [12] implementation
uses a small dataset. Considering these factors, we opt for the solution
provided by Ruppel et al. [10].

In the next Section II, we revisit the research conducted by Ruppel et al.
[10]. Within that section, we first delve into the details of their work.
Subsequently, in Section III, we thoroughly investigate two directions to
improve the simulation further.

##  II BioTac Sensor Simulation

The simulation provided by Ruppel et al. [10] is based on a neural network
trained with a real-world dataset to predict the BioTac sensor signals. This
dataset includes data from approximately one hour of recording sampled at 100
Hz. It consists of the complete BioTac output values, i.e., 19 electrode
voltage $e_{1}$, …, $e_{19}$, absolute and dynamic fluid pressure pdc and pac,
temperature tdc and heat flow tac. The dataset also includes the BioTac’s and
indenter’s position and orientation recorded with a vision tracking system
[30]. For the data pre-processing pipeline, the position of the indenter’s tip
is calculated using transformation matrices in the BioTac coordinate frame,
represented in Figure 1.

Figure 1: (a): Map of the electrodes on the BioTac sensor [31]. (b): BioTac
sensor and the position/ orientation of the coordinate frame [31].

These transformation matrices are calibrated to adjust for systematic offsets
by reducing the error between the probes’ positions and the BioTac surface. In
other words, probes with light touch contact at the end of a contact cycle and
with force measurements below $0.3;N$ are gathered, and their distance to the
skin surface of the BioTac is calculated with the transformation matrices. In
such cases, that distance should be equal to 0. Due to calibration issues,
this distance deviates from 0, $5.039\;mm$ on average. While iterating for
1000 steps and adjusting these transformation matrices, i.e., adding and
subtracting small values, the distance error is reduced to $0.444\;mm$ on
average for all these selected probes.

Two deep neural networks, $A$ and $B$, were proposed to simulate the sensor
outputs at timestep $t=T$. Both networks take as input the $x$, $y$, and $z$
contact position of the indenter’s tip in the BioTac coordinate frame at
timestep $t=T$, $F_{x}$, $F_{y}$, and $F_{z}$ values of the force sensor at
timesteps $T$, $T-10$, and $T+10$ and the temperature directly taken from
BioTac readings at timestep $t=T$, i.e., an input vector of length 13.

Network $A$ is a dense neural network with five hidden layers and over 6
million parameters, employing a pre-mapping fusion strategy [1]. On the other
hand, Network $B$ is a more compact neural network with around 800 thousand
parameters and adopting a midst-mapping fusion strategy [1]. Network $B$
consists of three separate dense layer columns that merge into a dense layer.
Positional, force, and temperature values undergo separate processing within
these columns.

Both networks were trained with the loss function presented in Equation 1,
which consists of the summation of the mean absolute error (MAE) and the mean
squared error (MSE).

$\text{$Loss$}=\frac{\sum_{i=1}^{n}(|y_{i}-\hat{y}_{i}|+(y_{i}-\hat{y}_{i})^{2})}{n}$
(1)

The authors report achieving 9.3% normalized mean absolute error (MAE) over
all channels calculated on the z-score normalized outputs of the neural
network, i.e., not in the original scale between 0 and 4095, with Network $B$.
However, Network $A$ reached 11.3% MAE. In both cases, the reported results
values consider only one of the two dynamic pressure channels ($pac_{0}$)
because both channels are nearly identical. The reported results also omitted
the heat flow channel (tac). During the neural network training phase, both
channels, i.e., $pac_{1}$ and tac, were utilized without justification for
this approach, resulting in an output vector comprising 23 channels.

Following their findings, Network $B$ with separate columns yielded better
results and had less trainable parameters. Thus, we use it in all our
subsequent investigations. In particular, we identify two possible
improvements: one related to the temperature used as an input, since this
information is unavailable in simulators. The second area of improvement
concerns the windowing size used for the force and position values used in the
original paper, which was not sufficiently justified.

We use a similar split used by Ruppel et al.[10] to analyze these possible
improvements. The test and validation sets consist of 30 chunks of data, each
comprising 1000 data points. This configuration ensures that each set,
independently, accounts for 10% of the entire dataset. We train the network
without early stopping for 50 epochs. We also conduct 10-fold cross-validation
to enhance the statistical robustness of our metric reporting. We also do not
consider the channels $pac_{1}$ and tac in the metric calculation.

###  II-A Investigating the Training with Temperature Readings

The authors used the current temperature readings (tdc) at timestep $t=T$ of
the BioTac sensor to predict the output values. Since the used simulation
environment Gazebo [23] lacks temperature information, they suggest fixing
that input value to a specific constant, specifically the mean temperature
value of the entire collected dataset. However, the BioTac sensor is sensitive
to temperature changes [13], and plotting the temperature values of the
dataset, Figure 2, it can be seen that the temperature readings of the BioTac
sensor continuously increase over time. This data suggests that fixing the
temperature to a specific constant value will decrease the model’s
performance.

Figure 2: Temperature values of the BioTac sensor in the Ruppel et al. [10]’s
dataset. The dashed black line represents the average temperature value.

To test the effect of fixing the temperature value, we train Network $B$ using
the same hyperparameters used by the authors. However, during the testing
phase, we first compute the network’s outputs using the correct temperature
values provided in the dataset. Second, we compute the network’s output using
a fixed temperature value. Here, we test all temperature values in the
dataset, including the mean value suggested by Ruppel et al. [10]. The results
are shown in Figure 3.

Figure 3: Average normalized MAE over all channels for ten folds for Network
$B$. The solid line represents the normalized MAE values when fixing the
temperature input value and probing all temperatures in the dataset. The
dashed line represents the normalized MAE, when using the appropriate correct
temperature values of the test set. The shaded area represents the upper and
lower limits calculated by the standard deviation. The vertical dash-dotted
line represents the mean temperature of the entire dataset.

Examining Figure 3, we observe the lowest error when using the correct
temperature as input, $9.6\%$. However, fixing the temperature in the best-
case scenario leads to an averaged normalized MAE of $22.8\%$. Hence, fixing
the input temperature value translates to a relative performance loss of
$24.4\%$. We determine the relative loss between the two scenarios by
determining the MAE of a naive model ($53.7\%$), which always outputs the mean
value of each channel as the prediction. These findings emphasize the need for
an improved solution focusing on learning the data outputs without relying on
temperature information since these are unavailable in the simulation
environment.

###  II-B Investigating the Windowing Size

The authors reported using the input vector force readings at three timesteps:
$T-10$, $T$, and $T+10$, omitting all intermediate timesteps and claiming
their inclusion would lead to overfitting. Moreover, the use of force values
from future timesteps, i.e., $T+10$, remains unexplained, and error reduction
by including future values was not quantified. Thus, we conducted a thorough
analysis of the choice of windowing size for force values and explored the
impact of incorporating force values from future timesteps. To investigate
this, we re-train Network $B$ across six different conditions, keeping the
current temperature and current position and varying in each experiment the
windowing size of the force values as follows:

  1. 1.

Forces $F_{t}$: $\forall\text{ }t\in\\{T,T-10,T+10\\}$.

  2. 2.

Forces $F_{t}$: $\forall\;t\in\\{T,T-10\\}$.

  3. 3.

Forces $F_{t}$: $\forall\;t\in\\{T\\}$.

  4. 4.

Forces $F_{t}$: $\forall\text{ }t\in\\{T,T-10,T-5,T+5,T+10\\}$.

  5. 5.

Forces $F_{t}$: $\forall\text{ }t\in[T-10,T+10],\quad t\in\mathbb{N}$.

  6. 6.

Forces $F_{t}$: $\forall\text{ }t\in[T-10,T],\quad t\in\mathbb{N}$.

The training and validation loss curves for all conditions are depicted in
Figure 4. Additionally, Table I illustrates these six conditions’ normalized
mean absolute errors.

Figure 4: Visualization of the training and validation loss values for all investigated conditions, showing no sign of overfitting. The solid line represents the mean training loss function over all ten folds. The dashed lines depict the mean validation loss over all ten folds. Table I: Results of all trained models over all ten folds. The value in parentheses represents the standard deviation. MAE is calculated on the output values in the original scale. Norm. MAE is calculated on the normalized output of the network; the lower, the better. | Nb. Param. | MAE |  | Norm.
---
MAE
MAE |  | Norm.
---
MAE
Over all Channels | Electrodes Only
1 | 806K | 18.977 | 0.228 | 17.147 | 0.232
(1.719) | (0.022) | (1.642) | (0.023)
2 | 805K | 23.220 | 0.237 | 18.941 | 0.239
(1.403) | (0.018) | (1.458) | (0.020)
3 | 804K | 25.739 | 0.245 | 19.012 | 0.240
(1.539) | (0.019) | (1.498) | (0.021)
4 | 807K | 21.944 | 0.233 | 19.003 | 0.240
(1.514) | (0.019) | (1.473) | (0.020)
5 | 819K | 21.809 | 0.234 | 19.194 | 0.242
(1.548) | (0.020) | (1.529) | (0.021)
6 | 812K | 22.319 | 0.234 | 19.038 | 0.240
(1.469) | (0.019) | (1.468) | (0.021)

Figure 4 shows no sign of overfitting when using shorter intervals for the
force values, as both training and validation loss curves converge over time.
Table I reveals no apparent decrease in the metrics between all six input
vector combinations when testing on a fixed temperature value equal to the
mean temperature of the entire dataset. However, relying solely on the mean
MAE over ten folds for comparison may be misleading. A more accurate approach
is to compute paired differences of normalized MAE for each fold separately
and average these differences. This ensures a fair comparison, as the same ten
folds are consistently used across all combination experiments.

We calculate the paired differences for the input combination used by the
authors against all other investigated combinations. Statistical assessment is
conducted using a left-tailed paired $t$-test. Since cross-validation with a
single dataset involves reusing different training data points in multiple
folds, it violates the independence assumptions of the paired $t$-test [32].
Hence, we employ the corrected paired $t$-test by Nadeau and Bengio [33],
which accounts for the number of training and test data points and addresses
the issue of high type I errors in the common paired $t$-test. The null
hypothesis assumes that the paired difference is equal to $0$, indicating no
differences. The alternative hypothesis suggests that the normalized MAE of
the first distribution is lower than that of the second distribution. We
report the paired normalized MAE difference and the $p$-value in Table II.

Table II: Significance Test with the corrected paired $t$-test [33] conducted on different input combinations pairs. The first value depicts the paired normalized MAE difference in percent over the ten folds, and the second value between parenthesis represents the $p$-value. | 1 vs 2 | 1 vs 3 | 1 vs 4 | 1 vs 5 | 1 vs 6
---|---|---|---|---|---
Paired MAE Diff. | -0.893% | -1.763% | -0.513% | -0.630% | -0.655%
$p$-value | (0.004) | (0.000) | (0.044) | (0.009) | (0.013)

Table II indicates that all comparisons are significant with $p<0.05$. These
results show that adding force values from past and future timesteps, 1 vs 2
and 1 vs 3, decreases the error by ca. $0.9\%$ and $1.8\%$, respectively. In
addition, using shorter intervals of $10\;ms$, 1 vs 5, or including force
values at timesteps $T-5$ and $T+5$, 1 vs 4, does not reduce the error values.
We will still consider these investigations and analyses and check different
input combinations for our proposed solutions.

##  III Methodology

Based on the results presented in the previous section, we identify two areas
of improvement for the BioTac simulation. The first concerns using or omitting
temperature as an input value. The second concerns the windowing size used for
the force and position values. This section describes the methodology used to
optimize the BioTac simulation.

###  III-A Training without Temperature Readings

Firstly, we train the baseline network $B$ provided by Ruppel et al.[10] and
use the mean temperature value of the dataset as input, and we use these
scores as baseline. Next, we implement three approaches: a classical method
using XGBoost [34], a gradient-boosting algorithm. XGBoost is fast and
suitable for predicting continuous output variables [35]. A feed-forward deep
neural network, given that this type of network is also commonly used in
regression tasks [36]. Finally, we implement a transformer network since these
proved beneficial in several time-series tasks [37]. Inspired by vision
transformers (ViT) [38], we reformulate the classification task used in the
ViT into a regression task by using a transformer encoder to predict the
output vector. We present in Figure 5 the used backbone architecture for the
transformer network [38].

Figure 5: Visualization of our used transformer architecture, based on
Dosovitskiy et al. [38]. L stands for the number of attention blocks used in
the transformer encoder.

We use Ruppel et al. [10]’s normalization, input, and output vectors with
minor changes. For instance, we use the same input vector but without the
temperature value as input resulting in an input vector of 12 values, i.e.,
$x$, $y$, and $z$ at timestep $T$, force values ($F_{xt}$, $F_{yt}$, and
$F_{zt}$) at timesteps $T-10$, $T$, and $T+10$. The output value used consists
of all 19 electrodes signals $e_{1}$…$e_{19}$, the pressure values $pdc$, and
the dynamic vibration $pac_{0}$ at timestep $T$, i.e., 21 output values. All
channels’ input and output values are normalized using the z-score
normalization.

In contrast to Ruppel et al., [10], we thoroughly fine-tune the
hyperparameters of the tested models, namely XGBoost, feed-forward deep neural
network, and transformer, using SMAC3 [39]. The normalized mean absolute error
(MAE) is used as a cost function. We randomly split the data into training and
validation sets (90% / 10%) with random seed values. To select the best
configuration, we utilize hyperband [40] as an intensifier for the feed-
forward network and the transformer, and we use a maximum budget of 30 epochs.
The hyperparameters and their corresponding search space are reported in the
supplementary material.

XGBoost does not support multi-variable output vectors. Thus, we fit one
regressor for each channel. All regressors use the same hyperparameters
determined with SMAC3. During training, we utilize the MAE loss function. The
feed-forward network and the transformer are optimized using Adam [41], and we
use the loss function presented in Equation 1.

Once the hyperparameters are determined, we use them to train our models from
scratch. We split the dataset into training, validation, and test set in 80% /
10% / 10%, respectively. We perform early stopping on the validation loss to
avoid overfitting. In contrast, k-fold cross-validation with ten folds is
executed for a robust test set evaluation.

###  III-B Varying the Windowing Size

We also thoroughly investigate the choice of the input window for all tested
approaches while omitting the temperature values. Beyond adjusting the force
values sampling interval proposed in Section II-B, we explore two additional
input combinations. For both input combinations 1 and 5, we include previous
and next positional values within the sampling window to provide the network
with contextual touch location information. These are respectively denoted as
input combinations 7 and 8:

  1. 7.

Current, last, and next position and forces $F_{t}$: $\forall\text{
}t\in\\{T,T-10,T+10\\}$.

  2. 8.

Current, last, and next position and forces $F_{t}$: $\forall\text{
}t\in[T-10,T+10],\quad t\in\mathbb{N}$.

We only test within the timestep range of $t-10$ to $t+10$ since the force
values’ sampling frequency of 100 Hz is a multiple of 10, which aligns with
the position’s sampling frequency of 10 Hz. We use the SMAC procedure
described in the previous Subsection III-A with the same hyperparameter search
space. More details are in the supplementary material.

##  IV Results

Here, we first report the results of our tested approaches that are trained
without temperature. Next, we delve into a detailed investigation of a better
input window.

###  IV-A Training without Temperature Readings

After executing SMAC for all approaches, the models are trained, and the mean
absolute error is calculated. We also determined the number of learnable
parameters. For XGBoost regressors, we calculate the number of parameters by
counting the number of nodes in all regressors across channels after training
them. We then average this count over the folds, accounting for variability
across each fold. The results are outlined in Table III.

Table III: Results of the all trained models over all ten folds. The value in parentheses represents the standard deviation. The metrics are calculated over all channels and over all electrodes. MAE is calculated on the output values in the original scale. Norm. MAE is calculated on the normalized output of the network; the lower, the better. | Nb. | MAE | Norm. | MAE | Norm.
---|---|---|---|---|---
| Param. | MAE | MAE
|  | Over all Channels | Electrodes Only
Ruppel et al. [10] | 806K | 18.977 | 0.228 | 17.147 | 0.232
(1.719) | (0.022) | (1.642) | (0.023)
Our XGBoost | 1584K | 13.368 | 0.150 | 11.446 | 0.150
Regressor | (1.340) | (0.014) | (1.204) | (0.015)
Our Feed-Forward | 2233K | 14.693 | 0.168 | 12.724 | 0.169
Neural Network | (1.386) | (0.015) | (1.252) | (0.016)
Our Transformer | 599K | 13.760 | 0.156 | 11.736 | 0.155
Encoder | (1.400) | (0.016) | (1.280) | (0.016)

From Table III, it can be seen that the XGBoost and transformer approaches
without temperature readings achieved a statistically significant improvement
compared to the baseline. The absolute improvement measured in the normalized
MAE ranges from $6.0\%$ to $7.8\%$. The transformer encoder stands out for its
compactness and reduced number of trainable parameters. While the XGBoost
reaches the lowest normalized mean absolute error with $15.0\%$.

We also present the normalized MAE for each channel in Figure 6. There, it can
be seen that the prediction errors differ substantially between channels.
Particularly, electrodes 7, 8, 9, and 10, located at the tip of the BioTac
sensor, result in the highest errors. This phenomenon is also seen for the
XGBoost regressors, which have a dedicated regressor for each channel. This
suggests a high non-linear dynamic of those electrodes.

Figure 6: Distribution of the normalized MAE for each channel and all models.

Another possible reason for the differences in prediction errors could be the
different number of samples available for each channel; i.e., the dataset
might be unbalanced. We assessed that by selecting all data points at the
beginning of a contact cycle having a force measurement higher than $0.3\;N$,
located close to the BioTac sensor surface within a distance of less than
$2\;mm$. Subsequently, we identify the nearest electrode to each contact point
using the 3D spatial positions of each electrode provided by Lin et al. [31].
The results are depicted in Figure 7.

Figure 7: Distribution of nearest contact points for each electrode.

Figure 7 highlights a non-uniform distribution of the number of touches near
each electrode, implying an unbalanced dataset. However, this may not be the
only reason contributing to the high MAE for specific electrodes. Based on the
work of Lin et al. [31], we know that the BioTac sensor does not have radial
symmetry. In addition, the fluid volume is not the same throughout the sensor.
Particularly near electrodes 7, 8, 9, and 10, there is a higher fluid volume,
which is partially visible in Figure 1. This supports the hypothesis that some
of the error is due to the non-linear dynamics of the sensor.

###  IV-B Varying the Windowing Size

After executing SMAC for all approaches, the models are trained, and the mean
absolute error is calculated. The best hyperparameters for each trained model
and input combination and the result metrics values are reported in the
supplementary material.

Figure 8 presents the distribution of the normalized MAE for all combinations
across all approaches. XGBoost attains the lowest MAE when trained with input
combination 8, incorporating force values with shorter intervals and the last
and next position, reaching a value of 14.8%, which represents an $8.0\%$
improvement over the baseline. Likewise, our transformer encoder achieves a
$7.8\%$ improvement over the baseline with the same input combination.

Figure 8: The Normalized MAE was calculated over all channels for the four
models and all input combinations.

We also conduct a corrected paired significance $t$-test [33] and report the
paired normalized MAE difference and $p$-value for three specific comparison
combinations: 5 vs 1, 7 vs 1 and 8 vs 1. The null hypothesis assumes that the
paired difference is equal to $0$, indicating no differences. The alternative
hypothesis suggests that the normalized MAE of the first distribution is lower
than that of the second distribution. Comparison 5 vs 1 assesses the impact of
using force values of shorter intervals on error reduction, 7 vs 1 examines
the effect of including the last and next different position, and 8 vs 1
evaluates the significance of incorporating more force values at shorter
intervals along with the last and next position in error reduction. The paired
normalized MAE difference, $t$-statistic and $p$-values for various
combination pairs across all models are reported in the supplementary
material.

Table IV indicates that adding more force values of shorter intervals does not
directly lead to a reduction in error (5 vs 1), and the same applies when
adding the last and next contact position (7 vs 1). However, combining both
information yielded statistically significant improvements for the XGBoost
regressor and the transformer encoder in the 8 vs 1 comparison. However, they
are small, ranging between 0.2% and 0.6%.

Table IV: Significance Test with the corrected paired $t$-test [33] conducted on different input combinations. The first value depicts the paired normalized MAE difference in percent over the ten folds, and the second value between parenthesis represents the $p$-value. | 5 vs 1 | 7 vs 1 | 8 vs 1
---|---|---|---
Our XGBoost | 0.282% | -0.055% | -0.208%
Regressor | (0.999) | (0.190) | (0.009)
Our Feed-Forward | 0.217% | 0.862% | 0.402%
Neural Network | (0.715) | (0.961) | (0.758)
Our Transformer | 0.077% | 0.188% | -0.662%
Encoder | (0.631) | (0.672) | (0.047)
Table V: Significance Test with the corrected paired $t$-test [33] conducted for all models. The first value depicts the paired normalized MAE difference value in percent, and the second value between parenthesis represents the $p$-value. vs | Our XGBoost | Our XGBoost | Our Transformer
---|---|---|---
Our FFNN | Our Transformer | Our FFNN
1 | -1.778% | -0.531% | -1.247%
(0.000) | (0.074) | (0.002)
2 | -1.360% | -0.259% | -1.101%
(0.001) | (0.250) | (0.000)
3 | -1.272% | 0.851% | -2.123%
(0.003) | (0.996) | (0.000)
4 | -1.989% | -0.032% | -1.956%
(0.000) | (0.459) | (0.000)
5 | -1.279% | -0.173% | -1.107%
(0.002) | (0.322) | (0.002)
6 | -1.106% | 0.092% | -1.197%
(0.003) | (0.616) | (0.002)
7 | -2.696% | -0.774% | -1.921%
(0.000) | (0.105) | (0.000)
8 | -2.389% | -0.077% | -2.312%
(0.001) | (0.439) | (0.002)

We also analyzed the significance between the tested approaches for all input
combinations and reported the results. Table V reveals that both the XGBoost
regressor and the transformer encoder outperform the feed-forward neural
network across all input combinations, demonstrating an error reduction
ranging between $1.1\%$ and $2.7\%$. However, a direct comparison between the
XGBoost and transformer yields no significant values, suggesting no difference
between the paired MAE.

Considering the importance of inference time during the deployment phase, we
compute the inference time over 100 random inputs for all approaches and
subsequently average the results. We run the inference calculation on a system
with an AMD Ryzen 7 pro 4750u CPU. Figure 9 shows the inference time against
the number of learnable parameters for all approaches and input combinations.
For Ruppel et al. [10] baseline network, we only plot results for the first
input combination, as all other combinations do not differ meaningfully in
inference time or numbers of parameters.

Figure 9: Visualization of inference time against the number of parameters for
each model over all input combinations.

Figure 9 reveals that the transformer network has the lowest number of
parameters compared to all other approaches, with an inference time ranging
from $1.0\;ms$ up to $3.0\;ms$, i.e., at least two to five times slower than
the baseline. This is attributed to the quadratic time complexity of the
attention blocks [42]. XGBoost is the fastest of the tested methods, with
approximately half the inference time of the baseline network.

Considering that the number of floating-point operations per second (FLOPS)
offers additional valuable insights into the energy consumption and inference
efficiency of neural networks [43], we include the FLOPS count for all the
neural networks, i.e., the baseline network, our feed-forward network, and our
transformer encoder, in our supplementary material.

##  V Conclusion

We present optimizations to simulate tactile signals from the BioTac 2P
sensor. Our contributions include a thorough analysis of Ruppel et al.’s work
[10] and a thorough analysis of three alternative solutions. We focused on two
areas: the use of temperature readings in the training pipeline of the
network, and the choice of other variables in the input vector and the
windowing.

We introduced three alternative models, i.e., an XGBoost regressor, a feed-
forward neural network, and an adapted transformer encoder. Our results
demonstrate that XGBoost and the transformer encoder can achieve significantly
lower error values than the baseline. The absolute error reduction stands at
$8.0\%$, while the relative error shows a statistically significant
improvement of $14.9\%$. Furthermore, our investigations reveal a
statistically significant improvement when incorporating force values from
future timesteps and previous and next position values. Albeit, there are
small between 0.2% up to 0.6% depending on the model used.

Notably, the XGBoost regressor has the lowest inference time, making it
preferable for scenarios where simulation inference time is critical. Despite
having less trainable parameters than all other models, the transformers
exhibit the highest inference time of the tested models.

We also report on the limitations of the dataset. Particularly, the dataset is
unbalanced and only includes a single indenter type. However, some of the
errors can be attributed to the non-linear dynamics of the sensor, probably
caused by the non-radial symmetry and non-uniform fluid volume throughout the
sensor.

As future work, the dataset needs to be improved and extended to include
different BioTac 2P sensors, varied surrounding temperatures, and several
indenter shapes to enhance model robustness and generalizability, mitigating
the absence of temperature information in the simulation environment and other
non-linear dynamics.

Due to the differentiated performance of different electrodes due to non-
linearities and unbalanced dataset, we also suggest training an ensemble of
transformer networks to better deal with the non-linear dynamic of the sensor.

## References

  * [1] N. Navarro-Guerrero, S. Toprak, J. Josifovski, and L. Jamone, “Visuo-Haptic Object Perception for Robots: An Overview,” _Autonomous Robots_ , vol. 47, no. 4, pp. 377–403, 2023.
  * [2] A. Parmiggiani, S. Dussoni, L. Natale, and G. Metta, “Tactile Sensors,” in _Encyclopedia of Robotics_. Springer, 2020, pp. 1–11.
  * [3] J. Josifovski, M. Malmir, N. Klarmann, B. L. Zagar, N. Navarro-Guerrero, and A. Knoll, “Analysis of Randomization Effects on Sim2Real Transfer in Reinforcement Learning for Robotic Manipulation Tasks,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022, pp. 10 193–10 200.
  * [4] S. Auddy, J. Hollenstein, M. Saveriano, A. Rodríguez-Sánchez, and J. Piater, “Continual Learning from Demonstration of Robotics Skills,” _Robotics and Autonomous Systems_ , vol. 165, p. 104427, 2023.
  * [5] S. Wang, M. Lambeta, P.-W. Chou, and R. Calandra, “TACTO: A Fast, Flexible, and Open-Source Simulator for High-Resolution Vision-Based Tactile Sensors,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 3930–3937, 2022.
  * [6] M. Lambeta, P.-W. Chou, S. Tian, B. Yang, B. Maloon, V. R. Most, D. Stroud, R. Santos, A. Byagowi, G. Kammerer, D. Jayaraman, and R. Calandra, “DIGIT: A Novel Design for a Low-Cost Compact High-Resolution Tactile Sensor With Application to In-Hand Manipulation,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 3, pp. 3838–3845, 2020.
  * [7] A. Padmanabha, F. Ebert, S. Tian, R. Calandra, C. Finn, and S. Levine, “OmniTact: A Multi-Directional High-Resolution Touch Sensor,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 618–624.
  * [8] A. Melnik, L. Lach, M. Plappert, T. Korthals, R. Haschke, and H. Ritter, “Using Tactile Sensing to Improve the Sample Efficiency and Performance of Deep Deterministic Policy Gradients for Simulated in-Hand Manipulation Tasks,” _Frontiers in Robotics and AI_ , vol. 8, no. 538773, 2021.
  * [9] J, Geukes, M. Nakatenus and R. Calandra. “Gazebo plugin for the icub skin,” 2017. Online Available: https://github.com/robertocalandra/icub-gazebo-skin.
  * [10] P. Ruppel, Y. Jonetzko, M. Görner, N. Hendrich, and J. Zhang, “Simulation of the SynTouch BioTac Sensor,” in _International Conference on Intelligent Autonomous Systems (IAS)_ , ser. AISC, vol. 867. Springer International Publishing, 2019, pp. 374–387.
  * [11] Y. S. Narang, B. Sundaralingam, K. Van Wyk, A. Mousavian, and D. Fox, “Interpreting and Predicting Tactile Signals for the SynTouch BioTac,” _The International Journal of Robotics Research_ , vol. 40, no. 12-14, pp. 1467–1487, 2021.
  * [12] B. S. Zapata-Impata, P. Gil, Y. Mezouar, and F. Torres, “Generation of Tactile Data From 3D Vision and Target Robotic Grasps,” _IEEE Transactions on Haptics_ , vol. 14, no. 1, pp. 57–67, 2021.
  * [13] N. Wettels, D. Popovic, V. J. Santos, R. S. Johansson, and G. E. Loeb, “Biomimetic Tactile Sensor for Control of Grip,” in _IEEE International Conference on Rehabilitation Robotics (ICORR)_ , 2007, pp. 923–932.
  * [14] N. Wettels, J. A. Fishel, and G. E. Loeb, “Multimodal Tactile Sensor,” in _The Human Hand as an Inspiration for Robot Hand Development_ , ser. Springer Tracts in Advanced Robotics. Springer International Publishing, 2014, no. 95, pp. 405–429.
  * [15] L. Natale and G. Cannata, “Tactile Sensing,” in _Humanoid Robotics: A Reference_. Springer Netherlands, 2019, pp. 2539–2561.
  * [16] Y. Chebotar, K. Hausman, Z. Su, G. Sukhatme, and S. Schaal, “Self-Supervised Regrasping Using Spatio-Temporal Tactile Features and Reinforcement Learning,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2016, pp. 1960–1966.
  * [17] D. Xu, G. E. Loeb, and J. A. Fishel, “Tactile Identification of Objects Using Bayesian Exploration,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2013, pp. 3056–3061.
  * [18] M. S. Arian, C. A. Blaine, G. E. Loeb, and J. A. Fishel, “Using the BioTac as a Tumor Localization Tool,” in _IEEE Haptics Symposium (HAPTICS)_ , 2014, pp. 443–448.
  * [19] C. Pacchierotti, D. Prattichizzo, and K. J. Kuchenbecker, “Cutaneous Feedback of Fingertip Deformation and Vibration for Palpation in Robotic Surgery,” _IEEE Transactions on Biomedical Engineering_ , vol. 63, no. 2, pp. 278–287, 2016.
  * [20] SynTouch, BioTac® Product Manual, August, 2018, V21: https://www.syntouchinc.com/wp-content/uploads/2018/08/BioTac-Manual-V.21.pdf.
  * [21] The Shadow Dexterous Hand: https://www.shadowrobot.com/dexterous-hand-series/.
  * [22] ATI Industrial Automation Inc.: F/T Sensor Nano17 https://www.ati-ia.com/products/ft/ft_models.aspx?id=Nano17.
  * [23] N. Koenig and A. Howard, “Design and Use Paradigms for Gazebo, an Open-Source Multi-Robot Simulator,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , vol. 3, 2004, pp. 2149–2154.
  * [24] C. R. Qi, L. Yi, H. Su, and L. J. Guibas, “PointNet++: Deep Hierarchical Feature Learning on Point Sets in a Metric Space,” in _International Conference on Neural Information Processing Systems (NIPS)_ , ser. NIPS’17. Curran Associates Inc., 2017, pp. 5105–5114.
  * [25] NVIDIA Isaac Gym Simulator: https://developer.nvidia.com/isaac-gym.
  * [26] Y. Narang, B. Sundaralingam, M. Macklin, A. Mousavian, and D. Fox, “Sim-to-Real for Robotic Tactile Sensing Via Physics-Based Simulation and Learned Latent Projections,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 6444–6451.
  * [27] T. Taunyazov, L. S. Song, E. Lim, H. H. See, D. Lee, B. C. Tee, and H. Soh, “Extended Tactile Perception: Vibration Sensing through Tools and Grasped Objects,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2021, pp. 1755–1762.
  * [28] SynTouch, BioTac® SP Product Manual, August, 2018, V3: https://www.syntouchinc.com/wp-content/uploads/2018/08/BioTac-SP-Manual-V.3.pdf.
  * [29] G. Olmschenk, Z. Zhu, and H. Tang, “Generalizing Semi-Supervised Generative Adversarial Networks to Regression Using Feature Contrasting,” _Computer Vision and Image Understanding_ , vol. 186, pp. 1–12, 2019.
  * [30] E. Olson, “AprilTag: A Robust and Flexible Visual Fiducial System,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2011, pp. 3400–3407.
  * [31] C.-H. Lin, J. A. Fishel, and G. E. Loeb, “Estimating point of contact, force and torque in a biomimetic tactile sensor with deformable skin,” in _”Technical report, SynTouch LLC”_ , 2013.
  * [32] T. G. Dietterich, “Approximate Statistical Tests for Comparing Supervised Classification Learning Algorithms,” _Neural Computation_ , vol. 10, no. 7, pp. 1895–1923, Oct. 1998.
  * [33] C. Nadeau and Y. Bengio, “Inference for the Generalization Error,” _Machine Learning_ , vol. 52, no. 3, pp. 239–281, Sep. 2003.
  * [34] T. Chen and C. Guestrin, “XGBoost: A Scalable Tree Boosting System,” in _ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD)_ , ser. KDD ’16. ACM, 2016, pp. 785–794.
  * [35] C. Bentéjac, A. Csörgő, and G. Martínez-Muñoz, “A Comparative Analysis of Gradient Boosting Algorithms,” _Artificial Intelligence Review_ , vol. 54, no. 3, pp. 1937–1967, 2021.
  * [36] L. Das, A. Sivaram, and V. Venkatasubramanian, “Hidden Representations in Deep Neural Networks: Part 2. Regression Problems,” _Computers & Chemical Engineering_, vol. 139, p. 106895, 2020.
  * [37] Q. Wen, T. Zhou, C. Zhang, W. Chen, Z. Ma, J. Yan, and L. Sun, “Transformers in Time Series: A Survey,” in _International Joint Conference on Artificial Intelligence (IJCAI)_ , vol. 6, 2023, pp. 6778–6786.
  * [38] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, J. Uszkoreit, and N. Houlsby, “An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale,” in _International Conference on Learning Representations (ICLR)_ , 2021.
  * [39] M. Lindauer, K. Eggensperger, M. Feurer, A. Biedenkapp, D. Deng, C. Benjamins, T. Ruhkopf, R. Sass, and F. Hutter, “SMAC3: A Versatile Bayesian Optimization Package for Hyperparameter Optimization,” _Journal of Machine Learning Research_ , vol. 23, no. 54, pp. 1–9, 2022.
  * [40] L. Li, K. Jamieson, G. DeSalvo, A. Rostamizadeh, and A. Talwalkar, “Hyperband: A Novel Bandit-Based Approach to Hyperparameter Optimization,” _Journal of Machine Learning Research_ , vol. 18, no. 185, pp. 1–52, 2018.
  * [41] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” in _International Conference on Learning Representations (ICLR)_ , ser. 3rd, 2015, p. 15.
  * [42] Y. Tay, M. Dehghani, D. Bahri, and D. Metzler, “Efficient Transformers: A Survey,” _ACM Computing Surveys_ , vol. 55, no. 6, pp. 109:1–109:28, 2022.
  * [43] R. Tang, W. Wang, Z. Tu, and J. Lin, “An Experimental Analysis of the Power Consumption of Convolutional Neural Networks for Keyword Spotting,” in _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , Calgary, AB, Canada, 2018, pp. 5479–5483.

Table VI: Summary of the hyperparameter search space used for the XGBoost Regressor, when running SMAC3 [39] and the selected configuration for all input combinations. Hyperparameters | Our XGBoost Regressor
---|---
Search Space | Input 1 | Input 2 | Input 3 | Input 4 | Input 5 | Input 6 | Input 7 | Input 8
eta |  | Uniform Float
---
$\in$ [0.0001, 0.5]
0.0431 | 0.06818 | 0.04942 | 0.04457 | 0.07178 | 0.06647 | 0.04901 | 0.04798
gamma |  | Uniform Integer
---
$\in$ [0, 10]
1 | 7 | 2 | 1 | 7 | 3 | 3 | 2
number estimate |  | Uniform Integer
---
$\in$ [100, 1000]
972 | 155 | 932 | 880 | 230 | 913 | 715 | 949
max depth |  | Uniform Integer
---
$\in$ [1, 10]
10 | 10 | 10 | 10 | 10 | 10 | 10 | 10
min child weight |  | Uniform Integer
---
$\in$ [1, 100]
95 | 6 | 87 | 67 | 59 | 83 | 84 | 84
max delta step |  | Uniform Integer
---
$\in$ [0, 10]
7 | 10 | 0 | 1 | 5 | 3 | 10 | 1
subsample |  | Uniform Float
---
$\in$ [0.5, 1]
0.647 | 0.9632 | 0.5042 | 0.5068 | 0.6114 | 0.5486 | 0.507 | 0.5038
colsample bytree |  | Uniform Float
---
$\in$ [0.5, 1]
0.9825 | 0.9517 | 0.8671 | 0.868 | 0.7449 | 0.8242 | 0.9927 | 0.7814
colsample bylevel |  | Uniform Float
---
$\in$ [0.5, 1]
0.9819 | 0.9223 | 0.8572 | 0.839 | 0.9264 | 0.7642 | 0.9479 | 0.8407
colsample bynode |  | Uniform Float
---
$\in$ [0.5, 1]
0.8042 | 0.9155 | 0.6173 | 0.7416 | 0.9802 | 0.924 | 0.8621 | 0.9989
Table VII: Table summarizing the hyperparameter search space used for the feed-forward deep neural network, when running SMAC3 [39] and the selected configuration for all input combinations. Hyperparameters | Our Feed-Forward Neural Network
---|---
Search Space | Input 1 | Input 2 | Input 3 | Input 4 | Input 5 | Input 6 | Input 7 | Input 8
Batch Size |  | Categorical
---
[256, 512]
256 | 512 | 256 | 512 | 512 | 512 | 512 | 512
Learning Rate |  | Categorical
---
$a\times e^{-c}$
for $a\in\mathbb{N}^{+}$ and $\in{[}1,9{]}$
$c\in\mathbb{N}^{+}$ and $\in{[}2,5{]}$
0.0003 | 0.0002 | 0.0005 | 0.0006 | 0.0003 | 0.0006 | 0.0004 | 0.0005
| Number
---
of Layers ($L$)
| Uniform Int
---
Lower: 4
Upper: 12
7 | 10 | 7 | 7 | 7 | 9 | 7 | 8
| Numebr of Neurons
---
in Layer i
for $i\in{[}0,L{]}$
| Uniform Int
---
Lower: 50
Upper: 1000
Step: 10
| [860, 670,
---
160, 580,
900, 1000,
440]
| [590,300,
---
820,520,
90,670,
850,120,
330,570]
| [620, 470,
---
120, 620,
830, 890,
350]
| [730, 390,
---
120, 630,
770, 720,
380]
| [820,740,
---
190,740,
1000,850,
430]
| [580,710,
---
600,170,
270,350,
70,780,
690]
| [580,580,
---
160,450,
920,920,
380]
| [470,600,
---
440,830,
790,900,
190,270]
| Activation function
---
in Layer i
for $i\in{[}0,L{]}$
| Categorical
---
[sigmoid, relu,
hardtanh, tanh,
leakyrelu, elu]
| [hardtanh, tanh,
---
elu, relu,
leakyrelu, relu,
leakyrelu]
| [elu, relu,
---
hardtanh, elu,
leakyrelu, relu,
elu, elu,
elu, leakyrelu]
| [tanh, elu,
---
elu,
leakyrelu,
leakyrelu,
elu,
leakyrelu]
| [tanh, tanh,
---
elu,
hardtanh,
sigmoid,
leakyrelu,
leakyrelu]
| [tanh, relu,
---
elu, leakyrelu,
elu, relu,
leakyrelu]
| [tanh, relu,
---
leakyrelu, elu,
leakyrelu,
leakyrelu,
leakyrelu,
hardtanh,
leakyrelu]
| [tanh, relu,
---
leakyrelu,
hardtanh,
sigmoid, relu,
leakyrelu]
| [tanh, relu,
---
elu, leakyrelu,
hardtanh,
leakyrelu,
elu, leakyrelu]
| Negative Slope
---
for Leakyrelu
| Categorical
---
$a\times e^{-1}$
for $a\in\mathbb{N}^{+}$ and $\in{[}1,9{]}$
0.5 | 0.7 | 0.4 | 0.2 | 0.1 | 0.4 | 0.3 | 0.6
Table VIII: Table summarizing the hyperparameter search space used for the transformer encoder, when running SMAC3 [39] and the selected configuration for all input combinations. Hyperparameters | Our Transformer Encoder
---|---
Search Space | Input 1 | Input 2 | Input 3 | Input 4 | Input 5 | Input 6 | Input 7 | Input 8
Batch Size |  | Categorical
---
[256, 512]
512 | 512 | 512 | 256 | 512 | 256 | 512 | 256
Learning Rate |  | Categorical
---
$a\times e^{-c}$
for $a\in\mathbb{N}^{+}$ and $\in{[}1,9{]}$
$c\in\mathbb{N}^{+}$ and $\in{[}2,5{]}$
0.00009 | 0.00004 | 0.00009 | 0.0004 | 0.004 | 0.0003 | 0.00004 | 0.002
| Number
---
of Layers ($L$)
| Uniform Int
---
Lower: 2
Upper: 8
3 | 8 | 3 | 2 | 7 | 4 | 7 | 2
| Number of
---
Multi-Heads
| Categorical
---
[1, 2, 4, 8]
4 | 8 | 4 | 8 | 1 | 4 | 1 | 8
Dropout Rate |  | Categorical
---
[0.0, 0.1, 0.2, 0.3, 0.4, 0.5]
0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
| Embedding
---
Dimensions
| Categorical
---
[32, 64, 128, 256, 512, 1024]
128 | 128 | 128 | 128 | 64 | 32 | 128 | 64
| Hidden
---
Dimension
| Categorical
---
[64, 128, 256, 512, 1024]
512 | 512 | 512 | 512 | 128 | 256 | 128 | 256
Table IX: Results of all trained networks over all 10 folds for all input
combinations. The value in parentheses represents the standard deviation. The
metrics are calculated over all channels and over all electrodes. MAE is
calculated on the output values in the original scale. Norm. MAE is calculated
on the normalized output of the network; the lower, the better.

|  Nb. Param. | MSE |  | Norm.
---
MSE
MSE |  | Norm.
---
MSE
Over all Channels | Electrodes Only
Ruppel et al.[10] | 1 | 806K | 18.977 | 0.228 | 17.147 | 0.232
(1.719) | (0.022) | (1.642) | (0.023)
2 | 805K | 23.220 | 0.237 | 18.941 | 0.239
(1.403) | (0.018) | (1.458) | (0.020)
3 | 804K | 25.739 | 0.245 | 19.012 | 0.240
(1.539) | (0.019) | (1.498) | (0.021)
4 | 807K | 21.944 | 0.233 | 19.003 | 0.240
(1.514) | (0.019) | (1.473) | (0.020)
5 | 819K | 21.809 | 0.234 | 19.194 | 0.242
(1.548) | (0.020) | (1.529) | (0.021)
6 | 812K | 22.319 | 0.234 | 19.038 | 0.240
(1.469) | (0.019) | (1.468) | (0.021)
7 | 809K | 22.497 | 0.234 | 18.811 | 0.238
(1.466) | (0.019) | (1.449) | (0.020)
8 | 822K | 21.555 | 0.231 | 18.927 | 0.239
(1.504) | (0.019) | (1.489) | (0.020)
Our XGBoost Regressor | 1 | 1584K | 13.368 | 0.150 | 11.446 | 0.150
(1.340) | (0.015) | (1.204) | (0.015)
2 | 1041K | 14.372 | 0.159 | 12.002 | 0.158
(1.446) | (0.016) | (1.298) | (0.016)
3 | 1209K | 15.994 | 0.168 | 12.373 | 0.163
(1.622) | (0.017) | (1.340) | (0.017)
4 | 1694K | 13.275 | 0.152 | 11.638 | 0.153
(1.430) | (0.016) | (1.335) | (0.016)
5 | 840K | 13.407 | 0.153 | 11.774 | 0.154
(1.393) | (0.015) | (1.277) | (0.016)
6 | 1046K | 13.785 | 0.156 | 11.866 | 0.156
(1.385) | (0.015) | (1.277) | (0.016)
7 | 1239K | 13.299 | 0.150 | 11.375 | 0.150
(1.330) | (0.015) | (1.182) | (0.015)
8 | 1331K | 12.873 | 0.148 | 11.385 | 0.150
(1.255) | (0.014) | (1.156) | (0.014)

|  Nb. Param. | MSE |  | Norm.
---
MSE
MSE |  | Norm.
---
MSE
Over all Channels | Electrodes Only
Our Feed-Forward Neural Network | 1 | 2233K | 14.693 | 0.168 | 12.724 | 0.169
(1.386) | (0.015) | (1.252) | (0.016)
2 | 2881K | 15.624 | 0.173 | 12.965 | 0.171
(1.536) | (0.017) | (1.388) | (0.018)
3 | 1701K | 17.010 | 0.181 | 13.265 | 0.175
(1.314) | (0.014) | (1.105) | (0.014)
4 | 1478K | 14.712 | 0.171 | 13.109 | 0.173
(1.317) | (0.015) | (1.247) | (0.016)
5 | 2554K | 14.056 | 0.166 | 12.721 | 0.169
(1.213) | (0.013) | (1.100) | (0.014)
6 | 1124K | 14.476 | 0.167 | 12.676 | 0.168
(1.150) | (0.013) | (1.043) | (0.013)
7 | 1794K | 15.231 | 0.177 | 13.351 | 0.178
(1.344) | (0.016) | (1.242) | (0.016)
8 | 2490K | 14.517 | 0.172 | 13.204 | 0.175
(1.450) | (0.017) | (1.404) | (0.017)
Our Transformer Encoder | 1 | 599K | 13.760 | 0.156 | 11.736 | 0.155
(1.400) | (0.016) | (1.280) | (0.016)
2 | 203K | 14.534 | 0.162 | 12.102 | 0.160
(1.482) | (0.017) | (1.374) | (0.017)
3 | 598K | 15.229 | 0.160 | 11.595 | 0.154
(1.617) | (0.017) | (1.328) | (0.017)
4 | 401K | 13.204 | 0.152 | 11.564 | 0.153
(1.363) | (0.016) | (1.266) | (0.016)
5 | 237K | 13.441 | 0.155 | 11.824 | 0.156
(1.241) | (0.014) | (1.147) | (0.015)
6 | 114K | 13.669 | 0.155 | 11.719 | 0.155
(1.505) | (0.018) | (1.376) | (0.018)
7 | 701K | 13.927 | 0.158 | 11.838 | 0.157
(1.614) | (0.019) | (1.468) | (0.019)
8 | 103K | 12.984 | 0.149 | 11.334 | 0.150
(1.552) | (0.018) | (1.473) | (0.019)

Table X: Extended Significance Test with the corrected paired $t$-test [33] conducted on different input combinations for all networks. The first value depicts the paired normalized MAE difference in percent over the ten folds, the second value represents $t$-statistic, and the third value between parenthesis represents the $p$-value. | 1 vs 2 | 1 vs 3 | 1 vs 4 | 1 vs 5 | 1 vs 6 | 7 vs 1 | 8 vs 1 | 5 vs 6 | 8 vs 5
---|---|---|---|---|---|---|---|---|---
Ruppel et al. [10] | -0.893% | -1.763% | -0.513% | -0.630% | -0.655% | 0.578% | 0.353% | -0.025% | -0.277%
(0.004) | (0.000) | (0.044) | (0.009) | (0.013) | (0.967) | (0.915) | (0.398) | (0.007)
Our XGBoost | -0.859% | -1.773% | -0.113% | -0.282% | -0.539% | -0.055% | -0.208% | -0.258% | -0.490%
Regressor | (0.000) | (0.000) | (0.129) | (0.001) | (0.000) | (0.190) | (0.009) | (0.009) | (0.001)
Our Feed-Forward | -0.441% | -1.267% | -0.324% | 0.217% | 0.133% | 0.862% | 0.402% | -0.084% | 0.620%
Neural Network | -2.222 | -5.436 | -2.184 | 0.589 | 0.373 | 1.995 | 0.729 | -0.226 | 1.296
| (0.027) | (0.000) | (0.028) | (0.715) | (0.641) | (0.961) | (0.758) | (0.413) | (0.886)
Our Transformer | -0.587% | -0.391% | 0.385% | 0.077% | 0.083% | 0.188% | -0.662% | 0.007% | -0.585%
Encoder | -2.479 | -1.400 | 1.074 | 0.344 | 0.438 | 0.461 | -1.870 | 0.019 | -1.537
| (0.018) | (0.097) | (0.845) | (0.631) | (0.664) | (0.672) | (0.047) | (0.507) | (0.079)
Table XI: Extended Significance Test with the corrected paired $t$-test [33] conducted for all network pairs. The first value depicts the paired normalized MAE difference in percent over the ten folds, the second value represents $t$-statistic, and the third value between parenthesis represents the $p$-value. vs | Our XGBoost | Our FFNN | Our Transformer | Our XGBoost | Our XGBoost | Our Transformer
---|---|---|---|---|---|---
Ruppel et al. [10] | Ruppel et al. [10] | Ruppel et al. [10] | Our FFNN | Our Transformer | Our FFNN
1 | -7.740% | -5.962% | -7.209% | -1.778% | -0.531% | -1.247%
-9.785 | -7.741 | -8.670 | -5.188 | -1.581 | -3.690
(0.000) | (0.000) | (0.000) | (0.000) | (0.074) | (0.002)
2 | -7.774% | -6.414% | -7.515% | -1.360% | -0.259% | -1.101%
-9.962 | -8.136 | -8.956 | -4.284 | -0.704 | -5.338
(0.000) | (0.000) | (0.000) | (0.001) | (0.250) | (0.000)
3 | -7.730% | -6.457% | -8.581% | -1.272% | 0.851% | -2.123%
-9.936 | -8.509 | -11.401 | -3.623 | 3.407 | -6.570
(0.000) | (0.000) | (0.000) | (0.003) | (0.996) | (0.000)
4 | -8.140% | -6.151% | -8.108% | -1.989% | -0.032% | -1.956%
-9.546 | -7.352 | -8.958 | -5.653 | -0.106 | -11.260
(0.000) | (0.000) | (0.000) | (0.000) | (0.459) | (0.000)
5 | -8.089% | -6.809% | -7.916% | -1.279% | -0.173% | -1.107%
-9.960 | -8.364 | -10.944 | -3.888 | -0.477 | -3.810
(0.000) | (0.000) | (0.000) | (0.002) | (0.322) | (0.002)
6 | -7.856% | -6.750% | -7.947% | -1.106% | 0.092% | -1.197%
-10.032 | -7.725 | -9.859 | -3.610 | 0.303 | -3.918
(0.000) | (0.000) | (0.000) | (0.003) | (0.616) | (0.002)
7 | -8.373% | -5.678% | -7.599% | -2.696% | -0.774% | -1.921%
-10.568 | -7.579 | -7.886 | -5.235 | -1.347 | -5.899
(0.000) | (0.000) | (0.000) | (0.000) | (0.105) | (0.000)
8 | -8.301% | -5.912% | -8.224% | -2.389% | -0.077% | -2.312%
-10.363 | -5.268 | -9.773 | -4.242 | -0.157 | -4.004
(0.000) | (0.000) | (0.000) | (0.001) | (0.439) | (0.002)
Table XII: Number of parameters in thousands, inference time in milliseconds, and floating-point operations per second (FLOPS) in millions for all approaches across all input combinations. The number of FLOPS is only calculated for the neural networks. | Ruppel et al.’s Network $B$ [10] | Our XGBoost Regressor | Our FF Neural Network | Our Transformer Encoder
---|---|---|---|---
| Numb. | Inference | Numb. | Numb. | Inference | Numb. | Numb. | Inference | Numb. | Numb. | Inference | Numb.
| Param. | ($ms$) | FLOPS | Param. | ($ms$) | FLOPS | Param. | ($ms$) | FLOPS | Param. | ($ms$) | FLOPS
1 | 806K | 0.628 | 1.61M | 1584K | 0.422 | - | 2233K | 0.876 | 4.46M | 599K | 1.618 | 5.97M
2 | 805K | 0.539 | 1.61M | 1041K | 0.272 | - | 2881K | 1.173 | 5.75M | 203K | 2.982 | 1.63M
3 | 804K | 0.539 | 1.61M | 1209K | 0.358 | - | 1701K | 0.741 | 3.40M | 598K | 1.396 | 3.58M
4 | 807K | 0.555 | 1.61M | 1694K | 0.432 | - | 1478K | 0.757 | 2.95M | 401K | 1.005 | 5.59M
5 | 819K | 0.572 | 1.65M | 840K | 0.304 | - | 2554K | 0.833 | 5.10M | 237K | 2.995 | 11.64M
6 | 812K | 0.552 | 1.62M | 1046K | 0.356 | - | 1124K | 0.613 | 2.24M | 114K | 1.427 | 2.98M
7 | 809K | 0.572 | 1.62M | 1239K | 0.387 | - | 1794K | 0.760 | 3.58M | 701K | 2.939 | 9.89M
8 | 822K | 0.607 | 1.65M | 1331K | 0.412 | - | 2490K | 1.002 | 4.98M | 103K | 0.988 | 5.30M