# An Ultra-Wideband Leaky Lens Antenna for Broadband Spectroscopic Imaging
Applications

Sebastian Hähnle, Ozan Yurduseven, , Sven van Berkel, , Nuria Llombart, , Juan
Bueno, Stephen J. C. Yates, Vignesh Murugesan, David J. Thoen, Andrea Neto and
Jochem J. A. Baselmans
This work is supported by the ERC COG 648135 MOSAIC. The contribution from
Nuria Llombart is supported by ERC Starting Grant ERC-2014- StG Grant LAA-THz-
CC, No. 639749.S. Hähnle, J. Bueno, S. J. C. Yates, V. Murugesan and J. J. A.
Baselmans are with the Netherlands Institute of Space Research (SRON),
Sorbonnelaan 2, Utrecht, 3584 CA Utrecht, The Netherlands (email:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>j.baselmans@sron.nl).O. Yurduseven is with Huawei Technologies Duesseldorf
GmbH, Leprosenweg 1, 82362 Weilheim in Oberbayern, Germany (e-mail:
ozan.yurduseven@ieee.orgin).S. van Berkel, N. Llombart, D. J. Thoen and A.
Neto are with the Microelectronics Department of the Electrical Engineering,
Mathemathics and Computer Science Faculty, Delft University of Technology,
Mekelweg 4, 2628 CD Delft, The Netherlands (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>a.neto@tudelft.nl).

###### Abstract

We present the design, fabrication and characterisation of a broadband leaky
lens antenna for broadband, spectroscopic imaging applications. The antenna is
designed for operation in the 300-900 GHz band. We integrate the antenna
directly into an Al-NbTiN hybrid MKID to measure the beam pattern and absolute
coupling efficiency at three frequency bands centred around 350, 650 and 850
GHz, covering the full antenna band. We find an aperture efficiency
$\eta_{ap}\approx 0.4$ over the whole frequency band, limited by lens
reflections. We find a good match with simulations for both the patterns and
efficiency, demonstrating a 1:3 bandwidth in the sub-mm wavelength range for
future on-chip spectrometers.

###### Index Terms:

Broadband antennas, Submillimeter wave devices, superconducting devices,
Ultrawideband antennas,

##  I Introduction

The emergence of on-chip spectrometers for far-infrared and sub-millimeter
astronomy, such as DESHIMA[1], SuperSpec[2] and Micro-Spec[3], emphasises the
need for efficient broadband radiation coupling in the frequency range of
$100-1000\ $GHz. These spectrometers offer the ability of measuring medium
resolution spectra ($F/\Delta F\sim 500$) with an instantaneous relative
bandwidth of up to 1:3, complementing high-resolution but bandwidth limited
heterodyne instruments. Future versions of these instruments will move towards
multi-pixel focal plane arrays, creating sub-mm imaging spectrometers. These
large format focal plane arrays (FPA) are ideally based on broad-band antenna
systems with a single beam per feed, allowing a tight sampling of the focal
plane. These FPA are typically coupled to reflector systems with large Focal
distance to Diameter ratio (F/D) $>$3\. Such antennas can also be used as an
alternative for CMB missions using multi-color pixels, such as PolarBear which
is currently using a sinuous antenna [4, 5].

We demonstrate the leaky lens antenna as an ideal candidate for these
applications. This antenna, first demonstrated by Neto et al. [6, 7] at
frequencies up to 70 GHz, is characterised by a 1:3 bandwidth, frequency
independent, linearly polarized beams and a high, frequency independent
aperture efficiency [8]. The first experimental demonstration at 1.55 THz was
presented in 2017 by Bueno et al,[9]. These experiments demonstrated super-THz
radiation detection with excellent detector sensitivity, but with a limited
aperture efficiency $\eta_{ap}=0.24$, evaluated only at a single frequency.
The low aperture efficiency was due to a combination of i) radiation loss in
the antenna-detector ground plane, ii) a slight misalignment between the feed
and the lens and iii) a mismatch, due to fabrication errors, between
transmission line and feed. In this paper we demonstrate the fabrication and
experimental validation of a leaky lens antenna coupled MKID (Microwave
Kinetic Inductance Detector [10]), which does not suffer from these
constraints, over a broad frequency band of 300-900 GHz. We measure the
absolute coupling efficiency and far field patterns at 3 frequencies: 350, 650
and 850 GHz using narrow band-pass quasioptical filters. From these we obtain
an aperture efficiency $\eta_{ap}\approx 0.4$ using a Si lens without anti-
reflection coating, which is in good agreement with the predicted performance.

Figure 1: (a) Image of the chip center with the tapered antenna slot in the
top and the MKID running along the center, coupled to the readout line at the
bottom. Additional blind MKIDs are located on the left and right. (b) The
leaky wave slot of the antenna on the SiN membrane. (c) Coupling structure to
the narrow section of the MKID. The Al center line is shorted via a galvanic
connection to the NbTiN ground plane. (d) Narrow MKID section at the membrane
edge, marked by the dashed blue line. The $1\ \mu$m thick SiN membrane is
transparent and therefore appears slightly blue in the back and front
illuminated picture. Below the dashed line, the MKID continues on the SiN
layer with Si substrate below. (e) The transition between narrow and wide
section of the MKID. (f) A schematic of the full 3D assembly including the
chip, spacer chip and lens. The total stack height is 4.325 mm. The radiation
launched in the vacuum gap is refracted at the Si/vacuum interface. (g) The
simulated radiation efficiency and its component efficiencies, with
$\eta_{rad}=\eta_{mismatch}\eta_{FtB}\eta_{pol}\eta_{reflection}$.

##  II Device and Antenna Design

The leaky lens antenna [6] is a lens feed consisting of a leaky wave slot on a
thin membrane that is kept at an electrically small distance from the
dielectric lens, in order to obtain a directive radiation pattern inside the
dielectric which illuminates efficiently the top part of the lens. The leaky
wave slot antenna is optimized via simulations in CST with an infinite
dielectric medium [11], following the same approach as in [8]. The resulting
tapered slot combines a high aperture efficiency which is almost constant over
a 1:3 bandwidth and is coupled to a single planar feed. The optimization of
the slot to create a frequency independent aperture efficiency results in a
tapered slot with a length of $990\mu$m, slot tapering angle of $15.8^{\circ}$
and a vacuum-gap between antenna and silicon lens of $h=13.8\ \mu$m. The
tapering angle is a compromise between polarization purity and antenna beam
width, while the slot length is determined by the lowest frequency. The lens
is optimized using a custom Physical Optics (PO) code to create diffraction
limited beams with maximized directivity. We use the custom PO code, because
the lens size (62 $\lambda$ inside silicon) is too large to allow for full-
wave analysis in CST Microwave Studio. This also prevented us from treating
multiple reflections inside the lens. This yields a truncated, hyper-
hemispherical lens shape with radius of curvature $R=2.95\ $mm, extension
length $L=1.0\ $mm and truncation angle $\theta=53.1^{\circ}$ without anti-
reflection coating. The diameter of the truncated lens is $5.5\ $mm. We
achieve a simulated radiation efficiency $\eta_{rad}$ between $0.45$ and
$0.55$ across the whole frequency band as shown in Fig. 1g. The dominant
effect is the reflection efficiency due to the lack of anti-reflection coating
at the Si-lens surface. Secondly, the polarization efficiency, representing
the fraction of the total emitted radiation which is in the co-polarized
state, ranges from $94\%$ at low frequencies to $82\%$ at high frequencies due
to the relatively large thickness of the vacuum gap which constitutes an
inherent tradeoff between polarization purity and directivity. All other
contributions to the radiation efficiency, i.e. the feed match
($\eta_{match}$) and front-to-back ratio ($\eta_{FtB}$), are high. The
simulated taper efficiency $\eta_{tap}$and aperture efficiency
$\eta_{ap}=\eta_{rad}\cdot\eta_{tap}$ are also given in Fig. 1 g. We find that
$\eta_{ap}\approx$ 0.4 over the whole frequency band. This is slightly worse
than that of corrugated feedhorns, but here we achieve a much larger
bandwidth, and comparable to Sinuous antenna’s. The development of a broad
band anti-reflection coating would increase the aperture efficiency to values
close to 0.6.

We couple the leaky wave slot to a narrow coplanar waveguide (CPW) feed with a
NbTiN ground and an Al central line (Fig. 1 a-e)). This narrow CPW line is
also the last section of an antenna coupled Microwave Kinetic Inductance
detector (MKID), which is a $\lambda$/4 CPW resonator with a length of several
mm and a resonant frequency of 4.5 GHz. For details on this specific MKID
design we refer to Ref.[12]. Radiation coupled to the antenna (panel b) will
propagate into the NbTiN-Aluminium CPW line (panel c), where it will be
absorbed only in aluminium central line. This is because the frequency of the
radiation (300-900 GHz) exceeds the gap frequency of Al ($F_{\Delta,Al}\approx
93\ $GHz for $T_{c}\approx 1.28\ $K) but does not exceed the gap frequency of
the NbTiN ground plane ($F_{\Delta}\approx 1.1\ $ THz). The NbTiN thereby
provides a lossless ground plane for the antenna as well as for the NbTiN-
Aluminium CPW line, eliminating ground plane loss as in Ref. [9]. The
radiation absorbed in the aluminium changes the complex surface impedance
which results in a slight shift in resonant frequency of the MKID. This shift
is detected by a sending a single frequency readout signal at the unperturbed
resonant frequency of the MKID through the readout line visible in the bottom
of Fig. 1a. The frequency shift causes a change in transmitted phase in the
readout signal that is measured using a homodyne measurement scheme [10]. The
narrow CPW section has a line width $w=2\ \mu$m and slot width $2.4\ \mu$m, a
thickness of 50 nm and a length of 1.5mm. It is fabricated on top of $1\ \mu$m
SiN. Underneath part of this line, the Si is removed to create a free standing
membrane, needed for a correct operation of the leaky lens antenna (see
Fig.1f)). Due to its small width, this section has low radiation loss and a
high kinetic inductance fraction, relevant for efficient detection of incoming
THz radiation. The remaining resonator therefore consists of a wide CPW of
$w=6\ \mu$m and $s=16\ \mu$m made of NbTiN on bare Si-substrate, with a
transition between the sections and step from SiN to Si as shown in Fig. 1e).
The wide section of the MKID then ends in a coupling structure that couples
the resonator weakly to the CPW readout line. We choose a very low coupling
quality factor of $Q_{c}=5000$, because the power absorbed in the device
during the beam pattern measurements is $\sim 1\ $nW, resulting in a very low
internal Q factor of the MKID. Additionally, two meshed layers of low-Tc
($0.65\ $K) $\beta$-phase Ta are located on the front and backside of the
chip, to reduce stray radiation inside the chip[13]. Two blind MKIDs are
located close to the antenna-coupled MKID (see Fig. 1a)) as reference
detectors.

##  III Fabrication

We process the chip from a 375 um thick 4-inch Si wafer coated on both sides
with a $1\ \mu$m thick, low tensile stress ($\sim 250\ $MPa) SiN layer,
deposited using low pressure chemical vapor deposition (LPCVD) (see Fig 1f)).
The fabrication is similar to the method presented by Bueno et al. [9] with
steps as follows: i) We etch the SiN on the chip device side with reactive ion
etching (RIE) using $35\%$ SF6 and $65\%$ O2 to create a sloped edge for the
MKID wide section (Fig. 1c)). ii) We etch the SiN on the wafer backside, to
create a hard mask for later Si removal. iii) We deposit a $350\ $nm thick,
low tensile stress NbTiN layer on the device side using reactive sputtering of
a NbTi target in a Nitrogen-Argon atmosphere, with the sloped edge of the SiN
providing good step coverage for the NbTiN [14]. We pattern the NbTiN layer
defining the antenna slot, the MKID resonator wide section and narrow section
ground plane as well as the microwave readout line. We then etch the NbTiN
using the recipe of step (i). iv) We sputter a 40 nm thick $\beta$-phase Ta
layer on the device side of the wafer, which is then patterned and etched to
form an absorbing mesh [13]. v) We spin coat $1\ \mu$m polyimide LTC9505 on
the wafer, which is patterned and cured into dielectric supports for
superconducting bridges along the readout line. vi) We sputter a $50\ $nm
thick Al layer on the device side, which is then patterned and wet etched with
a commercially available Al etchant [15] to define the central line of the
narrow MKID section as well as the bridges along the readout line. vii) The Si
wafer is etched in a KOH bath, with the device side protected with a
commercial protection tool, to create the membrane opening. viii) After the
KOH etch, we deposit 40 nm $\beta$-phase Ta on the wafer backside analogous to
step (iv).

To create the 3D structure of the lens-antenna shown in 1g), we also need a
spacer chip and Si-lens. The spacer chip is fabricated on a separate 350 um
thick 4-inch Si wafer with SiN on both sides, which is patterned. The wafer is
subsequently etched twice in a KOH bath: Once on the lens facing side to
create 4 square holes for the alignment pins, which will ensure the alignment
between the lens and the antenna, and once on the antenna facing side to
create the $14\ \mu$m deep features of the vacuum gap. The Si-lens was made by
a commercial partner (Tydex) according to the specifications of our design.

We assemble the device chip, spacer chip and lens in a 2-step procedure:
First, the chip and spacer chip are aligned and clamped together using
Aluminum alignment pins. The lens is then aligned using a microscope based,
home-made tool and glued on top of the spacer wafer using Loctite 406
cyanoacrylate glue, where we use markers on the spacerwafer from the KOH etch
for alignment. This glue is chosen for its extremely low viscosity, resulting
in a glue gap $<1\ \mu$m.

Figure 2:  (a) The CAD-model of the cryostat optics and the hot-source in the
configuration used for beampattern measurements. (b) A close-up of the source
assembly. (c) Filter transmission as function of frequency for the 3
individual filters used in the experiment.

##  IV Measurements and Results

###  IV-A Beam Patterns

Measurements of the lens-antenna beam pattern were carried out in the setup
shown in Fig. 2a), where the chip carrying the detector coupled lens-antenna
is located on the cold stage of a cryogen free cryocooler equipped with a
3He/4He sorption fridge with direct optical access to a source located on an
optical table at room temperature [16]. The cryostat reaches a base
temperature of $T\sim 265\ $mK with a holdtime of $\sim 32\ $ hours. The
opening angle of $\theta=\pm 20\deg$ is limited by the rotatable polarizer
mounted on the outside of the cryostat’s vacuum window.

The source assembly consists of three main parts: First, a glowbar (Scitec
Hawkeye IR-19) operated at $T=1150^{\circ}$ is used as a multimoded blackbody
source, located at one focal point of an elliptical mirror with a small
aperture located at the second focus to create a uniform source pattern. We
independently validated the source uniformity using a commercial bolometer and
rotation stage. Secondly, an eight-bladed chopper, coated with radiation
absorber[17, 18], is positioned between the two foci. This allows a
temperature modulation between the hot source and chopper blade at room
temperature. Finally, the glowbar/chopper combination is mounted on an XY-
scanner, allowing the measurement of 2D beam patterns of the antenna in the
plane of the scanner. Additionally, the cryostat facing side of the source
assembly has multiple blackened[18] plates mounted to it, reducing unwanted
reflections.

We perform 3 measurements, each with a different band pass filter on the 1K
stage of the cooler. The band pass characteristics are shown in Fig 2c). They
have a sufficiently small bandwidth so that the patterns are not smeared
significantly with respect to the step-size of our XY-scanner.

Figure 3: The normalized beam patterns in E- and H-plane for all measured
frequencies. Insets show $20^{\circ}\times 20^{\circ}$ maps of the measured 2D
beam patterns (left to right: $350$, $650$, $850\ $GHz) with a color scale
from $-20\ $dB (dark) to $0\ $dB (light) and a $-3\ $dB contour line.
Additionally, the orientation of the antenna slot compared to the measurement
planes is shown in the top-left.

We extract the MKID response as a function of the source position using time-
domain Fourier analysis, extracting the peak at the chopper frequency of
$f=254\ $Hz for all scanner positions in a step-and-integrate scan mode. The
resulting beam patterns are shown in Fig. 3 and compared to simulations
carried out at the band pass center frequency. The beam patterns, as simulated
using the in-house PO tool, do not take into account multiple reflections that
occur inside the lens surface because of the absence of an anti-reflection
coating. Multiple reflections can result in a decrease in directivity and an
increase in sidelobe-level which is difficult to predict and can change
significantly with frequency.

All three frequency bands show excellent agreement between experiment and
simulation for the main lobe down to the $-10\ $dB level. The beams are
diffraction limited as predicted by simulation, with a full width at half
maximum (FWHM) of $3.7^{\circ}$ at $850\ $GHz and $8.3^{\circ}$ at $350\ $GHz
due to the small lens diameter. The leaky-lens antenna is designed for future
multi-pixel spectroscopic instruments with closely packed Focal Plane Arrays.
Inherent to such designs is that the leaky-lenses are diffraction limited with
maximized directivity, as mentioned in Section II. A diffraction limited beam
can be described with an Airy pattern with a maximum theoretical Gaussicity of
82%. The Gaussicity of the beams is about 72% for all three frequencies and
are close to the theoretical maximum.

The measured side lobe levels are around $2$ to $5\ $dB higher than simulated,
which we tentatively attribute to multiple reflections inside the lens
increasing the effective detected power at large incident angles. An
additional possible contribution might be the detection of residual stray-
light via surface waves in the chip[13].

We measure the optical efficiency in a separate cryostat which provides a
highly controlled environment and is described in detail in Ref. [19]. For
this experiment the chip is mounted in a light-tight box on the cryostats cold
stage, an adiabatic dilution refrigerator (ADR), operated at a stable $120\
$mK. An opening in the light-tight box provides a field-of-view through
multiple optical filters to a thermal calibration source with a temperature
range of $2.7\ $K - $40\ $K. The optical filterstacks correspond to the bands
of the beam pattern measurements (see Fig. 4). However, their out-of-band
rejection is much larger, which is needed to sufficiently reject low frequency
radiation from the black body source operating at low temperatures.

The angular throughput to the source is limited by an aperture plate to an
angle of $\pm 10^{\circ}$. This is large enough to couple almost to the whole
beam at $350\ $GHz, but avoids straylight from large angles.

Figure 4:  (a) The transmission curves for the three optical filterstacks; (b)
The measured optical coupling efficiency, together with the simulated values
of the optical coupling efficiency and its contributions, the radiation
efficiency and spillover efficiency between the lens-antenna and the
calibration source.

###  IV-B Optical Efficiency

We determine the optical coupling efficiency $\eta_{opt}(\nu_{0})$ of the MKID
coupled lens-antenna at a frequency $\nu_{0}$ experimentally by measuring the
noise-equivalent power $NEP_{exp}$ at a source temperature $T_{s}$:

$\eta_{opt}(\nu_{0})=\frac{\int 2P_{s,\nu}h\nu d\nu+\int\frac{4\Delta
P_{s,\nu}}{\eta_{pb}}d\nu}{NEP_{exp}^{2}-\int 2P_{s,\nu}h\nu
F_{\nu}O_{\nu}d\nu}$ (1)

with the superconductor bandgap $\Delta$, the pair-breaking efficiency
$\eta_{pb}\sim 0.4$ for our film [20], the filter transmission $F_{\nu}$,
occupation number $O_{\nu}$ and the single-moded source power after filter
transmission $P_{s,\nu}$. A detailed description and derivation of this
equation and its components can be found in Ferrari et al. [21].

This method is valid for narrow frequency bands and photon-noise limited
MKIDs, and allows us to obtain precise absolute measurements of the optical
efficiency at multiple frequencies across the whole band of the lens-antenna.
For our Al based hybrid MKIDs, photon noise limited performance is proven
experimentally, if the noise spectrum is white with a roll-off consistent with
the quasiparticle lifetime. This is true for the detector presented at
absorbed power levels $P>10\ $fW. While it is in principle possible to measure
the response over the whole band using a Fourier-transform spectrometer (FTS)
[9], this only provides a relative power measurement and not an absolute
validation of the optical coupling efficiency of the detector .

The optical coupling efficiency measurements are carried out at multiple
$T_{s}$, corresponding to a power range between $10\ $fW and $1\ $pW, which
are averaged to obtain the values shown in Fig. 4 in comparison with
simulation. for a detailed description of the experimental procedure we refer
to Ref.[12]. As shown in bottom panel of Fig. 4 we find excellent agreement
between the measured and simulated optical coupling efficiency
$\eta_{opt}=\eta_{rad}\cdot\eta_{Spillover}$ for all frequency bands. Here
$\eta_{spillover}$ describes the spillover between the lens-antenna beam and
the source, which is close to unity except for the lowest frequency. Hence
this result demonstrates that the measured $\eta_{rad}$ is in good agreement
with our model calculation. Together with the good agreement of the patterns
which validate the model calculation of $\eta_{tap}$, we can conclude that the
aperture efficiency $\eta_{ap}=\eta_{rad}\cdot\eta_{tap}$ is well described by
our model, and given by $\eta_{ap}\approx 0.4$ over the whole frequency band
(see Fig.1). Only at $850\ $GHz we see a slightly higher coupling, possibly
due to a small stray light contribution: The beam width is much smaller than
the source aperture at this frequency.

##  V Conclusion

In conclusion, we have measured beam patterns and optical efficiency of a
leaky lens antenna over a 1:3 bandwidth around $350\ $GHz, $650\ $GHz and
$850\ $GHz. We infer an aperture efficiency of the lens-antenna $\eta_{ap}$ =
0.4 over the whole frequency band, based on measurements of the absolute
coupling efficiency between the detector and a thermal calibration source
together with a good agreement between measured and simulated beam patterns.
The aperture efficiency is limited primarily by the antenna radiation
efficiency, which is low due to the absence of an anti-reflection coating. The
agreement between simulation and experiment show the validity of the design
process, allowing for the design of optimised lens-antennas for on-chip
spectrometer or CMB applications. The development of a broadband anti-
reflection coating would allow for a significant increase in optical
efficiency by up to 30%.f Given the wide angular extent of the primary field
in the dielectric an anti-reflection coating based upon meta-materials is
extremely difficult to machine due to the steepness of the lens. An
interesting alternative would be the use of stacks of foam as proposed by Ref.
[22] or the use of thermal spray coatings [4], which allow for a more
practical approach towards a conformal coating. Changing the lens geometry and
feed position could allow for frequency stable patterns[23], with lower
directivity, which would allow for a more efficient broad-band coupling to a
reflector system.

## References

  * [1] A. Endo, P. van der Werf, R. M. J. Janssen, P. J. de Visser, T. M. Klapwijk, J. J. A. Baselmans, L. Ferrari, A. M. Baryshev, and S. J. C. Yates, “Design of an integrated filterbank for deshima: On-chip submillimeter imaging spectrograph based on superconducting resonators,” _Journal of Low Temperature Physics_ , vol. 167, no. 3, pp. 341–346, May 2012. [Online]. Available: https://doi.org/10.1007/s10909-012-0502-1
  * [2] E. Shirokoff, P. S. Barry, C. M. Bradford, G. Chattopadhyay, P. Day, S. Doyle, S. Hailey-Dunsheath, M. I. Hollister, A. Kovács, H. G. Leduc, C. M. McKenney, P. Mauskopf, H. T. Nguyen, R. O’Brient, S. Padin, T. J. Reck, L. J. Swenson, C. E. Tucker, and J. Zmuidzinas, “Design and performance of superspec: An on-chip, kid-based, mm-wavelength spectrometer,” _Journal of Low Temperature Physics_ , vol. 176, no. 5, pp. 657–662, Sep 2014. [Online]. Available: https://doi.org/10.1007/s10909-014-1122-8
  * [3] G. Cataldo, W.-T. Hsieh, W.-C. Huang, S. H. Moseley, T. R. Stevenson, and E. J. Wollack, “Micro-spec: an ultracompact, high-sensitivity spectrometer for far-infrared and submillimeter astronomy,” _Appl. Opt._ , vol. 53, no. 6, pp. 1094–1102, Feb 2014. [Online]. Available: http://ao.osa.org/abstract.cfm?URI=ao-53-6-1094
  * [4] Y. Inoue _et al._ , “POLARBEAR-2: an instrument for CMB polarization measurements,” _Proc. SPIE Int. Soc. Opt. Eng._ , vol. 9914, p. 99141I, 2016\.
  * [5] R. O’Brient, P. Ade, K. Arnold, J. Edwards, G. Engargiola, W. L. Holzapfel, A. T. Lee, M. J. Myers, E. Quealy, G. Rebeiz, P. Richards, and A. Suzuki, “A dual-polarized broadband planar antenna and channelizing filter bank for millimeter wavelengths,” _Applied Physics Letters_ , vol. 102, no. 6, p. 063506, 2013. [Online]. Available: https://doi.org/10.1063/1.4791692
  * [6] A. Neto, “Uwb, non dispersive radiation from the planarly fed leaky lens antenna 2014; part 1: Theory and design,” _IEEE Transactions on Antennas and Propagation_ , vol. 58, no. 7, pp. 2238–2247, July 2010.
  * [7] A. Neto, S. Monni, and F. Nennie, “Uwb, non dispersive radiation from the planarly fed leaky lens antenna 2014;part ii: Demonstrators and measurements,” _IEEE Transactions on Antennas and Propagation_ , vol. 58, no. 7, pp. 2248–2258, July 2010.
  * [8] O. Yurduseven, N. L. Juan, and A. Neto, “A dual-polarized leaky lens antenna for wideband focal plane arrays,” _IEEE Transactions on Antennas and Propagation_ , vol. 64, no. 8, pp. 3330–3337, Aug 2016.
  * [9] J. Bueno, O. Yurduseven, S. J. C. Yates, N. Llombart, V. Murugesan, D. J. Thoen, A. M. Baryshev, A. Neto, and J. J. A. Baselmans, “Full characterisation of a background limited antenna coupled kid over an octave of bandwidth for thz radiation,” _Applied Physics Letters_ , vol. 110, no. 23, p. 233503, 2017. [Online]. Available: https://doi.org/10.1063/1.4985060
  * [10] P. K. Day, H. G. LeDuc, B. A. Mazin, A. Vayonakis, and J. Zmuidzinas, “A broadband superconducting detector suitable for use in large arrays,” _Nature_ , vol. 425, p. 817, Oct 2003. [Online]. Available: http://dx.doi.org/10.1038/nature02037
  * [11] N. Llombart and A. Neto, “Thz time-domain sensing: The antenna dispersion problem and a possible solution,” _IEEE Transactions on Terahertz Science and Technology_ , vol. 2, no. 4, pp. 416–423, July 2012.
  * [12] R. M. J. Janssen, J. J. A. Baselmans, A. Endo, L. Ferrari, S. J. C. Yates, A. M. Baryshev, and T. M. Klapwijk, “High optical efficiency and photon noise limited sensitivity of microwave kinetic inductance detectors using phase readout,” _Applied Physics Letters_ , vol. 103, no. 20, p. 203503, 2013\. [Online]. Available: https://doi.org/10.1063/1.4829657
  * [13] S. J. C. Yates, A. M. Baryshev, O. Yurduseven, J. Bueno, K. K. Davis, L. Ferrari, W. Jellema, N. Llombart, V. Murugesan, D. J. Thoen, and J. J. A. Baselmans, “Surface wave control for large arrays of microwave kinetic inductance detectors,” _IEEE Transactions on Terahertz Science and Technology_ , vol. 7, no. 6, pp. 789–799, Nov 2017.
  * [14] D. J. Thoen, B. G. C. Bos, E. A. F. Haalebos, T. M. Klapwijk, J. J. A. Baselmans, and A. Endo, “Superconducting nbtin thin films with highly uniform properties over a ${varnothing}$ 100 mm wafer,” _IEEE Transactions on Applied Superconductivity_ , vol. 27, no. 4, pp. 1–5, June 2017\.
  * [15] TechniEtch Al80 Al etchant from Microchemicals GmbH.
  * [16] S. Hähnle, J. Bueno, R. Huiting, S. J. C. Yates, and J. J. A. Baselmans, “Large angle optical access in a sub-kelvin cryostat,” _Journal of Low Temperature Physics_ , May 2018. [Online]. Available: https://doi.org/10.1007/s10909-018-1940-1
  * [17] T. O. Klaassen, M. C. Diez, J. H. Blok, C. Smorenburg, K. J. Wildeman, and G. Jakob, “Optical characterization of absorbing coatings for sub-millimeter radiation,” in _Proceedings of the 12th International Symposium on Space THz Technology, Feb. 2001_ , 2001, pp. 400–409.
  * [18] Stycast2850FT loaded with 1mm diameter SiC grains.
  * [19] J. Baselmans, S. Yates, P. Diener, and P. de Visser, “Ultra low background cryogenic test facility for far-infrared radiation detectors,” _Journal of Low Temperature Physics_ , vol. 167, no. 3, pp. 360–366, May 2012. [Online]. Available: https://doi.org/10.1007/s10909-012-0511-0
  * [20] T. Guruswamy, D. J. Goldie, and S. Withington, “Quasiparticle generation efficiency in superconducting thin films,” _Superconductor Science Technology_ , vol. 27, no. 5, p. 055012, May 2014.
  * [21] L. Ferrari, O. Yurduseven, N. Llombart, S. J. C. Yates, J. Bueno, V. Murugesan, D. J. Thoen, A. Endo, A. M. Baryshev, and J. J. A. Baselmans, “Antenna coupled mkid performance verification at 850 ghz for large format astrophysics arrays,” _IEEE Transactions on Terahertz Science and Technology_ , vol. 8, no. 1, pp. 127–139, Jan 2018.
  * [22] A. Nadolski, A. M. Kofman, J. D. Vieira, P. A. R. Ade, Z. Ahmed, A. J. Anderson, J. S. Avva, R. Basu Thakur, A. N. Bender, B. A. Benson, J. E. Carlstrom, F. W. Carter, T. W. Cecil, C. L. Chang, J. F. Cliche, A. Cukierman, T. de Haan, J. Ding, M. A. Dobbs, D. Dutcher, W. Everett, A. Foster, J. Fu, J. Gallichio, A. Gilbert, J. C. Groh, S. T. Guns, R. Guyser, N. W. Halverson, A. H. Harke-Hosemann, N. L. Harrington, J. W. Henning, W. L. Holzapfel, N. Huang, K. D. Irwin, O. B. Jeong, M. Jonas, A. Jones, T. S. Khaire, M. Korman, D. L. Kubik, S. Kuhlmann, C. L. Kuo, A. T. Lee, A. E. Lowitz, S. S. Meyer, D. Michalik, J. Montgomery, T. Natoli, H. Nguyen, G. I. Noble, V. Novosad, S. Padin, Z. Pan, J. Pearson, C. M. Posada, W. Quan, A. Rahlin, J. E. Ruhl, J. T. Sayre, E. Shirokoff, G. Smecher, J. A. Sobrin, A. A. Stark, K. T. Story, A. Suzuki, K. L. Thompson, C. Tucker, K. Vanderlinde, G. Wang, N. Whitehorn, V. Yefremenko, K. W. Yoon, and M. R. Young, “Broadband anti-reflective coatings for cosmic microwave background experiments,” in _Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX_ , ser. Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, vol. 10708, Jul 2018, p. 1070843.
  * [23] D. F. Filipovic, S. S. Gearhart, and G. M. Rebeiz, “Double-slot antennas on extended hemispherical and elliptical silicon dielectric lenses,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 41, no. 10, pp. 1738–1749, Oct 1993.


# Integrated Ray-Tracing and Coverage Planning Control
using Reinforcement Learning

Savvas Papaioannou, Panayiotis Kolios, Theocharis Theocharides,
Christos G. Panayiotou and Marios M. Polycarpou The authors are with the KIOS
Research and Innovation Centre of Excellence (KIOS CoE) and the Department of
Electrical and Computer Engineering, University of Cyprus, Nicosia, 1678,
Cyprus. {papaioannou.savvas, pkolios, ttheocharides, christosp,
<EMAIL_ADDRESS>
###### Abstract

In this work we propose a coverage planning control approach which allows a
mobile agent, equipped with a controllable sensor (i.e., a camera) with
limited sensing domain (i.e., finite sensing range and angle of view), to
cover the surface area of an object of interest. The proposed approach
integrates ray-tracing into the coverage planning process, thus allowing the
agent to identify which parts of the scene are visible at any point in time.
The problem of integrated ray-tracing and coverage planning control is first
formulated as a constrained optimal control problem (OCP), which aims at
determining the agent’s optimal control inputs over a finite planning horizon,
that minimize the coverage time. Efficiently solving the resulting OCP is
however very challenging due to non-convex and non-linear visibility
constraints. To overcome this limitation, the problem is converted into a
Markov decision process (MDP) which is then solved using reinforcement
learning. In particular, we show that a controller which follows an optimal
control law can be learned using off-policy temporal-difference control (i.e.,
Q-learning). Extensive numerical experiments demonstrate the effectiveness of
the proposed approach for various configurations of the agent and the object
of interest.

##  I Introduction

The scientific advancements of the last decade in the areas of robotics,
control, and machine learning have spurred an increased interest in the
utilization of autonomous systems in various application domains including
search-and-rescue [1, 2, 3], security [4, 5, 6], and monitoring [7, 8, 9].

Coverage planning [10] is amongst the most important tasks found in many
applications including infrastructure inspection, automated maintenance, area
search, and surveillance, and thus plays a pivotal role in designing and
executing automated missions using autonomous mobile agents. In coverage
planning we are interested in finding a trajectory which allows an autonomous
agent (e.g., a mobile robot) to observe (or cover) with its sensor every
point/region within a specified area of interest. More specifically, during an
automated coverage mission, the mobile agent must autonomously plan its
coverage trajectory (e.g., determine its control inputs; possibly over a
planning horizon), which allow the efficient coverage of the area of interest,
while satisfying certain kinematic and sensing constraints. As discussed in
more detail in Sec. II, over the years, a plethora of coverage planning
approaches have been proposed in the literature. However, there is still room
for improvement until this technology reaches the required level of maturity.
Specifically, the vast majority of coverage planning approaches mainly
accounts for agents with uncontrollable and fixed sensors, and without
considering the agent’s kinematic and sensing constraints. This allows the
coverage planning problem to be reduced to a conventional path planning
problem which can further be transformed to some form of the traveling
salesman problem i.e., seeking to find the shortest path that passes through a
finite set of points.

Complementary to the existing literature, in this work we propose an
integrated ray-tracing and coverage planning control approach where a mobile
agent equipped with a controllable camera sensor with limited sensing domain
(i.e., finite sensing range and angle of view), autonomously generates its
trajectory such that the total surface area of a known object/area of interest
is optimally covered. In the proposed approach we integrate ray-tracing into
the planning process in order to determine the visible parts of the scene
through the agent’s camera sensor. The problem of integrated ray-tracing and
coverage planning control is first formulated as constrained optimal control
problem (OCP), which however is very challenging to be solved due to the
existence of non-convex and non-linear constraints. For this reason, the
problem is transformed into a finite Markov decision process (MDP) and is
solved using reinforcement learning. Specifically, the contributions of this
work are the following:

  * •

We propose a coverage planning control approach which allows an autonomous
mobile agent to efficiently cover the total surface area of an object of
interest by jointly optimizing its kinematic and camera control inputs. The
proposed approach integrates ray-tracing into the coverage planning process in
order to simulate the physical behavior of light-rays, thus enabling the agent
to identify the visible parts of the scene through its camera field-of-view.

  * •

The problem of integrated ray-tracing and coverage control is initially
formulated as a constrained optimal control problem and subsequently converted
into a finite Markov decision process (MDP) which can be solved efficiently
using reinforcement learning. Specifically, we show that off-policy temporal-
difference control (i.e., Q-learning) can be utilized to learn a coverage
controller which follows an optimal control law.

  * •

Extensive numerical experiments demonstrate the effectiveness of the proposed
approach for various configurations of the agent and the object of interest.

The rest of the paper is organized as follows. Section II summarizes the
related work on coverage planning control. Then, Section III develops the
system model, Section IV formulates the problem tackled, and Section V
discusses the details of the proposed coverage planning approach. Finally,
Section VI evaluates the proposed approach and Section VII concludes the paper
and discusses future work.

##  II Related Work

Initial works on coverage planning control [11, 12, 13, 14, 15] proceed by
decomposing the environment into a finite number of disjoint cells, and then
utilize a path-planning algorithm [16] to find the best path that passes from
every cell. The works in [17, 18, 19, 20] investigate the coverage planning
problem with multiple mobile agents which however exhibit fixed and
uncontrollable sensor footprints. Different variations of the coverage
planning problem are investigated in [21, 22, 23] i.e., in [21, 22] the
objective is to derive a decentralized, control law which guides a team of
agents to an optimal configuration that maximizes coverage, and in [23] the
coverage planning problem is investigated with multiple heterogeneous agents.
The coverage planning problem is also investigated with the use of unmanned
aerial vehicles (UAVs) [24, 25, 26, 27]. Specifically, the authors in [24]
propose an exact cellular decomposition method to plan the coverage path of
UAVs equipped with fixed sensors in a polygon area, whereas in [25, 26] the
problem is investigated under similar assumptions with fixed-wing UAVs. In
[27] the coverage planning problem is investigated in the presence of
visibility constraints. The authors in [28, 29, 30] propose a search planning
framework based on mixed integer quadratic programming (MIQP), which allows an
autonomous UAV agent to cover in 3D specific objects of interest. Finally, in
[31, 32] learning based coverage planning techniques are investigated.
Specifically, in [31], an end-to-end deep reinforcement learning coverage
planning approach is proposed for an autonomous UAV agent with battery
constraints, whereas in [32] the authors design a multi-agent coverage
planning approach based on supervised imitation learning, for searching a
finite number of cells within a bounded environment.

##  III System Model

###  III-A Agent Kinematic Model

We assume that an autonomous mobile agent, which in this work is represented
by a point-mass object, evolves inside a bounded planar workspace
$\mathcal{W}\subset\mathbb{R}^{2}$ according to a discrete-time kinematic
model of the form $x_{t}=f(x_{t-1},u_{t-1})$, given by:

$x_{t}=x_{t-1}+d_{R}\begin{bmatrix}\cos(\vartheta)\\\
\sin(\vartheta)\end{bmatrix},~{}t\geq 1,~{}x_{0}=\bar{x}$ (1)

where $x_{t}=[\text{x},\text{y}]^{\top}\in\mathbb{R}^{2}$ is the kinematic
state of the agent at time $t$, which is given by the agent’s position in 2D
cartesian coordinates, and $\bar{x}$ is the agent’s initial state. The agent’s
kinematic control input is given by $u_{t}=[d_{R},\vartheta]$, where
$d_{R}\in\mathbb{R}$ is the radial displacement, and $\vartheta\in[0,2\pi)$
denotes the agent’s heading. Finally, we assume that the agent’s state is
bounded i.e., $x_{t}\in\mathcal{X},\forall t$ and the control input is
constrained according the agent’s kinematic capabilities i.e.,
$u_{t}\in\mathcal{U},\forall t$.

###  III-B Agent Sensor Model

Figure 1: The figure illustrates the agent’s camera FOV state for 5 rotation
angles, where $x_{t}$ is the agent’s kinematic state, $h$ is the camera
sensing range, $\varphi$ is the FOV’s angle of view, $s_{\ell}$ and $b_{\ell}$
are the FOV side length and base length respectively, and
$\mathcal{R}_{x_{t},\theta_{t}}$ denotes the set of camera rays at time $t$
when the camera FOV has been rotated by $\theta_{t}$.

The autonomous mobile agent is equipped with a controllable camera sensor with
limited sensing domain (i.e., finite sensing range and angle of view), which
uses for observing the surrounding environment. The projected camera field-of-
view (FOV) is modeled in this work as an isosceles triangle parameterized by
its angle $\varphi$ at the apex and its height $h$, which are used to model
the FOV’s angle of view and sensing range respectively, as illustrated in Fig.
1. The FOV’s side length $(s_{\ell})$ and base length $(b_{\ell})$ are
computed according to:

$\displaystyle s_{\ell}$
$\displaystyle=h\times\text{cos}(\varphi/2)^{-1},\quad
b_{\ell}=2s_{\ell}\times\text{sin}(\varphi/2)$ (2)

and thus the matrix $\mathcal{F}_{o}$ which contains the vertices of a
downward facing FOV with the apex at the origin is given by:

$\mathcal{F}_{o}=\begin{bmatrix}0&-b_{\ell}/2&b_{\ell}/2\\\ 0&-h&-h\\\
\end{bmatrix}$ (3)

The agent’s camera FOV is controllable i.e., it can be rotated around the
agent’s position $x_{t}$ by an angle $\theta_{t}\in\Theta\subseteq[0,2\pi)$ at
time $t$ according to: $\mathcal{F}_{t}=Q(\theta_{t})\mathcal{F}_{o}+x_{t}$,
where $\mathcal{F}_{t}$ is the rotated camera FOV state, and $Q(\theta_{t})$
is a 2D rotation matrix which is given by:

$Q(\theta_{t})=\begin{bmatrix}\text{cos}(\theta_{t})&\text{sin}(\theta_{t})\\\
-\text{sin}(\theta_{t})&\text{cos}(\theta_{t})\end{bmatrix}$ (4)

Moreover, we assume that at each time-step $t$ a finite set of light-rays
(which are assumed to be straight), and which indicate the direction of the
propagation of light, enter the camera’s optical sensor and cause the matter
inside the agent’s FOV to be imaged. The set of light-rays entering the
camera’s optical sensor at time $t$ when the agent is at state $x_{t}$ and the
camera FOV has been rotated by an angle $\theta_{t}$ is denoted as
$\mathcal{R}_{x_{t},\theta_{t}}=\\{R_{1},..,R_{|\mathcal{R}_{x_{t},\theta_{t}}|}\\}$,
with $|.|$ to denote the set cardinality. The individual ray
$R_{i}\in\mathcal{R}_{x_{t},\theta_{t}}$ is given by the line-segment:
$R_{i}=\\{\rho_{i}+s(x_{t}-\rho_{i})|s\in[0,1]\\}$, where $x_{t}$ is the ray’s
end point on the optical sensor given by the agent’s position at time $t$ and
$\rho_{i}\in\mathbb{R}^{2}$ is a fixed point on the base of the camera’s FOV
and determines the ray’s origin, as illustrated in Fig. 1.

To summarize, given an initial kinematic state $x_{0}$, the agent’s trajectory
(i.e., a sequence of kinematic and FOV states) over a finite planning horizon
of length $T$ time-steps can be optimized to meet certain objective by
appropriately selecting the control inputs $\\{\hat{u}_{t}|t=0,..,T-1\\}$,
where $\hat{u}_{t}=[u_{t},\theta_{t}]$.

###  III-C Object/Area of Interest Model

The objective of our autonomous mobile agent is to optimally decide its
kinematic and camera FOV inputs $\\{\hat{u}_{t}|t=0,..,T-1\\}$ over a finite
planning horizon such that the total surface area of a bounded convex object
or area of interest $\mathcal{C}\subset\mathcal{W}$ is covered by its sensor’s
FOV. More precisely, we are interested in the coverage of the object’s
boundary $\partial\mathcal{C}$, as illustrated in Fig. 2, since we are
operating in a 2D environment. However, the problem formulation and the
proposed approach can be generalized in 3D space as well.

Figure 2: In the scenario above we are interested in finding the agent’s
optimal control inputs which result in the coverage of the object’s
($\mathcal{C}$) boundary $\partial\mathcal{C}$. $\Delta\mathcal{C}$ is the
piece-wise linear approximation of $\partial\mathcal{C}$, which is composed of
a finite number of line segments $L_{p}$ connecting the points
$p\in\mathcal{P}$.

The boundary $\partial\mathcal{C}$ of the object of interest $\mathcal{C}$ is
approximated in this work by a piece-wise linear model denoted as
$\Delta\mathcal{C}$, which is obtained by connecting together a number of
points $p\in\partial\mathcal{C}$ (uniformly) sampled from the boundary. The
number of points sampled is chosen according to the mission requirements
(e.g., required coverage precision). Specifically, we denote the finite set of
points sampled from the object’s boundary as
$\mathcal{P}=\\{p_{1},..,p_{|\mathcal{P}|}\\}\subset\partial\mathcal{C}$.
Subsequently, the piece-wise linear approximation $\Delta\mathcal{C}$ of the
boundary $\partial\mathcal{C}$ is composed of line segments:
$L_{p}=\\{p+r(\hat{p}-p)|r\in[0,1],\hat{p}\neq p\in\mathcal{P}\\}$, which have
been generated by connecting the points
$(p,\hat{p})_{p\neq\hat{p}}\in\mathcal{P}$, such that the resulting line
segments $L_{p}\forall p$ belong to the boundary of the convex hull of
$\mathcal{P}$, as shown in Fig. 2. As we discuss in more detail in Sec. IV,
the agent’s objective now becomes the coverage of all points $p\in\mathcal{P}$
on the object’s boundary.

##  IV Problem Formulation

The problem of integrated ray-tracing and coverage planning control tackled in
this work can now be formulated as an optimal control problem (OCP) shown in
Problem (P1).

Problem (P1):Optimal Controller
$\displaystyle~{}~{}~{}~{}\underset{\mathbf{\hat{U}}_{T}}{\arg\max}~{}\mathcal{J}_{\text{coverage}}$
(5a) subject to: $t=[1,..,T]$ $\displaystyle x_{t}=f(x_{t-1},u_{t-1})$
$\displaystyle\forall t$ (5b) $\displaystyle x_{0}=\bar{x}$ (5c)
$\displaystyle\exists~{}\tau\in[1,..,T]:\left(p\in\mathcal{F}_{\tau}\right)\wedge\mathds{1}_{G_{\tau}}(p)$
$\displaystyle\forall p\in\mathcal{P}$ (5d) $\displaystyle
x_{0},x_{t}\notin\mathcal{C}$ $\displaystyle\forall t$ (5e) $\displaystyle
x_{0},x_{t}\in\mathcal{X}$ $\displaystyle\forall t$ (5f)
$\displaystyle\hat{u}_{t}=[u_{t},\theta_{t}]\in\mathcal{U}\times\Theta$
$\displaystyle\forall t$ (5g)

The objective is to find the optimal control inputs
$\mathbf{\hat{U}}_{T}=\\{\hat{u}_{0},..,\hat{u}_{T-1}\\}$, subject to the
constraints in Eqn. (5b)-(5g), over a finite planning horizon of length $T$
time-steps, which will result in the maximization of the coverage performance
objective $\mathcal{J}_{\text{coverage}}$ i.e., Eqn. (5a). As we have already
discussed in Sec. III-A, the constraints in Eqn. (5b)-(5c) are due to the
agent’s kinematic model. Then, the constraint in Eqn. (5d) implements the
coverage functionality. That is, for every point $p\in\mathcal{P}$ on the
object’s boundary $\partial\mathcal{C}$, it seeks to find a time-step
$\tau\in[1,..,T]$ inside the planning horizon, where the applied control
inputs $\\{\hat{u}_{0},..,\hat{u}_{\tau}\\}$ generate a trajectory (i.e., a
sequence of kinematic and FOV states) which result in the coverage of point
$p$ at time $\tau$ (i.e., $p\in\mathcal{F}_{\tau}$). Subsequently, the second
part of the constraint in Eqn. (5d) requires that at time-step $\tau$, at
which point $p$ resides inside the agent’s FOV, the same point $p$ must also
be visible i.e., overall Eqn. (5d) requires that point $p$ must reside within
the visible FOV of the agent. This constraint is required to account for the
fact that the agent’s FOV at some time-step $\tau\in[1,..,T]$ might be blocked
by obstacles, prohibiting the light-rays from reaching the camera’s optical
sensor and therefore rendering certain parts of the scene unobservable. In
essence the indicator function $\mathds{1}_{G_{\tau}}(p)$, utilizes the
function $G_{\tau}(p)$ in order to trace light-rays back to their source, and
thus determining the existence of a light-ray which traces-back to point $p$,
consequently determining $p$’s visibility on the agent’s camera. Thus
$\mathds{1}_{G_{\tau}}(p)$ is defined in this work as:

$\mathds{1}_{G_{\tau}}(p)=\begin{cases}1,&\quad\text{if
}G_{\tau}(p)\neq\emptyset,\\\ 0,&\quad\text{o.w. }\end{cases}$ (6)

Next, the constraint Eqn. (5e) implements a collision avoidance constraint
with the object of interest $\mathcal{C}$ i.e., the kinematic state $x_{t}$ of
the agent must always reside outside of the object’s boundary
$\partial\mathcal{C}$. Finally, the constraints in Eqn. (5f) and (5g) restrict
the agent’s state and control inputs within the desired operating bounds.
Next, we provide a more detailed discussion on how the formulation discussed
above can be tackled.

###  IV-A Visible FOV Determination

In this section we describe in more detail how the visible parts of the
agent’t FOV can be determined using ray-tracing. Essentially, we are
interested in determining whether the propagation of light-rays have been
blocked by obstacles and thus rendering parts of the scene unobservable. This
procedure is also useful in determining which of the points $p\in\mathcal{P}$
on the object’s boundary belong to the foreground and which points belong to
the background (and thus are unobservable from a specific agent state).

Recall that the set of all light-rays entering the camera’s optical sensor at
time $t$ is given by $\mathcal{R}_{x_{t},\theta_{t}}$, where $x_{t}$ and
$\theta_{t}$ is the agent’s kinematic state and FOV rotation angle at time $t$
respectively. Additionally, as we have already mentioned in Sec. III-C, the
boundary $\partial\mathcal{C}$ of the object of interest has been piece-wise
linearly approximated with a finite number of line segments
$\Delta\mathcal{C}=\\{L_{p}|p\in\mathcal{P}\\}$. We can now determine whether
point $p\in\mathcal{P}$ belongs to the visible FOV at time $t$, with respect
to the light-rays in $\mathcal{R}_{x_{t},\theta_{t}}$ by finding at least one
light-ray $\hat{R}\in\mathcal{R}_{x_{t},\theta_{t}}$ which traces back to
point $p$ as follows:

$G_{t}(p)=\\{\hat{R}\in\mathcal{R}_{x_{t},\theta_{t}}:\hat{R}\otimes\Delta\mathcal{C}=\hat{L}_{p}\\}$
(7)

Figure 3: The figure illustrates the ray-tracing procedure utilized in the
proposed coverage planning approach for identifying the visible parts of the
scene through the agent’s camera. In the above figure point $p^{\prime}$ is
visible, whereas point $p$ is occluded. The propagation of rays are from
$\rho\rightarrow x_{t}$.

where the operation $R\otimes\Delta\mathcal{C}$ defines the intersection of
the light-ray $R$, with all line segments $L_{p}\in\Delta\mathcal{C}$, and
returns the line segment $\hat{L}_{p}$ (if exists) which is the closest to the
camera’s optical sensor where the light-ray $R$ converges (i.e., see Fig.3).
Stated differently, $\hat{L}_{p}$ is the line-segment which light-ray $R$
intersects last. Since, light-rays propagate towards the camera’s optical
sensor where they converge, the line-segment encountered the last by a light-
ray is the line-segment $\hat{L}_{p}$ which belongs to the visible FOV. For
instance this is the line segment on the boundary of the object of interest
which can be imaged through the agent’s camera. On the other hand the line-
segment which first intersects with the light-ray (assuming the light-ray
intersects with more than one line-segments) is not visible e.g., this line-
segment might belong to the invisible side of the object as depicted in Fig.
3. Consequently, the point $p$ on the visible line-segment $\hat{L}_{p}$ is
also visible.

Let now the camera-ray $R$ given by $R=\\{\rho+s(x_{t}-\rho)|s\in[0,1]\\}$, to
have $x$ and $y$ cartesian coordinates given by $\rho(x)+s[x_{t}(x)-\rho(x)]$
and $\rho(y)+s[x_{t}(y)-\rho(y)]$ respectively. Also, let the $x$ and $y$
cartesian coordinates of the a line-segment
$L_{p}=\\{p+r(\hat{p}-p)|r\in[0,1],\hat{p}\neq p\\}$ to be
$p(x)+r[\hat{p}(x)-p(x)]$ and $p(y)+r[\hat{p}(y)-p(y)]$ respectively. The
intersection $R\otimes L_{p}$ of the light-ray $R$ with the line segment
$L_{p}$ can be computed by solving the following system of linear equations
for the two unknowns i.e., $(s,r)$:

$\begin{bmatrix}x_{t}(x)-\rho(x)&p(x)-\hat{p}(x)\\\
x_{t}(y)-\rho(y)&p(y)-\hat{p}(y)\end{bmatrix}\begin{bmatrix}s\\\
r\end{bmatrix}=\begin{bmatrix}p(x)-\rho(x)\\\ p(y)-\rho(y)\end{bmatrix}$ (8)

An intersection between $R$ and $L_{p}$ exists only when the solution
$(s,r)\in[0,1]$. It should be noted here that the event of no intersection
implies no visibility. Summarizing, the agent with kinematic state $x_{t}$,
and FOV state $\mathcal{F}_{t}$ (parameterized by angle $\theta_{t}$) observes
and covers point $p\in\mathcal{P}$ at time $t$ when: a) $p$ resides inside the
agent’s camera FOV $p\in\mathcal{F}_{t}$ and b) there exists a light-ray
$R\in\mathcal{R}_{x_{t},\theta_{t}}$ which traces back to point $p$ i.e., the
last intersection of $R$ is with the line-segment $L_{p}$ which contains point
$p$. This is exactly what the constraint in Eqn. (5d) describes.

###  IV-B Collision Avoidance Constraints

Collision avoidance between the agent and the object of interest is
implemented in this work by making sure that the agent’s kinematic state
$x_{t},\forall t$ always resides outside the convex hull of the object’s
$\mathcal{C}$ polygonal region approximation defined by its boundary
$\Delta\mathcal{C}$ as:

$\exists p\in\mathcal{P}:\alpha^{\top}_{p}x_{t}>\beta_{p}$ (9)

where $(a_{p},\beta_{p})$ are the coefficients (i.e., $a_{p}$ is the outward
normal vector on $L_{p}$, and $\beta_{p}$ is a constant) of the line equation
which contains the line segment $L_{p},p\in\mathcal{P}$. Since the area inside
the object of interest is modeled as a convex polygonal region with boundary
$\Delta\mathcal{C}$ defined by $|\mathcal{P}|$ linear equations containing the
line segments $L_{p}$, a collision occurs at time $t$ when
$\alpha^{\top}_{p}x_{t}\leq\beta_{p},\forall p$. Thus, a collision is avoided
if there exists $p\in\mathcal{P}$ such that
$\alpha^{\top}_{p}x_{t}>\beta_{p}$. This can be implemented with the following
set of constraints:

$\displaystyle\alpha^{\top}_{p}x_{t}>\beta_{p}-Mb^{p}_{t},~{}\forall t,p$
(10a)
$\displaystyle\sum_{p=1}^{|\mathcal{P}|}b^{p}_{t}\leq(|\mathcal{P}|-1),~{}\forall
t,$ (10b)

where the variable $b^{p}_{t}\in\\{0,1\\}$ is a binary decision variable, and
$M$ is a large positive constant. Essentially, Eqn. (10a) first checks if for
some time $t\in[1,..,T]$, the $p_{\text{th}}$ equation is true i.e.,
$\alpha^{\top}_{p}x_{t}>\beta_{p}$ and if so, it sets the binary decision
variable $b^{p}_{t}=0$. On the other hand, if
$\alpha^{\top}_{p}x_{t}\leq\beta_{p}$, then $b^{p}_{t}$ is activated in order
to make sure that Eqn. (10a) is satisfied. Subsequently, Eqn. (10b) requires
that the number of times $b^{p}_{t}$ is activated at each time $t$ is less or
equal to $(|\mathcal{P}|-1)$ which signifies that the agent’s kinematic state
resides outside the area enclosed by the object of interest i.e.,
$x_{t}\notin\mathcal{C},\forall t$.

###  IV-C Coverage Performance Objective

Finally the coverage performance objective i.e.,
$\mathcal{J}_{\text{coverage}}$ is defined in this work as:

$\mathcal{J}_{\text{coverage}}=\sum_{t=1}^{T}J(\mathcal{F}_{t},G_{t})\sigma(t)$
(11)

where $\mathcal{F}_{t}$ is the agent’s FOV state at time $t$ (associated with
the agent’s kinematic state $x_{t}$), $G_{t}$ is the ray-tracing function, and
$J(\mathcal{F}_{t},G_{t})$ is defined as:

$J(\mathcal{F}_{t},G_{t})=\sum_{p\in\mathcal{P}}\ell(\mathcal{F}_{t},G_{t},p)$
(12)

The function $\ell(\mathcal{F}_{\tau},G_{\tau},p)\in\\{0,1\\}$ is further
given by:

$\ell(\mathcal{F}_{t},G_{t},p)=\begin{cases}1,&\exists\tau\leq
T:\left(p\in\mathcal{F}_{\tau}\right)\wedge\mathds{1}_{G_{\tau}}(p),\\\
0,&\text{o.w. }\end{cases}$ (13)

In essence $\mathcal{J}_{\text{coverage}}$ is maximized when all points
$p\in\mathcal{P}$ become visible at some point in time inside the planning
horizon and they are covered by the agent’s camera. Finally, the time-
dependent term $\sigma(t)$ is used here to penalize points that are covered
later in the planning horizon e.g., $\sigma(t)=(T-t)T^{-1}$, thus encouraging
the agent to finish the coverage mission as soon as possible.

##  V RL-based Integrated Ray-tracing and Coverage Control

The optimal control problem i.e., Problem (P1), presented in the previous
section is quite challenging to be solved efficiently in real-time, mainly due
to the non-linear and non-convex constraints introduced when someone attempts
to solve Eqn. (8), which is required by the constraint in Eqn. (5d). It should
also be noted here that the computational complexity of the OCP in (P1) grows
at best linearly with the size of the planning horizon; which can introduce
additional challenges for large problems. Finally, the size of the planning
horizon depends on the coverage scenario (e.g., the size of the object of
interest, the initial location of the agent, etc.), and thus cannot be easily
determined and tuned a-priori. If the planning horizon is chosen too short the
problem becomes infeasible and the coverage mission fails. On the other hand,
if the size of the planning horizon is overestimated, unnecessary
computational complexity will be introduced. In order to tackle some of the
challenges mentioned above, the optimal control problem in (P1) is
approximated and transformed into a finite Markov decision process (MDP) for
which an optimal control policy can be learned from experience using
reinforcement learning (RL).

###  V-A Problem re-formulation using MDP

The problem of integrated ray-tracing and coverage planning control is re-
formulated in this section as a Markov decision process, which is then solved
using reinforcement learning. A finite Markov Decision Process (MDP) [33] is a
discrete-time stochastic control process which can be defined with a set of
components i.e., a tuple
$\langle\mathcal{S},\mathcal{A},\mathcal{T},R_{w}\rangle$; where $\mathcal{S}$
represents the state-space of the system (i.e., a finite set describing all
possible states of the system), $\mathcal{A}$ is the action space (i.e., the
finite set of actions available to the agent), $\mathcal{T}$ is a discrete-
time state transition function i.e.,
$\mathcal{T}:\mathcal{S}\times\mathcal{A}\mapsto\mathcal{S}$, which describes
the evolution of the system due to actions, and finally
$R_{w}:\mathcal{S}\times\mathcal{A}\mapsto\mathbb{R}$ is the reward function
which returns the immediate reward when action $a$ is executed at state $s$
and the system transitions to some new state $s^{\prime}$. Reinforcement
learning [34] can be used to solve an MDP by finding a control policy (i.e., a
mapping from states to actions) which maximizes the cumulative reward. In
order to do that, RL methods track the interactions of the agent with the
environment, and learn to optimize the control policy using reinforcement,
i.e., a feedback in the form of reward or punishment. RL methods are
particularly useful for solving a) large MDP problems, which are otherwise
intractable with exact dynamic programming methods, and b) problems where the
exact mathematical model of the MDP is not known. We can now go ahead and
discuss how we have transformed our problem to an MDP. More specifically, the
MDP components i.e., $\langle\mathcal{S},\mathcal{A},\mathcal{T},R_{w}\rangle$
are constructed for our problem as follows:

####  V-A1 State-space $\mathcal{S}$

In order to construct the system’s state-space we first decompose the
environment $\mathcal{W}$ into a finite number of non-overlapping cells
$\\{c_{1},..,c_{|\hat{\mathcal{W}}|}\\}$, essentially creating a uniform grid
$\hat{\mathcal{W}}=\bigcup_{i=1}^{|\hat{\mathcal{W}}|}c_{i}$, which within the
agent can evolve by moving from one cell to another. Let $\bar{p}$ to denote
the center of mass (i.e., centroid) of the object of interest $\mathcal{C}$,
given by $\bar{p}=|\mathcal{P}|^{-1}\sum_{i=1}^{|\mathcal{P}|}p_{i}$. Let
$\hat{d}_{t}=q(\left\lVert x_{t}-\bar{p}\right\rVert_{2})$ to denote the
quantized Euclidian distance between the agent’s kinematic state and the
object’s center of mass $\bar{p}$, where the quantizer function $q()$ maps the
distance values into the finite set $\hat{D}$. Subsequently the MDP state-
space can now be defined as:

$\mathcal{S}=\hat{\mathcal{W}}\times\hat{P}\times\hat{D}$ (14)

where $\hat{P}=[0,..,|\mathcal{P}|]$, with $\mathcal{P}$ denoting the set of
points sampled from the object’s boundary. Therefore, a particular state
$s_{t}=[s_{a},s_{b},s_{c}]\in\mathcal{S}$ at time $t$ describes the agent’s
kinematic state i.e., position $s_{a}\in\tilde{\mathcal{W}}$, the total number
of points $s_{b}\in\hat{P}$ that are observed and covered for the first time
from state $s_{a}$, and the agent’s distance $s_{c}\in\hat{D}$ from the
object’s centroid.

####  V-A2 Action-space $\mathcal{A}$

The agent’s kinematic model and control inputs are discretized in space
according to:

$\hat{x}_{t}=\hat{x}_{t-1}+l_{R}\Delta_{R}\begin{bmatrix}\cos(l_{\vartheta}\Delta_{\vartheta})\\\
\sin(l_{\vartheta}\Delta_{\vartheta})\end{bmatrix},\begin{array}[]{l}l_{\vartheta}=0,...,N_{\vartheta}\\\
l_{R}=0,...,N_{R}\end{array}$ (15)

where the agent’s kinematic control input has now become
$u^{d}_{t}=[l_{R}\Delta_{R},l_{\vartheta}\Delta_{\vartheta}]$, where
$\Delta_{R}$ is the radial step size, $\Delta_{\vartheta}=2\pi/N_{\vartheta}$,
and the parameters $(N_{\vartheta},N_{R})$ are used to specify the total
number of admissible control actions, denoted as $\mathcal{U}^{d}$.
Specifically, the parameters $(\Delta_{R},N_{\vartheta},N_{R})$ are used to
ultimately determine the evolution of the agent’s kinematic state
$\hat{x}_{t}$ inside the discretized representation of the world. The camera’s
FOV rotation signal is also discretized to take its values from within a
finite set of rotation angles $\theta^{d}_{t}\in\Theta^{d}$. Consequently, the
camera FOV state $\mathcal{F}_{t}$ at time $t$, also takes its values from
within a finite set of possible FOV configurations. The agent’s control
signal, referred to as the agent’s action hereafter, is denoted as
$a_{t}=[u^{d}_{t},\theta^{d}_{t}]\in\mathcal{U}^{d}\times\Theta^{d}$

####  V-A3 Reward $R_{w}$

Finally, we have designed a reward function $R_{w}(s_{t},a_{t})$ which closely
resembles the constraints and objectives of the optimal control problem
presented in Sec. IV, and which is given by:

$R_{w}(s_{t},a_{t})=-w_{1}-w_{2}\mathds{1}_{s_{t+1}\in\mathcal{C}}+w_{3}J(s_{t+1})$
(16)

where the notation $(s_{t},a_{t},s_{t+1})$ refers to current state of the
agent $s_{t}$, the action $a_{t}$ applied at state $s_{t}$, and the transition
of the agent to the new state $s_{t+1}$ in the next time-step. The parameters
$w_{1},w_{2}$ and $w_{3}$ are tuning weights which control the system’s
reinforcement behavior. More specifically, $w_{1}$ is a punishment term which
used in Eqn. (16) in order to penalize the agent for every time-step that
passes and the object of interest remains not fully covered. In other words
this term encourages the agent to finish the coverage mission as soon as
possible. Next, the indicator function $\mathds{1}_{s_{t+1}\in\mathcal{C}}$,
which is weighted by $w_{2}$, uses the procedure described in Sec. IV-B to
identify whether the state-action pair $(s_{t},a_{t})$ results in a collision
between the agent and the object of interest at the next time-step, and if so
the agent is penalized with $-w_{2}$. Finally, the last term rewards the agent
for its successful coverage efforts. Specifically, in the last term we make
use of the function $J()$ described in Eqn. (12) to count the number of points
which have been covered for the first time when the agent transitions to the
new state $s_{t+1}$; and we generate a reward signal weighted by $w_{3}$ as
shown.

###  V-B Off-policy Temporal-Difference Control

Based on the MDP formulation discussed above the problem of integrated ray-
tracing and coverage planning control, shown in (P1), can now be written as:

$Q_{\pi}(s,a)=\mathbb{E}_{\pi}\left[\sum_{t=0}^{\infty}\gamma^{t}R_{w}(s_{t},a_{t})|s_{0}=s,a_{0}=a\right]$
(17)

where $\gamma\in[0,1]$ is a discount factor which determines the importance of
future rewards. The state-action value function $Q_{\pi}(s,a)$ measures the
expected infinite horizon discounted cumulative reward when the agent starts
from state $s$, takes action $a$ and then follows the control policy $\pi$.
Subsequently, the optimal state-action value function is obtained as
$Q_{\star}(s,a)=\max_{\pi}Q_{\pi}(s,a)$ by following an optimal control policy
$\pi_{\star}$ which always chooses the optimal action in every state. The
problem re-formulation discussed so far is the MDP equivalent of the
maximization of the objective in Eqn. (5a) under the constraints (5b)-(5g),
i.e., Problem (P1). In this work, we use off-policy temporal-difference
control or otherwise known as Q-learning [35], in order to learn an optimal
control policy $\pi_{\star}$ from experience. In particular, Q-learning is an
iterative bootstrapping process which uses temporal-differencing in order to
learn the state-action value function as: $Q(s_{t},a_{t})=$

$Q(s_{t},a_{t})+\alpha\left[R_{w}(s_{t},a_{t})+\gamma\max_{a}Q(s_{t+1},a)-Q(s_{t},a_{t})\right]$

where $\alpha$ is the learning rate. Being an off-policy method, Q-learning is
able to learn the state-action value function $Q(s,a)$, which directly
approximates $Q_{\star}(s,a)$, i.e., the optimal state-action value function,
independent of the policy being followed. When the optimal state-action value
function $Q_{\star}(s,a)$ is obtained, the optimal policy $\pi_{\star}$ can be
derived by taking at each time-step the actions that maximize
$Q_{\star}(s,a)$. Finally, Q-learning guarantees that an optimal policy can be
found for any finite MDP given enough exploration time and a partly-random
policy [36].

##  VI Evaluation

###  VI-A Simulation Setup

The simulation setup used for the evaluation of the proposed approach is as
follows: The surveillance area $\mathcal{W}$ of size $20$m by $20$m is
converted into a uniform grid $\hat{\mathcal{W}}$ of size 10 by 10. At any
time-step $t$, the quantized distance $\hat{d}_{t}$ between the agent and the
centroid of the object of interest takes values from within the set
$\hat{D}=\\{0,0.5,1,1.5,...,\lceil\sqrt{2}\times 20\rceil\\}$m. The agent
kinematics follow Eqn. (15) with $\Delta_{R}=2$m, $N_{\vartheta}=8$, and
$N_{R}=1$. The agent’s camera angle of view angle $\varphi$ is set to 40deg,
and the sensing range $h=10$m. In total we consider 5 camera rotation angles
$\theta^{d}$ which take values in the set
$\Theta^{d}=\\{-85,-42.5,0,42.5,85\\}$deg as shown in Fig. 1. At each time-
step $t$ and for each camera FOV configuration 5 light-rays enter the optical
sensor i.e., $|\mathcal{R}_{\hat{x}_{t},\theta^{d}_{t}}|=5$. The region/object
of interest $\mathcal{C}$ is represented by a bell-shaped curve (as
illustrated in Fig. 5), given by
$f(x)=a\times\text{exp}\left(\frac{-(x-b)^{2}}{2c^{2}}\right)$ where $a,b$ and
$c$ are free parameters. The objective is to to cover a total of 11 points
$\mathcal{P}=\\{p_{1},..,p_{11}\\}$ uniformly sampled from the object’s
boundary $\partial\mathcal{C}$. The punishment and reward parameters
$w_{1},w_{2}$ and $w_{3}$ in Eqn. (16) are set as 1, 100, and 2 respectively.
The Q-learning parameters $\alpha$ and $\gamma$ (i.e., learning-rate and
discount factor) are set to 0.1 and 0.8 respectively. We use a decreasing
$\epsilon$-greedy exploration strategy with an initial $\epsilon=0.9$ which
decays in every time-step as $\epsilon_{t+1}=0.9999\epsilon_{t}$. Finally, for
each episode the agent is randomly initialized inside $\hat{\mathcal{W}}$, and
interacts with the environment until the terminal condition is met i.e., all
points $|\mathcal{P}|$ have been covered or the time limit of 100 steps has
been reached.

###  VI-B Performance Evaluation

Figure 4: The figure shows the cumulative reward obtained in each episode
during training. (a) The episodic reward obtained during training with a
single configuration of the object of interest, and (b) the corresponding
coverage time. (c) The evolution of the episodic reward during the training
with multiple random configurations of the object of interest. Figure 5: The
figure shows an illustrative example of four different random realizations of
the object of interest that have been sampled during training. Figure 6: The
figure illustrates 3 different scenarios (a), (b) and (c) for the problem of
integrated ray-tracing and coverage planning control, where the agent follows
the learned optimal control law $\pi_{\star}$.

We begin our evaluation by investigating the learning performance of the
proposed approach for a single and multiple objects of interest (i.e.,
different realizations from the same object class). Specifically, Fig. 4(a)
shows the infinite horizon discounted cumulative reward per episode when the
proposed system is trained over 5000 episodes for a single object of interest
with parameters $(a,b,c)=(8,8,2)$. As the figure shows the episodic reward
increases over time as the agent learns the optimal policy which maximizes the
state-action value function, in which case it reaches a plateau as shown in
the figure. Since the total cumulative reward obtained in each episode depends
on the initial agent position, and because the agent in each episode is
randomly initialized, the optimal control policy results in a cumulative
reward which takes its values within a small range of values as shown in Fig.
4(a) after episode 2000. This is also evident in Fig. 4(b) which shows the
total coverage time during training. As we can observe once the optimal policy
is identified i.e., after 2000 episodes the agent requires on average 12 time-
steps to finish the mission (depending on its initial state), whereas
initially during the early episodes the mission is not finished not even
within 100 time-steps. Then Fig. 4(c) shows the learning curve also in terms
of the discounted cumulative reward when the agent is trained over 100000
episodes for multiple realizations (i.e., parameter configurations) of the
object of interest whose parameters $a,b$, and $c$ are uniformly sampled
within the intervals $[1,18]$, $[5,15]$ and $[1,4]$ respectively. In this
scenario, for each episode a new realization of the object of interest is
sampled and the agent is randomly initialized inside the surveillance area. An
illustrative example of different object realizations sampled in each episode
during training is shown in Fig. 5. Figure 4(c) indicates that the optimal
policy can be identified and the state-action value function can be learned
for multiple realizations of the object of interest. Figure 6 shows in more
detail the agent’s coverage plan under the learned policy for 3 random
realizations of the object of interest. More specifically, Fig. 6 depicts 3
different parameter configurations of the object of interest i.e., with the
parameters $(a,b,c)$ set as $(10,10,2.5)$, $(5,6,1.4)$, and $(14,15,1)$ for
the objects depiceted in Fig. 6(a), Fig. 6(b), and Fig. 6(c) respectively. In
Fig. 6 the first row illustrates the agent’s kinematic states ($x_{t}$) under
the learned control policy. The second row in addition to the agent’s
kinematic states, also shows the camera’s FOV states $\mathcal{F}_{t}$ under
the learned policy. The points $p_{1},...,p_{11}$ on the object’s boundary to
be covered are shown as $\circ$. The agent’s kinematic states $x_{t}$ (denoted
with $\square$) are colored-coded and numbered based on their time-index $t$
as shown in the figure (i.e., from black to deep red, indicating start and end
times respectively).

In addition, when the agent with kinematic state $x_{t}$ covers point $p$ at
time $t$, we mark point $p$ with the same color as $x_{t}$, so that we can
distiguish at which agent states each point is observed and covered. For
instance in Fig. 6(a), the agent at time-step $t=1$, with state $x_{1}$
(colored black) covers points $p_{2}$ and $p_{3}$, and therefore these points
are colored black as well. At the next time-step the agent is at $x_{2}$
(colored blue) and observes point $p_{1}$ (note that $p_{1}$ resides inside
the agent’s camera FOV), thus $p_{1}$ is marked blue as shown. At $t=3$ the
agent’s state is colored cyan and so does point $p_{4}$ which is observed at
that time, and so on and so forth. In order to make the figure easier to read,
we only display the FOV configuration at time-steps for which a coverage event
occurs e.g, in Fig. 6(b) for instance the FOV configuration at time-steps
$t=[1,2,5,6]$ is ommitted for clarity.

Note that the control policy that the agent follows integrates ray-tracing.
This is evident in Fig.6(a) where at time-step $t=4$ the agent is at state
$x_{4}$ (colored green) and points $p_{5}$ and $p_{7}$ reside inside the
agent’s triangular FOV. However, the agent only observes point $p_{5}$ as
indicated by its color-coding i.e., point $p_{5}$ is colored green, whereas
point $p_{7}$ is marked yellow and is not observed at $t=4$ since it is not
visible from the agent’s state $x_{4}$. Observe that the light-rays at $x_{4}$
are blocked by the object, and cannot be traced back to $p_{7}$. The point
$p_{7}$ becomes visible and is covered when the agent moves at state $x_{6}$.
This behavior can also be observed in Fig.6(b) when the agent is at $x_{3}$,
in Fig. 6(c) when the agent is at state $x_{5}$, etc. Finally, note that the
agent follows an optimal control policy which generalizes for multiple
realizations of the object of interest, i.e., in all scenarios, all points are
covered, the coverage time is minimized, and collisions between the agent and
the object of interest are avoided. It is also worth noting that the optimal
control policy has learned the optimal planning horizon length i.e., the
coverage mission finishes in 9, 8 and 12 times-steps for the scenarios shown
in Fig. 6(a), Fig. 6(b), and Fig. 6(c) respectively, as opposed to the OCP
formulation in which the optimal horizon length must be fined-tuned a-priori.

##  VII Conclusion

In this work we have presented an integrated ray-tracing and coverage planning
control approach which enables a mobile agent to jointly decide its kinematic
and camera control inputs in such a way so that the total surface area of an
object of interest is covered in the minimum amount of time. The proposed
approach integrates ray-tracing into the coverage planning problem, in order
to trace the propagation of light-rays back to their source thus ultimately
identifying which parts of the scene are visible through the agent’s camera.
The problem is posed as a Markov decision process (MDP) and an optimal control
law is obtained using reinforcement learning i.e., Q-learning. Future works
include the investigation of this problem in 3D environments, its extension to
multiple learning agents, and the adaptation of the proposed approach to
continuous state-action spaces using deep reinforcement learning.

## Acknowledgments

This work is supported by the European Union’s Horizon 2020 research and
innovation programme under grant agreement No 739551 (KIOS CoE), by the
European Union Civil Protection under grant agreement No 873240 (AIDERS), and
from the Republic of Cyprus through the Directorate General for European
Programmes, Coordination and Development.

## References

  * [1] S. Papaioannou, P. Kolios, T. Theocharides, C. G. Panayiotou, and M. M. Polycarpou, “Decentralized search and track with multiple autonomous agents,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_. IEEE, 2019, pp. 909–915.
  * [2] C. Liu and J. K. Hedrick, “Model predictive control-based target search and tracking using autonomous mobile robot with limited sensing domain,” in _2017 American Control Conference (ACC)_. IEEE, 2017, pp. 2937–2942.
  * [3] S. Papaioannou, P. Kolios, T. Theocharides, C. G. Panayiotou, and M. M. Polycarpou, “A cooperative multiagent probabilistic framework for search and track missions,” _IEEE Transactions on Control of Network Systems_ , vol. 8, no. 2, pp. 847–858, 2021.
  * [4] S. Papaioannou, P. Kolios, C. G. Panayiotou, and M. M. Polycarpou, “Cooperative simultaneous tracking and jamming for disabling a rogue drone,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, pp. 7919–7926.
  * [5] P. Valianti, S. Papaioannou, P. Kolios, and G. Ellinas, “Multi-agent coordinated close-in jamming for disabling a rogue drone,” _IEEE Transactions on Mobile Computing_ , 2021, doi:10.1109/TMC.2021.3062225.
  * [6] S. Papaioannou, P. Kolios, and G. Ellinas, “Downing a rogue drone with a team of aerial radio signal jammers,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, pp. 2555–2562.
  * [7] M. Corah and N. Michael, “Distributed submodular maximization on partition matroids for planning on large sensor networks,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 6792–6799.
  * [8] N. Zhou, C. G. Cassandras, X. Yu, and S. B. Andersson, “Decentralized event-driven algorithms for multi-agent persistent monitoring tasks,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_ , 2017, pp. 4064–4069.
  * [9] J. Moon, S. Papaioannou, C. Laoudias, P. Kolios, and S. Kim, “Deep reinforcement learning multi-uav trajectory control for target tracking,” _IEEE Internet of Things Journal_ , vol. 8, no. 20, pp. 15 441–15 455, 2021\.
  * [10] E. Galceran and M. Carreras, “A survey on coverage path planning for robotics,” _Robotics and Autonomous Systems_ , vol. 61, no. 12, pp. 1258–1276, 2013.
  * [11] H. Choset and P. Pignon, “Coverage path planning: The boustrophedon cellular decomposition,” in _Field and Service Robotics_. Springer, 1998, pp. 203–209.
  * [12] E. U. Acar, H. Choset, A. A. Rizzi, P. N. Atkar, and D. Hull, “Morse decompositions for coverage tasks,” _The International Journal of Robotics Research_ , vol. 21, no. 4, pp. 331–344, 2002.
  * [13] E. Acar, H. Choset, and J. Y. Lee, “Sensor-based coverage with extended range detectors,” _IEEE Transactions on Robotics_ , vol. 22, no. 1, pp. 189–198, 2006.
  * [14] W. H. Huang, “Optimal line-sweep-based decompositions for coverage algorithms,” in _Proceedings 2001 ICRA. IEEE International Conference on Robotics and Automation_ , vol. 1. IEEE, 2001, pp. 27–32.
  * [15] R. Mannadiar and I. Rekleitis, “Optimal coverage of a known arbitrary environment,” in _2010 IEEE International Conference on Robotics and Automation_. IEEE, 2010, pp. 5525–5530.
  * [16] M. N. A. Wahab, S. Nefti-Meziani, and A. Atyabi, “A comparative review on mobile robot path planning: Classical or meta-heuristic methods?” _Annual Reviews in Control_ , vol. 50, pp. 233–252, 2020.
  * [17] W. Li and C. Cassandras, “Distributed cooperative coverage control of sensor networks,” in _Proceedings of the 44th IEEE Conference on Decision and Control_ , 2005, pp. 2542–2547.
  * [18] A. Janchiv, D. Batsaikhan, B. Kim, W. G. Lee, and S.-G. Lee, “Time-efficient and complete coverage path planning based on flow networks for multi-robots,” _International Journal of Control, Automation and Systems_ , vol. 11, no. 2, pp. 369–376, 2013.
  * [19] N. Karapetyan, K. Benson, C. McKinney, P. Taslakian, and I. Rekleitis, “Efficient multi-robot coverage of a known environment,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 1846–1852.
  * [20] G. Cannata and A. Sgorbissa, “A minimalist algorithm for multirobot continuous coverage,” _IEEE Transactions on Robotics_ , vol. 27, no. 2, pp. 297–312, 2011.
  * [21] M. Schwager, D. Rus, and J.-J. Slotine, “Decentralized, adaptive coverage control for networked robots,” _The International Journal of Robotics Research_ , vol. 28, no. 3, pp. 357–375, 2009.
  * [22] J. Cortes, S. Martinez, T. Karatas, and F. Bullo, “Coverage control for mobile sensing networks,” _IEEE Transactions on Robotics and Automation_ , vol. 20, no. 2, pp. 243–255, 2004.
  * [23] L. C. Pimenta, V. Kumar, R. C. Mesquita, and G. A. Pereira, “Sensing and coverage for a network of heterogeneous robots,” in _47th IEEE Conference on Decision and Control_. IEEE, 2008, pp. 3947–3952.
  * [24] Y. Li, H. Chen, M. J. Er, and X. Wang, “Coverage path planning for UAVs based on enhanced exact cellular decomposition method,” _Mechatronics_ , vol. 21, no. 5, pp. 876–885, 2011.
  * [25] A. Xu, C. Viriyasuthee, and I. Rekleitis, “Efficient complete coverage of a known arbitrary environment with applications to aerial operations,” _Autonomous Robots_ , vol. 36, no. 4, pp. 365–381, 2014.
  * [26] L. Paull, C. Thibault, A. Nagaty, M. Seto, and H. Li, “Sensor-driven area coverage for an autonomous fixed-wing unmanned aerial vehicle,” _IEEE Transactions on Cybernetics_ , vol. 44, no. 9, pp. 1605–1618, 2014.
  * [27] S. Papaioannou, P. Kolios, T. Theocharides, C. G. Panayiotou, and M. M. Polycarpou, “Integrated guidance and gimbal control for coverage planning with visibility constraints,” _IEEE Transactions on Aerospace and Electronic Systems_ , 2022, doi:10.1109/TAES.2022.3199196.
  * [28] ——, “Towards Automated 3D Search Planning for Emergency Response Missions,” _Journal of Intelligent & Robotic Systems_, vol. 103, no. 1, pp. 1–19, 2021.
  * [29] ——, “3D Trajectory Planning for UAV-based Search Missions: An Integrated Assessment and Search Planning Approach,” in _2021 International Conference on Unmanned Aircraft Systems (ICUAS)_. IEEE, 2021, pp. 517–526.
  * [30] ——, “Uav-based receding horizon control for 3d inspection planning,” in _2022 International Conference on Unmanned Aircraft Systems (ICUAS)_ , 2022, pp. 1121–1130.
  * [31] M. Theile, H. Bayerlein, R. Nai, D. Gesbert, and M. Caccamo, “UAV coverage path planning under varying power constraints using deep reinforcement learning,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 1444–1449.
  * [32] G. Sanna, S. Godio, and G. Guglieri, “Neural Network Based Algorithm for Multi-UAV Coverage Path Planning,” in _2021 International Conference on Unmanned Aircraft Systems (ICUAS)_. IEEE, 2021, pp. 1210–1217.
  * [33] M. L. Puterman, “Markov decision processes,” _Handbooks in operations research and management science_ , vol. 2, pp. 331–434, 1990.
  * [34] R. S. Sutton, “A special issue of machine learning on reinforcement learning,” _Machine learning_ , vol. 8, 1992.
  * [35] R. S. Sutton and A. G. Barto, “Reinforcement learning: An introduction,” _Robotica_ , vol. 17, no. 2, pp. 229–235, 1999.
  * [36] C. J. Watkins and P. Dayan, “Q-learning,” _Machine learning_ , vol. 8, no. 3, pp. 279–292, 1992.


* [82] W. Porod, F. Staub and A. Vicente, _A Flavor Kit for BSM models_ , _Eur. Phys. J. C_ 74 (2014) 2992 [1405.1434].
  * [83] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _MicrOMEGAs: A Program for calculating the relic density in the MSSM_ , _Comput. Phys. Commun._ 149 (2002) 103 [hep-ph/0112278].
  * [84] G. Belanger, F. Boudjema, C. Hugonie, A. Pukhov and A. Semenov, _Relic density of dark matter in the NMSSM_ , _JCAP_ 09 (2005) 001 [hep-ph/0505142].
  * [85] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _MicrOMEGAs 2.0: A Program to calculate the relic density of dark matter in a generic model_ , _Comput. Phys. Commun._ 176 (2007) 367 [hep-ph/0607059].
  * [86] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _micrOMEGAs: A Tool for dark matter studies_ , _Nuovo Cim. C_ 033N2 (2010) 111 [1005.4133].
  * [87] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _micrOMEGAs $\\_$3: A program for calculating dark matter observables_, _Comput. Phys. Commun._ 185 (2014) 960 [1305.0237].
  * [88] D. Barducci, G. Belanger, J. Bernon, F. Boudjema, J. Da Silva, S. Kraml et al., _Collider limits on new physics within micrOMEGAs $\\_$4.3_, _Comput. Phys. Commun._ 222 (2018) 327 [1606.03834].
  * [89] A. Fowlie and M. H. Bardsley, _Superplot: a graphical interface for plotting and analysing MultiNest output_ , _Eur. Phys. J. Plus_ 131 (2016) 391 [1603.00555].
  * [90] S. Matsumoto, S. Mukhopadhyay and Y.-L. S. Tsai, _Effective Theory of WIMP Dark Matter supplemented by Simplified Models: Singlet-like Majorana fermion case_ , _Phys. Rev. D_ 94 (2016) 065034 [1604.02230].
  * [91] P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak and G. Weiglein, _$HiggsSignals$ : Confronting arbitrary Higgs sectors with measurements at the Tevatron and the LHC_, _Eur. Phys. J. C_ 74 (2014) 2711 [1305.1933].
  * [92] O. Stål and T. Stefaniak, _Constraining extended Higgs sectors with HiggsSignals_ , _PoS_ EPS-HEP2013 (2013) 314 [1310.4039].
  * [93] P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak and G. Weiglein, _Probing the Standard Model with Higgs signal rates from the Tevatron, the LHC and a future ILC_ , _JHEP_ 11 (2014) 039 [1403.1582].
  * [94] P. Bechtle, S. Heinemeyer, T. Klingl, T. Stefaniak, G. Weiglein and J. Wittbrodt, _HiggsSignals-2: Probing new physics with precision Higgs measurements in the LHC 13 TeV era_ , _Eur. Phys. J. C_ 81 (2021) 145 [2012.09197].
  * [95] P. Bechtle, O. Brein, S. Heinemeyer, G. Weiglein and K. E. Williams, _HiggsBounds: Confronting Arbitrary Higgs Sectors with Exclusion Bounds from LEP and the Tevatron_ , _Comput. Phys. Commun._ 181 (2010) 138 [0811.4169].
  * [96] P. Bechtle, O. Brein, S. Heinemeyer, G. Weiglein and K. E. Williams, _HiggsBounds 2.0.0: Confronting Neutral and Charged Higgs Sector Predictions with Exclusion Bounds from LEP and the Tevatron_ , _Comput. Phys. Commun._ 182 (2011) 2605 [1102.1898].
  * [97] P. Bechtle, O. Brein, S. Heinemeyer, O. Stal, T. Stefaniak, G. Weiglein et al., _Recent Developments in HiggsBounds and a Preview of HiggsSignals_ , _PoS_ CHARGED2012 (2012) 024 [1301.2345].
  * [98] P. Bechtle, O. Brein, S. Heinemeyer, O. Stål, T. Stefaniak, G. Weiglein et al., _$\mathsf{HiggsBounds}-4$ : Improved Tests of Extended Higgs Sectors against Exclusion Bounds from LEP, the Tevatron and the LHC_, _Eur. Phys. J. C_ 74 (2014) 2693 [1311.0055].
  * [99] P. Bechtle, D. Dercks, S. Heinemeyer, T. Klingl, T. Stefaniak, G. Weiglein et al., _HiggsBounds-5: Testing Higgs Sectors in the LHC 13 TeV Era_ , _Eur. Phys. J. C_ 80 (2020) 1211 [2006.06007].
  * [100] L. M. Carpenter, R. Colburn, J. Goodman and T. Linden, _Indirect Detection Constraints on s and t Channel Simplified Models of Dark Matter_ , _Phys. Rev. D_ 94 (2016) 055027 [1606.04138].
  * [101] X.-J. Huang, C.-C. Wei, Y.-L. Wu, W.-H. Zhang and Y.-F. Zhou, _Antiprotons from dark matter annihilation through light mediators and a possible excess in AMS-02 $\bar{p}/p$ data_, _Phys. Rev. D_ 95 (2017) 063021 [1611.01983].
  * [102] Particle Data Group collaboration, M. Tanabashi et al., _Review of Particle Physics_ , _Phys. Rev. D_ 98 (2018) 030001.
  * [103] B. O’Leary and J. E. Camargo-Molina, “VevaciousPlusPlus.” https://github.com/JoseEliel/VevaciousPlusPlus, 2014.
  * [104] J. E. Camargo-Molina, B. O’Leary, W. Porod and F. Staub, _$\mathbf{Vevacious}$ : A Tool For Finding The Global Minima Of One-Loop Effective Potentials With Many Scalars_, _Eur. Phys. J. C_ 73 (2013) 2588 [1307.1477].
  * [105] CMS collaboration, A. M. Sirunyan et al., _Searches for pair production of charginos and top squarks in final states with two oppositely charged leptons in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _JHEP_ 11 (2018) 079 [1807.07799].
  * [106] CMS collaboration, A. M. Sirunyan et al., _Search for supersymmetric partners of electrons and muons in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _Phys. Lett. B_ 790 (2019) 140 [1806.05264].
  * [107] CMS collaboration, A. M. Sirunyan et al., _Combined search for electroweak production of charginos and neutralinos in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _JHEP_ 03 (2018) 160 [1801.03957].
  * [108] CMS collaboration, A. M. Sirunyan et al., _Search for supersymmetry with Higgs boson to diphoton decays using the razor variables at $\sqrt{s}=$ 13 TeV_, _Phys. Lett. B_ 779 (2018) 166 [1709.00384].
  * [109] CMS collaboration, A. M. Sirunyan et al., _Search for electroweak production of charginos and neutralinos in multilepton final states in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _JHEP_ 03 (2018) 166 [1709.05406].
  * [110] CMS collaboration, A. M. Sirunyan et al., _Search for new phenomena in final states with two opposite-charge, same-flavor leptons, jets, and missing transverse momentum in pp collisions at $\sqrt{s}=13$ TeV_, _JHEP_ 03 (2018) 076 [1709.08908].
  * [111] ATLAS collaboration, M. Aaboud et al., _Search for electroweak production of supersymmetric particles in final states with two or three leptons at $\sqrt{s}=13\,$TeV with the ATLAS detector_, _Eur. Phys. J. C_ 78 (2018) 995 [1803.02762].
  * [112] ATLAS collaboration, M. Aaboud et al., _Search for chargino and neutralino production in final states with a Higgs boson and missing transverse momentum at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Phys. Rev. D_ 100 (2019) 012006 [1812.09432].
  * [113] ATLAS collaboration, M. Aaboud et al., _Search for chargino-neutralino production using recursive jigsaw reconstruction in final states with two or three charged leptons in proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Phys. Rev. D_ 98 (2018) 092012 [1806.02293].
  * [114] ATLAS collaboration, G. Aad et al., _Search for chargino-neutralino production with mass splittings near the electroweak scale in three-lepton final states in $\sqrt{s}$=13 TeV $pp$ collisions with the ATLAS detector_, _Phys. Rev. D_ 101 (2020) 072001 [1912.08479].
  * [115] ATLAS collaboration, G. Aad et al., _Search for electroweak production of charginos and sleptons decaying into final states with two leptons and missing transverse momentum in $\sqrt{s}=13$ TeV $pp$ collisions using the ATLAS detector_, _Eur. Phys. J. C_ 80 (2020) 123 [1908.08215].
  * [116] ATLAS collaboration, G. Aad et al., _Search for direct production of electroweakinos in final states with one lepton, missing transverse momentum and a Higgs boson decaying into two $b$-jets in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Eur. Phys. J. C_ 80 (2020) 691 [1909.09226].
  * [117] CMS collaboration, A. M. Sirunyan et al., _Search for supersymmetry in final states with two oppositely charged same-flavor leptons and missing transverse momentum in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _JHEP_ 04 (2021) 123 [2012.08600].
  * [118] ATLAS collaboration, G. Aad et al., _Search for chargino–neutralino pair production in final states with three leptons and missing transverse momentum in $\sqrt{s}=13$ TeV pp collisions with the ATLAS detector_, _Eur. Phys. J. C_ 81 (2021) 1118 [2106.01676].
  * [119] ATLAS collaboration, G. Aad et al., _Search for direct production of electroweakinos in final states with one lepton, missing transverse momentum and a Higgs boson decaying into two $b$-jets in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Eur. Phys. J. C_ 80 (2020) 691 [1909.09226].
  * [120] ATLAS collaboration, G. Aad et al., _Search for charginos and neutralinos in final states with two boosted hadronically decaying bosons and missing transverse momentum in $pp$ collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector_, _Phys. Rev. D_ 104 (2021) 112010 [2108.07586].
  * [121] ATLAS collaboration, M. Aaboud et al., _Search for photonic signatures of gauge-mediated supersymmetry in 13 TeV $pp$ collisions with the ATLAS detector_, _Phys. Rev. D_ 97 (2018) 092006 [1802.03158].
  * [122] ATLAS collaboration, G. Aad et al., _Search for supersymmetry in events with four or more charged leptons in 139 fb -1 of $\sqrt{s}$ = 13 TeV pp collisions with the ATLAS detector_, _JHEP_ 07 (2021) 167 [2103.11684].
  * [123] ATLAS collaboration, G. Aad et al., _Searches for electroweak production of supersymmetric particles with compressed mass spectra in $\sqrt{s}=$ 13 TeV $pp$ collisions with the ATLAS detector_, _Phys. Rev. D_ 101 (2020) 052005 [1911.12606].
  * [124] ATLAS collaboration, M. Aaboud et al., _Search for electroweak production of supersymmetric states in scenarios with compressed mass spectra at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Phys. Rev. D_ 97 (2018) 052010 [1712.08119].
  * [125] CMS collaboration, A. M. Sirunyan et al., _Search for new physics in events with two soft oppositely charged leptons and missing transverse momentum in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _Phys. Lett. B_ 782 (2018) 440 [1801.01846].
  * [126] F. Feroz, M. P. Hobson and M. Bridges, _MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics_ , _Mon. Not. Roy. Astron. Soc._ 398 (2009) 1601 [0809.3437].
  * [127] F. Feroz, M. P. Hobson, E. Cameron and A. N. Pettitt, _Importance Nested Sampling and the MultiNest Algorithm_ , _Open J. Astrophys._ 2 (2019) 10 [1306.2144].
  * [128] C. K. Khosa, S. Kraml, A. Lessa, P. Neuhuber and W. Waltenberger, _SModelS Database Update v1.2.3_ , _LHEP_ 2020 (2020) 158 [2005.00555].
  * [129] W. Beenakker, R. Hopker and M. Spira, _PROSPINO: A Program for the production of supersymmetric particles in next-to-leading order QCD_ , hep-ph/9611232.
  * [130] J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer, _MadGraph 5 : Going Beyond_ , _JHEP_ 06 (2011) 128 [1106.0522].
  * [131] E. Conte, B. Fuks and G. Serret, _MadAnalysis 5, A User-Friendly Framework for Collider Phenomenology_ , _Comput. Phys. Commun._ 184 (2013) 222 [1206.1599].
  * [132] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159 [1410.3012].
  * [133] DELPHES 3 collaboration, J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens et al., _DELPHES 3, A modular framework for fast simulation of a generic collider experiment_ , _JHEP_ 02 (2014) 057 [1307.6346].
  * [134] M. Drees, H. Dreiner, D. Schmeier, J. Tattersall and J. S. Kim, _CheckMATE: Confronting your Favourite New Physics Model with LHC Data_ , _Comput. Phys. Commun._ 187 (2015) 227 [1312.2591].
  * [135] D. Dercks, N. Desai, J. S. Kim, K. Rolbiecki, J. Tattersall and T. Weber, _CheckMATE 2: From the model to the limit_ , _Comput. Phys. Commun._ 221 (2017) 383 [1611.09856].
  * [136] J. S. Kim, D. Schmeier, J. Tattersall and K. Rolbiecki, _A framework to create customised LHC analyses within CheckMATE_ , _Comput. Phys. Commun._ 196 (2015) 535 [1503.01123].
  * [137] J. Cao, J. Lian, Y. Pan, D. Zhang and P. Zhu, _Improved (g $-$ 2)μ measurement and singlino dark matter in $\mu$-term extended $\mathbb{Z}$3-NMSSM_, _JHEP_ 09 (2021) 175 [2104.03284].
  * [138] A. Berlin, P. Gratia, D. Hooper and S. D. McDermott, _Hidden Sector Dark Matter Models for the Galactic Center Gamma-Ray Excess_ , _Phys. Rev. D_ 90 (2014) 015032 [1405.5204].
  * [139] T. Gherghetta, B. von Harling, A. D. Medina, M. A. Schmidt and T. Trott, _SUSY implications from WIMP annihilation into scalars at the Galactic Center_ , _Phys. Rev. D_ 91 (2015) 105004 [1502.07173].