# Possible detection of the progenitor of the Type II supernova SN2023ixf

Joanne L Pledger Jeremiah Horrocks Institute, University of Central
Lancashire, Preston, PR1 2HE, UK  Michael M Shara  American Museum of Natural
History, New York City, NY 10024, USA

###### Abstract

Stellar evolution theory predicts multiple pathways to the explosive deaths of
stars as supernovae. Locating and characterizing the progenitors of well-
studied supernovae is important to constrain the theory, and to justify and
design future surveys to improve on progenitor detections. Here we report the
serendipitous pre-explosion imaging, by the Hubble Space Telescope, of
SN2023ixf, one of the nearest extragalactic supernovae ever discovered, in the
galaxy M101. The extremely red color and absolute magnitude MF814W =
–5.11${}^{+0.65}_{-0.47}$ mag suggest that the progenitor was a red
supergiant. Comparison with stellar evolutionary isochrones suggests it is
within the relatively low initial mass range of $\sim$8-10M⊙ and that there is
likely a lot of dust present at the supernova site.

supernovae: individual: 2023ixf

††facilities: HST(ACS/WFC3), Gemini(GMOS-N)††software: MULTIDRIZZLE
(https://www.stsci.edu/ koekemoe/multidrizzle/), IRAF (Tody, 1986), PARSEC
isochrones are freely available from http://stev.oapd.inaf.it/

##  1 Introduction

Supernovae are responsible for the production of many of the chemical elements
and dust in galaxies and throughout the universe (Burbidge et al., 1957;
Cameron, 1957; Hoyle & Wickramasinghe, 1970; Johnson et al., 2020), the
heating of the interstellar medium and galactic fountains (McKee & Cowie,
1975; Corbelli & Salpeter, 1988), cosmic rays (Baade & Zwicky, 1934; Blandford
& Ostriker, 1978), long-duration gamma-ray bursts (Cano et al., 2017) and
possibly even the formation of our own Solar System (Clayton, 1982). They are
the most luminous distance indicators known, responsible for the detection of
the dark energy (Riess et al., 1998; Perlmutter et al., 1999). Understanding
how supernovae come about is thus important to almost every field in
astrophysics. The most basic question that must be answered to make progress
in understanding supernovae is: What stars become supernovae?

Nearly 50 years ago Beatrice Tinsley (Tinsley, 1975) posed that question, and
summarized everything then known about the subject in an 11 page review
article. In the interim over 5600 refereed journal articles have greatly
advanced both observations and stellar evolution theory to yield some answers,
but half a century later there are still multiple types of supernovae for
which the answer remains “We still don’t know for certain”. This ignorance is
a major impediment to testing theories of the late stages of stellar evolution
in both single and binary stars (Smartt, 2009; Langer, 2012; Sukhbold et al.,
2016; Williams et al., 2019; Rodríguez, 2022; Liu et al., 2023).

To date there have been $\sim$30 direct detections of Type II SN progenitors
(see Van Dyk (2017) for a review), the first being SN1961V in NGC 1058
(Bertola, 1964; Zwicky, 1964) and the closest being SN1987A in the LMC
(Gilmozzi et al., 1987). Many pre-SN observations only yield an upper limit;
successful detections of progenitors are typically limited to $\sim$30 Mpc. At
a distance of just 6.4 Mpc (Shappee & Stanek, 2011) SN2023ixf in M101 is one
of the closest and brightest extra-galactic supernovae ever detected.
Intensive followup studies at many wavelengths are underway, in the X-ray
(ATel #16049; Grefenstette (2023) and ATel #16044; Kawai et al. (2023).),
optical (ATel #16045; Villafane et al. (2023)) and also a hunt for the
associated neutrinos (ATel #16043; Thwaites et al. (2023). This will be one of
the most intensively studied supernovae ever, thus characterizing its
progenitor is clearly of importance.

SN2023ixf in M101 was discovered by K. Itagaki on 19th May 2023 with a
position of $\alpha$(J2000) = 14:03:38.564 and $\delta$(J2000) = 54:18:42.02,
(Itagaki, 2023) and was classified as a Type II SN shortly after by Perley et
al. (2023). In this paper we use archival Hubble Space Telescope ACS/WFC
imaging to identify the possible progenitor of this event and use evolutionary
models to predict an initial mass.

##  2 Observations of the SN site & Data Reduction

As part of our ongoing survey of M101 for massive stars we had already
downloaded archival Hubble Space Telescope (HST) Advanced Camera for Surveys
(ACS) images taken with its F435W, F555W and F814W filters to combine with our
HST Wide Field Camera 3 (WFC3) data taken with a narrowband F469N filter,
tuned to the strong HeII emission line at 468.6 nm. The images were reduced
using the multidrizzle routine and mapped onto our WFC3/F469N narrow-band
imaging as described in Shara et al. (2013). Photometry was performed using
the standalone daophot code (Stetson, 1987) and a model point-spread function
(PSF) based on isolated point-like stars was built and applied to all other
stars in the field. Zero-points from the HST literature were applied to
convert our observed magnitudes into Vega magnitudes using the ACS Zeropoint
Calculator webtool 111The Zeropoint calculator is available at
https://www.stsci.edu/hst/instrumentation/acs/data-analysis/zeropoints

Using the methods outlined in Shara et al. (2013) we determine our photometric
detection limit for the F435W image to be mF435W = 28.5 mag, but although this
is the faintest object we can detect we do not always do so with 100%
efficiency. The magnitude at which we are confident we will always detect an
object, our 100$\%$ detection limit, is mF435W = 26.7 mag which. Using the
extinction law of Cardelli et al. (1989) and assuming an average extinction of
A(H$\alpha$) = 1.06 mag and a Milky Way foreground reddening of E(B–V) = 0.01
mag (Lee et al., 2009), corresponding to AF435W = 1.74 mag, we determine our
absolute magnitude 100% detection limit to be MF435W = -–4.07 mag. The mean
Galactic extinction curve does not work well for other galaxies, especially
with high star formation rates like M101 (Calzetti et al., 1994) so we use
values from Pang et al. (2016) for AF555W and AF814W which are presented along
with completeness limits for each filter in Table 1.

Table 1: Archival HST ACS/WFC data used for progenitor identification with the detection limits given in Vega magnitudes. All data were obtained under program ID 9490 (PI: Kuntz).  Filter | Date Obs | Exposure | Faintest | Aλ | 100$\%$ detection limit
---|---|---|---|---|---
|  | Time (sec) | Limit |  | Apparent | Absolute
F435W | 2002-11-15 | 1620 | 28.5 | 1.74 | 26.7 | –4.07
F555W | 2002-11-16 | 1440 | 28.0 | 1.36 | 26.6 | –3.79
F814W | 2002-11-16 | 1440 | 27.3 | 0.80 | 25.8 | –4.03

Whilst average extinction values are useful to determine our detection limit
for objects in M101, more local values of extinction are required for analysis
of the progenitor region and ultimately the determination of the initial mass.
The SN site lies $\sim$1 arcsecond from the HII region NGC 5461 which has
undergone analysis by Kennicutt & Garnett (1996) who analysed spectra of
regions 1105, 1098, 1086 and 1052 shown in their Figure 2b. In Table 2 we
present C(H$\beta$) and R23 values from Table 2. in Kennicutt & Garnett (1996)
and determine AF814W = 0.49${}^{+0.359}_{--0.158}$. We note that the R23
indicator is double–valued and thus we use the calibration from equation 8 in
Yin et al. (2007) to calculate the lower branch 12$+$log(O/H) and Pilyugin &
Thuan (2005) with an excitation parameter P=0.9 for the upper branch.
Metallicities derived from these equations for each of the four H II regions
local to the SN site are presented in Table 2. The average metallicity for the
lower branch is 12$+$log(O/H) = 7.485 and 12$+$log(O/H) = 8.583 for the upper
branch. We note that no H$\alpha$ fluxes are presented in Kennicutt & Garnett
(1996) so the O3N2 or N2 calibrators from e.g. Pettini & Pagel (2004) cannot
be used. However, metallicity determinations of other nearby H II regions,
albeit further from the SN site, presented in Pledger et al. (2018) using the
O3N2 calibration suggests the upper branch is more likely.

Table 2: Four HII regions close to the site of SN2023ixf. C(H$\beta$) and R23 are taken from Kennicutt & Garnett (1996) and we note ID names correspond to those in Hodge et al. (1990). ID | C(H$\beta$) | E(B-V) | AF814W | R23 | 12$+$log(O/H)lower | 12$+$log(O/H)upper
---|---|---|---|---|---|---
1105 | 0.36 | 0.277 | 0.413 | 6.01 | 7.58 | 8.53
1098 | 0.74 | 0.570 | 0.849 | 4.50 | 7.40 | 8.63
1086 | 0.32 | 0.246 | 0.362 | 5.02 | 7.47 | 8.59
1052 | 0.29 | 0.223 | 0.332 | 5.22 | 7.49 | 8.58
Average | 0.428 | 0.329 | 0.490 | 5.188 | 7.485 | 8.583

##  3 Identification of the SN progenitor

To identify any potential progenitors of SN2023ixf we must compare the
location of the SN with pre-SN imaging which in turn requires a consistent co-
ordinate system between the images. From our previous work our HST images in
the three different filters have already been aligned using a geometric
transformation from the GEOMAP routine within the HST MULTIDRIZZLE software
(Shara et al., 2013). We applied the same method to r–band Gemini/GMOS-N
observations of SN2023ixf (PI: Lotz, Program ID: GN-2003A-DD-105) taken on 05
June 2023. We compared the positions of 30 sources common to both images which
yielded a geometric transformation of the Gemini image onto the F814W/ACS
image with an RMS error of $\pm$24mas. We applied this transformation using
GEOTRAN, and then compared co-ordinates of our targets again, finding a
standard deviation of $\sigma$ = 0.03 arseconds. This allows us to determine
an accurate position for the SN in the archival imaging as shown in Figure 1.

Figure 1: HST/ACS F814W image taken in November 2002 (left) and Gemini/GMOS
r-band image taken in June 2023 (right) at the site of SN2023ixf in M101 which
became available after the initial submission of this paper. The location of
the SN is highlighted and is co-incident with a single source in the F814W
image. We note the SN is on the edge of one of the GMOS-N chips.

Using this position and error we detected no progenitor in the F435W image,
with the nearest source $\sim$0.44 arcseconds to the south-west of the
supernova location. The F555W image also shows no obvious photometric
detection, though one could argue at a hint of a source at the detection limit
but slightly off centred from the SN location as shown in Figure 2.
Unfortunately our WFC3/F469N images just miss the SN site, however a Type II
progenitor would not be expected to be a star with strong helium emission such
as a Wolf Rayet star.

The F814W image reveals a source co-incident with the SN location, within our
positional error limit (Figure 2). Inspection of the profile of the source
reveals an asymmetric profile as shown in Figure 3 suggesting that the source
is two partially resolved stars which would be consistent with the hint of an
offset source in the F555W image. This is supported further from comparison
with a nearby point-like star which does not show any evidence of an
asymmetric profile. At a distance of 6.4 Mpc, the 0.15 arcsecond spatial
resolution of our drizzled HST/ACS images corresponds to a physical size of
4.65 pc which could easily hide multiple stars. However, our positional error
suggests that the more prominent source in F814W is co-incident with the SN
site within 1$\sigma$ whereas the fainter, more easterly source is further
away.

Figure 2: Hubble Space Telescope ACS/WFC images of the $\sim$2 arcsecond2
region around the site of SN2023ixf in F435W (top left), F555W (top right) and
F814W (bottom) filters. The location of the two sources identified in the
F814W image are marked with a circle (with the size corresponding to 1$\sigma$
positional error) whilst the SN location is indicated by an X. We note that
North is up and East is left.

Figure 3: Mean X-profile of the source at the location of SN2023ixf (left).
The peak at x$\sim$3517 corresponds to the brightest part of the source and
another peak is suggested at x$\sim$3515 for the second, fainter source. For
context the profile of a nearby point-like star is shown (right) for
comparison.

Our photometry of these two partially unresolved sources shows that the
westerly source (located at $\alpha$ = 14:03:38.544 and $\delta$ = 54:18:41.90
) has an apparent Vega magnitude of mF814W = 24.41$\pm$0.06 mag, which
assuming our average extinction AF814W = 0.49${}^{+0.359}_{-0.158}$ mag from
Table 2 and a distance of 6.4$\pm$0.7 Mpc (Shappee & Stanek, 2011),
corresponds to an absolute magnitude of MF814W = –5.11${}^{+0.65}_{-0.47}$
mag. Following the initial submission of this paper we note that Soraisam et
al. (2023) (ATel $\\#16050$) find an apparent magnitude of mF814W =
24.39$\pm$0.08 mag consistent with ours and that Jacobson-Galan et al. (2023)
find E(B-V) = 0.033 mag for SN2023ixf based on optical spectra and fitting of
the NaI D absorption line which would be most consistent with the lower end of
our absolute magnitude. The other possible progenitor is the eastern source
(located at $\alpha$ = 14:03:38.552 and $\delta$ = 54:18:41.90). It is not
well-enough resolved for a magnitude to be determined.

Figure 4: Color-magnitude diagram plotting 80-200Myr PARSEC stellar
evolutionary isochrones along with our progenitor candidate at a metallicity
of log(O/H)$+$12 = 8.58. The V-I plotted is a lower limit using the 100%
detection limit of our F555W images and as such we expect the true value to be
further to the right on this plot. We also note that the star will should
explode at the end of the evolutionary track.

##  4 Mass Estimate

To determine a mass estimate for the progenitor we used PARSEC222PARSEC
isochrones are freely available from http://stev.oapd.inaf.it/ stellar
evolutionary isochrones with ages 50–200 Myr from Bressan et al. (2012). We
ran models based on both our upper and lower metallicity values using our
average extinction value of AF814W = 0.49 mag and use our absolute F814W
magnitude of –5.11 mag, assuming a lower limit of V – I $\sim$ 2.2 mag based
on our 100% detection limit in F555W of 26.6 mag to assess the initial mass of
the progenitor. Error bars shown in Figure 4 are based on errors in distance,
extinction and magnitude. Figure 4 show the stellar evolutionary tracks that
best match our observations assuming a metallicity of log(O/H)$+$12 = 8.58.
Within our error bars the progenitor star is $\sim$60-200 Myr old which
corresponds to an initial mass of $\sim$4–7M ⊙ given the star would explode at
the end of the track. We note that stars in the lower end of this mass range
would not be expected to explode as a RSG. Similarly, if we assume the lower
metallicity of log(O/H)$+$12 = 7.49 then the luminosity falls below that
expected to explode as a RSG, consistent with a comparison with BPASS models
(J. Eldridge; private communication; Stanway & Eldridge (2018)). We propose
that the likely explanation for this low mass is that there is a lot of dust
that we don’t account for and that the progenitor’s true initial mass is in
the lower mass range for RSG of 8-10 M⊙. The presence of dust is consistent
with the findings of Szalai & Dyk (2023)(ATel $\\#$16042) who suggest a mass
of $\sim$15M⊙ from Spitzer data. Further analysis of the model-dependent dust,
combined with light curve modelling, is required to fully assess the nature of
the progenitor.

##  5 Conclusions

We have detected a possible optical counterpart to the progenitor of the type
IIP supernova 2023ixf in archival Hubble Space Telescope images. The
counterpart is not seen in (blue) F435W or (visible) F555W filter images, but
is easily detected in (NIR) F814W images. The red color and brightness of the
progenitor are consistent with it being a red supergiant of at least 7M⊙,
likely higher at $\sim$8–10⊙ if there is significant dust present. There is
also an extremely faint, unresolved source that could also potentially be the
progenitor of SN2023ixf although our astrometric calibration suggests this is
unlikely. There is a suggestion of a detection in the F555W image offset from
the F814W source but this would suggest a bluer object not typical of Type II
SN progenitors. Post-SN imaging with HST will be crucial to confirm the
disappearance of the progenitor, as is modelling of the SN light-curve and the
dust surrounding the supernova.

##  6 Acknowledgments

MMS thanks James Garland for alerting him to this nearby supernova. We thank
David Zurek for advice on magnitude system conversions and Anne Sansom for
advice on stellar isochrones. We also thank Jan Eldridge for discussion of the
progenitor mass in the BPASS models. Some of the data presented in this paper
were obtained from the Mikulski Archive for Space Telescopes (MAST) at the
Space Telescope Science Institute. The data used from HST is available at
MAST: https://doi.org/10.17909/qbmp-nw13 (catalog 10.17909/qbmp-nw13). Some of
this work was supported by the generous support of the late Hilary and Ethel
Lipsitz and we gratefully acknowledge them for their support. This research is
based on NASA/ESA Hubble Space Telescope observations obtained at the Space
Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy Inc. under NASA contract NAS5-26555.
Based on observations obtained at the international Gemini Observatory, a
program of NSF’s NOIRLab, which is managed by the Association of Universities
for Research in Astronomy (AURA) under a cooperative agreement with the
National Science Foundation on behalf of the Gemini Observatory partnership:
the National Science Foundation (United States), National Research Council
(Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio
de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência,
Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space
Science Institute (Republic of Korea).

## References

  * Baade & Zwicky (1934) Baade, W., & Zwicky, F. 1934, Proceedings of the National Academy of Science, 20, 259, doi: 10.1073/pnas.20.5.259
  * Bertola (1964) Bertola, F. 1964, Annales d’Astrophysique, 27, 319
  * Blandford & Ostriker (1978) Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29, doi: 10.1086/182658
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Burbidge et al. (1957) Burbidge, E. M., Burbidge, G. R., Fowler, W. A., & Hoyle, F. 1957, Reviews of Modern Physics, 29, 547, doi: 10.1103/RevModPhys.29.547
  * Calzetti et al. (1994) Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582, doi: 10.1086/174346
  * Cameron (1957) Cameron, A. G. W. 1957, PASP, 69, 201, doi: 10.1086/127051
  * Cano et al. (2017) Cano, Z., Wang, S.-Q., Dai, Z.-G., & Wu, X.-F. 2017, Advances in Astronomy, 2017, 8929054, doi: 10.1155/2017/8929054
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Clayton (1982) Clayton, D. D. 1982, in NATO Advanced Study Institute (ASI) Series C, Vol. 90, Supernovae: A Survey of Current Research, ed. M. J. Rees & R. J. Stoneham, 535–564, doi: 10.1007/978-94-009-7876-8_34
  * Corbelli & Salpeter (1988) Corbelli, E., & Salpeter, E. E. 1988, ApJ, 326, 551, doi: 10.1086/166116
  * Gilmozzi et al. (1987) Gilmozzi, R., Cassatella, A., Clavel, J., et al. 1987, Nature, 328, 318, doi: 10.1038/328318a0
  * Grefenstette (2023) Grefenstette, B. 2023, The Astronomer’s Telegram, 16049, 1
  * Hodge et al. (1990) Hodge, P. W., Gurwell, M., Goldader, J. D., & Kennicutt, Robert C., J. 1990, ApJS, 73, 661, doi: 10.1086/191483
  * Hoyle & Wickramasinghe (1970) Hoyle, F., & Wickramasinghe, N. C. 1970, Nature, 226, 62, doi: 10.1038/226062a0
  * Itagaki (2023) Itagaki, K. 2023, Transient Name Server Discovery Report, 2023-1158, 1
  * Jacobson-Galan et al. (2023) Jacobson-Galan, W. V., Dessart, L., Margutti, R., et al. 2023, arXiv e-prints, arXiv:2306.04721, doi: 10.48550/arXiv.2306.04721
  * Johnson et al. (2020) Johnson, J. A., Fields, B. D., & Thompson, T. A. 2020, Philosophical Transactions of the Royal Society of London Series A, 378, 20190301, doi: 10.1098/rsta.2019.0301
  * Kawai et al. (2023) Kawai, N., Serino, M., Negoro, H., et al. 2023, The Astronomer’s Telegram, 16044, 1
  * Kennicutt & Garnett (1996) Kennicutt, Robert C., J., & Garnett, D. R. 1996, ApJ, 456, 504, doi: 10.1086/176675
  * Langer (2012) Langer, N. 2012, ARA&A, 50, 107, doi: 10.1146/annurev-astro-081811-125534
  * Lee et al. (2009) Lee, J. C., Gil de Paz, A., Tremonti, C., et al. 2009, ApJ, 706, 599, doi: 10.1088/0004-637X/706/1/599
  * Liu et al. (2023) Liu, Z.-W., Roepke, F. K., & Han, Z. 2023, arXiv e-prints, arXiv:2305.13305, doi: 10.48550/arXiv.2305.13305
  * McKee & Cowie (1975) McKee, C. F., & Cowie, L. L. 1975, ApJ, 195, 715, doi: 10.1086/153373
  * Pang et al. (2016) Pang, X., Pasquali, A., & Grebel, E. K. 2016, AJ, 151, 23, doi: 10.3847/0004-6256/151/2/23
  * Perley et al. (2023) Perley, D., Meynardie, W., Chu, M., & Fremling, C. 2023, Transient Name Server Classification Report, 2023-146, 1
  * Perlmutter et al. (1999) Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565, doi: 10.1086/307221
  * Pettini & Pagel (2004) Pettini, M., & Pagel, B. E. J. 2004, MNRAS, 348, L59, doi: 10.1111/j.1365-2966.2004.07591.x
  * Pilyugin & Thuan (2005) Pilyugin, L. S., & Thuan, T. X. 2005, ApJ, 631, 231, doi: 10.1086/432408
  * Pledger et al. (2018) Pledger, J. L., Shara, M. M., Wilde, M., et al. 2018, MNRAS, 473, 148, doi: 10.1093/mnras/stx2190
  * Riess et al. (1998) Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009, doi: 10.1086/300499
  * Rodríguez (2022) Rodríguez, Ó. 2022, MNRAS, 515, 897, doi: 10.1093/mnras/stac1831
  * Shappee & Stanek (2011) Shappee, B. J., & Stanek, K. Z. 2011, ApJ, 733, 124, doi: 10.1088/0004-637X/733/2/124
  * Shara et al. (2013) Shara, M. M., Bibby, J. L., Zurek, D., et al. 2013, AJ, 146, 162, doi: 10.1088/0004-6256/146/6/162
  * Smartt (2009) Smartt, S. J. 2009, ARA&A, 47, 63, doi: 10.1146/annurev-astro-082708-101737
  * Soraisam et al. (2023) Soraisam, M., Matheson, T., Andrews, J., et al. 2023, The Astronomer’s Telegram, 16050, 1
  * Stanway & Eldridge (2018) Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 479, 75, doi: 10.1093/mnras/sty1353
  * Stetson (1987) Stetson, P. B. 1987, PASP, 99, 191, doi: 10.1086/131977
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38, doi: 10.3847/0004-637X/821/1/38
  * Szalai & Dyk (2023) Szalai, T., & Dyk, S. V. 2023, The Astronomer’s Telegram, 16042, 1
  * Thwaites et al. (2023) Thwaites, J., Vandenbroucke, J., Santander, M., & IceCube Collaboration. 2023, The Astronomer’s Telegram, 16043, 1
  * Tinsley (1975) Tinsley, B. M. 1975, PASP, 87, 837, doi: 10.1086/129857
  * Tody (1986) Tody, D. 1986, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 627, Instrumentation in astronomy VI, ed. D. L. Crawford, 733, doi: 10.1117/12.968154
  * Van Dyk (2017) Van Dyk, S. D. 2017, Philosophical Transactions of the Royal Society of London Series A, 375, 20160277, doi: 10.1098/rsta.2016.0277
  * Villafane et al. (2023) Villafane, J., Nespral, D., Zamora, O., et al. 2023, The Astronomer’s Telegram, 16045, 1
  * Williams et al. (2019) Williams, B. F., Hillis, T. J., Blair, W. P., et al. 2019, ApJ, 881, 54, doi: 10.3847/1538-4357/ab2190
  * Yin et al. (2007) Yin, S. Y., Liang, Y. C., Hammer, F., et al. 2007, A&A, 462, 535, doi: 10.1051/0004-6361:20065798
  * Zwicky (1964) Zwicky, F. 1964, ApJ, 139, 514, doi: 10.1086/147779


# A random-key GRASP for combinatorial optimization

A.A. Chaves1,∗, M.G.C. Resende1,2, R.M.A. Silva3

1Departamento de Ciência e Tecnologia, Universidade Federal de São Paulo, São
José dos Campos, SP, Brazil
2Industrial and Systems Engineering, University of Washington, Seattle, WA,
U.S.A.
3Centro de Informática, Universidade Federal de Pernambuco, Recife, PE, Brazil

Abstract. This paper proposes a problem-independent GRASP metaheuristic using
the random-key optimizer (RKO) paradigm. GRASP (greedy randomized adaptive
search procedure) is a metaheuristic for combinatorial optimization that
repeatedly applies a semi-greedy construction procedure followed by a local
search procedure. The best solution found over all iterations is returned as
the solution of the GRASP. Continuous GRASP (C-GRASP) is an extension of GRASP
for continuous optimization in the unit hypercube. A random-key optimizer
(RKO) uses a vector of random keys to encode a solution to a combinatorial
optimization problem. It uses a decoder to evaluate a solution encoded by the
vector of random keys. A random-key GRASP is a C-GRASP where points in the
unit hypercube are evaluated employing a decoder. We describe random key GRASP
consisting of a problem-independent component and a problem-dependent decoder.
As a proof of concept, the random-key GRASP is tested on five NP-hard
combinatorial optimization problems: traveling salesman problem, tree of hubs
location problem, Steiner triple covering problem, node capacitated graph
partitioning problem, and job sequencing and tool switching problem.

Keywords. GRASP; combinatorial optimization; random-key optimizer.

11footnotetext: Corresponding author. E-mail addresses:
<EMAIL_ADDRESS>(A. Chaves<EMAIL_ADDRESS>(M.G.C. Resende),
<EMAIL_ADDRESS>(R.M.A. Silva).

##  1\. Introduction

In this paper, a problem-independent GRASP metaheuristic is proposed using the
random-key optimizer (RKO) paradigm [50]. GRASP (greedy randomized adaptive
search procedure) is a metaheuristic for combinatorial optimization that
repeatedly applies a semi-greedy construction procedure followed by a local
search procedure [19, 20, 49]. The best solution found over all iterations is
returned as the solution of the GRASP. Continuous GRASP (C-GRASP) is an
extension of GRASP for continuous optimization in the unit hypercube [26]. A
random-key optimizer (RKO) makes use of a vector of random keys [5] to encode
a solution to a combinatorial optimization problem. It uses a decoder to
evaluate a solution encoded by the vector of random keys. A random-key GRASP
(RK-GRASP) is a C-GRASP where points in the unit hypercube are evaluated by
means of a decoder. We describe the random-key GRASP metaheuristic consisting
of a problem-independent component and a problem-dependent decoder. As a proof
of concept, the random-key GRASP is tested on five NP-hard combinatorial
optimization problems: traveling salesman problem, tree of hubs location
problem, Steiner triple covering problem, node capacitated graph partitioning
problem, and job sequencing and tool switching problem.

The paper is organized as follows. In Section 2, the concept of Random-Key
Optimizer (RKO) and its application in a RK-GRASP are described. Section 3
describes the implementation of RK-GRASP, including constructive and local
search phases. Section 4 presents five combinatorial optimization problems to
which RK-GRASP is applied, simply by making use of different
decoders.Experimental results obtained by applying RK-GRASP to the five
combinatorial optimization problems are summarized in Section 5. Concluding
remarks are made in Section 6, where the key findings and contributions of the
paper are summarized.

##  2\. Random-Key optimizer

Random-key genetic algorithms (RKGA) were first introduced by Bean [5]. In an
RKGA, solutions are encoded as vectors of random keys, i.e., randomly
generated real numbers in the interval $[0,1)$. A population of $p$ random-key
vectors is evolved over a number of generations. The initial population
consists of randomly generated $n$-vectors. At each iteration, each random-key
vector is decoded and evaluated. The decoder is a deterministic algorithm,
usually a heuristic, that takes as input a vector of random keys and returns a
feasible solution of the problem being solved along the cost of the solution.
The population is partitioned into an elite set of $p_{e}$ solutions and a
non-elite set of $p-p_{e}$ solutions. The population is evolved as follows.
First, all elite solutions are copied to the next generation, and $p_{m}$
mutant solutions, each consisting of an $n$-vector of randomly generated keys,
are added to the next generation. Finally, the remaining $p-p_{e}-p_{m}$
elements are made of solutions produced by a crossover of pairs of solutions
drawn at random from the population of the current generation. Mating is done
with parameterized uniform crossover [53] where the child inherits the key of
one of the parents with probability $\rho$ and of the other with probability
$1-\rho$. We call $\rho$ the inheritance probability. Biased random-key
genetic algorithms [24] is an extension of RKGA where one of the parents is
always selected from the elite set of the population, and that parent is
associated with the inheritance probability $\rho>\frac{1}{2}$.

Both RKGA and BRKGA are problem-independent algorithms in the sense that there
is a clear separation between the solver and the problem being solved where
the connection of the solver with the problem is done by way of a decoder. For
each type of problem to be solved, a new decoder is implemented to serve as
the link between the solver and the problem. These algorithms are examples of
Random-Key Optimizers (RKO) [50]. An RKO is an optimization heuristic
algorithm that solves a discrete optimization problem indirectly in the
continuous unit $n$-dimensional hypercube $\mathcal{H}_{n}$. For each solution
$x\in\mathcal{H}_{n}$, a decoder $\mathcal{D}$ maps $x$ to the solution
$\mathcal{D}(x)$ in the solution space of the discrete optimization problem.
With such a separation of solver and decoder, one needs only to implement the
solver once and then it can be reused to solve a number of different problems
by simply devising a decoder for that problem. Examples of Application
Programming Interfaces (APIs) for BRKGA are Toso and Resende [55], Andrade et
al. [3], Oliveira et al. [41], and Chaves and Lorena [6].

The first paper to extend the concept of RKO to other algorithmic frameworks
was Schuetz et al. [50] where in addition to a BRKGA for robot motion
planning, the authors propose an RKO using dual annealing, an extension of
generalized simulated annealing [56, 57]. Recently, RKO has been implemented
in simulated annealing, iterated local search, and variable neighborhood
search for the tree hub location problem [33].

##  3\. Random-Key GRASP

This section introduces random-key GRASP, or RK-GRASP, a problem-independent
GRASP that solves discrete optimization problems through continuous
optimization and a decoder. An advantage of using this GRASP is that the user
only needs to implement a decoder since the algorithm-specific components of
GRASP are implemented as an API. When using a standard GRASP the user needs to
tailor both the semi-greedy construction and the local search for the problem
on hand. We begin by reviewing the basic GRASP for combinatorial optimization
and then consider C-GRASP or continuous GRASP for continuous global
optimization over the continuous $n$-dimensional hypercube.

###  3.1. GRASP

The metaheuristic Greedy Randomized Adaptive Search Procedure (GRASP) [49] was
introduced by Feo and Resende [19, 20]. A GRASP is a multi-start procedure in
which at each iteration a semi-greedy solution is constructed and local search
is applied to this solution. In greedy construction, one builds a solution one
element at a time. Among the set of elements whose inclusion in the solution
does not lead to infeasibility, the greedy algorithm selects an element that
best contributes to the objective function. In a semi-greedy algorithm,
instead of selecting a single element, a set of elements is selected, each
having the best or almost the best contribution to the objective function,
none of which leads to infeasibility. The set of elements from which to select
the single element is called the restricted candidate list, or simply RCL. By
adjusting the amount of randomness in the selection process, one can mimic a
totally random selection or a totally greedy process. Neither the greedy nor
semi-greedy construction process is guaranteed to produce a solution that is
globally optimal, or even locally optimal with respect to some local
neighborhood structure. Given a solution $x$ produced by a greedy or semi-
greedy construction procedure, let $\mathcal{N}(x)$ denote the local
neighborhood of $x$, i.e. the set of solutions obtained by making a small
perturbation (or change) to $x$. In the case of minimization (maximization),
local search stops at local minimum (maximum) $x$ if there is no solution
$y\in\mathcal{N}(x)$ such the cost $f(y)<f(x)$ (or $f(y)>f(x)$ for the case of
maximization). If such a solution exists, we set $x\leftarrow y$ and the
search resumes. The best locally optimal solution visited over all iterations
is returned as the GRASP solution.

###  3.2. Continuous GRASP

Continuous GRASP, or simply C-GRASP, is an extension of GRASP for solving
continuous optimization problems subject to box or simple bounding constraints
[26, 27],

$\min_{x\in\mathbb{R}^{n}}\\{f(x)\;|\;L_{n}\leq x\leq U_{n}\\},$

where $L_{n}$ and $U_{n}$ are vectors of lower and upper bounds on $x$, i.e.
$L_{n}(i)\leq x_{i}\leq U_{n}(i),$ for $i=1,\ldots,n$ and $f(x)$ is the cost
of solution $x$. Cost $f(x)$ can be evaluated in a multitude of ways, e.g.
analytically, through simulation, via a mathematical program, or with a
decoder. The objective is to find a global optimum. Like GRASP, C-GRASP is a
multi-start procedure in which each iteration consists of a construction, or
diversification, phase followed by a local search, or intensification, phase.
C-GRASP evaluates points on a dynamic grid, with grid size initially set to
$h=h_{0}$. Each construction phase starts at the current solution $x$
(initially a random point $x\in\mathbb{R}^{n}\;|\;L_{n}\leq x\leq U_{n}$). To
build an RCL, C-GRASP performs a line search on $f(x)$ in each direction
$e_{i}=(0,0,\ldots,0,1,0,\ldots,0,0)$, for $i=1,2,\ldots,n$, where the only
nonzero is a 1 in position $i$. For $i=1,2,\ldots,n$, the line search is
limited to evaluating $f(x+e_{i}\cdot h\cdot k)$ for all values of
$k\in\\{0,1,-1,2,-2,\ldots\\}$ such that $L_{i}\leq x_{i}+e_{i}\cdot h\cdot
k\leq U_{i}$. The result of line search $i$ is $z_{i}$ with cost
$g_{i}=f(z_{i})$. Let $g_{m}=\min\\{g_{i}\;|\;i=1,\ldots,n\\}$ and
$g_{M}=\max\\{g_{i}\;|\;i=1,\ldots,n\\}$. The best line search solutions, i.e.
those with $g_{i}\leq(1-\alpha)\cdot g_{m}+\alpha\cdot g_{M}$ for some
$\alpha\in[0,1]$ are placed in an RCL and an index $j$ is selected at random
to be fixed in $x$ with $x=z_{j}$. Direction $e_{j}$ is flagged to no longer
be explored in this construction iteration. This is repeated until a solution
is constructed.

Once a semi-greedy solution $x$ has been constructed, a local search, or
intensification, phase is applied around $x$. Several implementations of local
search have been described. Suppose the current semi-greedy solution is
$\bar{x}$. In the first paper, Hirsch et al. [26] suggest examining a given
maximum number of points MaxDirToTry of the form
$\bar{x}+h\cdot\\{-1,0,1\\}^{n}$. This local search examines only grid points.
Hirsch, Pardalos, and Resende [27] sample a user-defined number of feasible
grid points and project each one of them onto the surface of the hyper-sphere
of radius $h$, centered at $\bar{x}$. Each point is evaluated as it is
projected, and a first-improvement policy is applied. If an improving point is
found, then $\bar{x}$ is set to this improving solution. If no improving
solution is found, the grid size is reduced by a specified factor, and the
search process continues. See [26, 27, 52] for details.

###  3.3. Random-Key GRASP

We consider optimization problems of the form

$\min\\{f(\mathcal{D}(x))\mid x\in\mathcal{H}_{n}\\},$

where $\mathcal{H}_{n}$ denotes the unit hypercube in $\mathbb{R}^{n}$ and
$f(\mathcal{D}(x))$ represents the cost value returned by a decoder
$\mathcal{D}$ when given $x\in\mathcal{H}_{n}$ as input.

RK-GRASP, or random-key GRASP, embodies the core principles of both the GRASP
and C-GRASP, encompassing constructive and local search phases. This paper
introduces a constructive phase inspired by the line search strategy of
C-GRASP. Furthermore, we introduce three distinct methods for the local search
phase: Grid Search, Nelder-Mead Search, and Random Variable Neighborhood
Descent (RVND). Subsequently, we provide the pseudo-codes of these methods and
a summary of their concepts.

Algorithm 1 presents pseudo-code for RK-GRASP. Initially, a random-key vector
$x$ is generated randomly, and the parameter value $h$ is set as $h_{s}$. The
RK-GRASP search process continues while $h$ is greater than or equal to the
lower limit ($h_{e}$). Like GRASP, each iteration consists of constructive and
local search phases. Unlike GRASP, the constructive phase depends on the local
optimum found in the local search phase of the previous iteration. In the
constructive phase, a new solution $x^{\prime}$ is generated from the current
solution $x$ using partially greedy moves. The local search phase intensifies
the constructed solution and obtains a local optimal solution
$x^{\prime\prime}$. If the objective function value of the $x^{\prime\prime}$
is better than the best solution found so far, $x^{*}$ is updated. Otherwise,
the grid is reduced by halving the value of $h$. Finally, the algorithm
evaluates whether the new solution should be accepted. While a straightforward
acceptance criterion could be employed, such as accepting all solutions as in
C-GRASP, the RK-GRASP uses a simulated annealing acceptance criterion to
enhance the solution quality. The algorithm continues until a stopping
criterion is met.

Data: $n,hs,he$

Result: Best solution found during the search process ($x^{*}$)

1 $f(\mathcal{D}(x^{*}))\leftarrow\infty$;

2 while  _not stopping condition_ do

3 $x\leftarrow$ CreateInitialSolution($n$);

4 $h\leftarrow hs$;

5 while  _$h\geq he$_ do

6 $x^{\prime}\leftarrow$ ConstructGreedyRandomized($x,h,n$);

7 $x^{\prime\prime}\leftarrow$ LocalSearch($x^{\prime},h,n$);

8 if  _$f(\mathcal{D}(x^{\prime\prime})) <f(\mathcal{D}(x^{*}))$_ then

9 $x^{*}\leftarrow x^{\prime\prime}$;

10

11 else

12 $h\leftarrow h/2$;

13

14 if  _accept( $x^{\prime\prime},x$)_ then

15 $x\leftarrow x^{\prime\prime}$;

16

17

return $x^{*}$

Algorithm 1 RK-GRASP

The constructive phase involves perturbing a solution $x$ with semi-greedy
moves that are randomly selected. Algorithm 2 presents the pseudocode for the
constructive phase. The $\mathit{UnFixed}$ vector is initialized with the
indices of all random keys in the solution vector $x$. The constructive phase
continues as long as the $\mathit{UnFixed}$ vector is not empty. In each
iteration, a LineSearch method is applied to each random key whose index $i$
is in $\mathit{UnFixed}$, returning the best objective function value $g_{i}$
for that random key obtained through a linear search over the interval $[0,1)$
and its corresponding value $r_{i}$. Subsequently, the Restricted Candidate
List (RCL) is constructed, considering the indices of the random keys that
generated solutions with objective function values within the range defined by
$\alpha$. The parameter $\alpha$ is randomly set at the beginning of the
construction phase within the interval $[0,1]$, where 0 represents the
insertion of only the index of the best random key, and 1 represents the
insertion of all indices that are $\mathit{UnFixed}$. From the RCL, an index
$j$ is randomly selected, the random key value of this index is updated with
the value found by LineSearch, and this index is removed from
$\mathit{UnFixed}$. To speed up the construction process [27], we use the
$\mathit{Reuse}$ flag to bypass the application of the LineSearch in
subsequent iterations where the random key has not been modified (i.e., when
$x_{j}=r_{j}$).

Data: $x,h,n$

Result: A constructive semi-greedy solution.

1 $\mathit{UnFixed}\leftarrow\\{1,2,...,n\\}$;

2 $\alpha\leftarrow$ UnifRand(0,1);

3 $\mathit{Reuse}\leftarrow\mathit{false}$;

4 while  _$\mathit{UnFixed}\neq\emptyset$_ do

5 $min\leftarrow+\infty$; $max\leftarrow-\infty$;

6 for  _$i=1,...,n$_ do

7 if  _$i\in\mathit{UnFixed}$_ then

8 if  _$\mathit{Reuse}=\mathit{false}$_ then

9 $[r_{i},g_{i}]\leftarrow$ LineSearch($x,h,i$);

10

11 if  _$min >g_{i}$_ then

12 $min\leftarrow g_{i}$;

13 if  _$max <g_{i}$_ then

14 $max\leftarrow g_{i}$;

15

16

17 $RCL\leftarrow\emptyset$;

18 for  _$i=1,...,n$_ do

19 if  _$i\in\mathit{UnFixed}\wedge g_{i}\leq min+\alpha\cdot(max-min)$_ then

20 $RCL\leftarrow RCL\cup\\{i\\}$;

21

22

23 $j\leftarrow$ RandomlySelectElement($RCL$);

24 if  _$x_{j}=r_{j}$_ then

25 $Reuse\leftarrow true$;

26

27 else

28 $x_{j}\leftarrow r_{j}$;

29 $f(\mathcal{D}(x))\leftarrow g_{j}$;

30 $\mathit{Reuse}\leftarrow\mathit{false}$;

31

32 $\mathit{UnFixed}\leftarrow\mathit{UnFixed}\backslash\\{j\\}$;

33

return $x$

Algorithm 2 ConstructGreedyRandomized

The LineSearch begins by generating a set of integers $K$ within the range
$[-1/h,1/h]$. The next step involves generating a set of feasible random keys
$RK$ derived from the initial random key $x_{i}$ and the values $K_{j}$ and
$h$ ($0\leq x_{i}+K_{j}\cdot h<1$). The random key $i$ is then updated with
each value contained in $RK$, and the solution $x$ is decoded to evaluate its
objective function value. In this paper, we propose adapting the greedy
algorithm so that instead of selecting the best options, it only considers a
subset of $q<|RK|$ possible random-key values (chosen randomly) during each
iteration. LineSearch returns, from this subset, the random-key value that
generates the solution with the best objective function value
($\mathit{best}_{r}$ and $\mathit{best}_{f}$).

Data: $x,h,i$

Result: The best value for random key $i$.

1 $K\leftarrow\\{0,1,-1,2,-2,...,1/h,-1/h\\}$;

2 $RK\leftarrow\emptyset$;

3 for  _$j=1,...,|K|$_ do

4 if  _$0\leq(x_{i}+K_{j}\cdot h) <1$_ then

5 $RK\leftarrow RK\cup(x_{i}+K_{j}\cdot h)$;

6

7

8$q\leftarrow\lceil log_{2}(1/h)\rceil$ \+ 1;

9 $\mathit{best}_{f}\leftarrow\infty$; $\mathit{best}_{r}\leftarrow 0$;

10 for  _$t=1,...,q$_ do

11 $j\leftarrow$ UnifRand(1,$|RK|$);

12 $x_{i}\leftarrow RK_{j}$;

13 if  _$f(\mathcal{D}(x)) <\mathit{best}_{f}$_ then

14 $\mathit{best}_{f}\leftarrow f(\mathcal{D}(x))$;

15 $\mathit{best}_{r}\leftarrow x_{i}$;

16

17 $RK\leftarrow RK\backslash\\{j\\}$

return $\mathit{best}_{r},\mathit{best}_{f}$

Algorithm 3 LineSearch

The Local Search phase intensifies the search in the neighborhood of the
constructed solution, seeking improved solutions. RK-GRASP can incorporate
different local search strategies. In this paper, we developed three local
search heuristics: Grid Search, Nelder-Mead Search, and Random Variable
Neighborhood Descent (RVND). These local search heuristics are problem-
independent and intensify the search within the random-key solution space,
utilizing the decoder to evaluate neighboring solutions.

Algorithm 4 presents pseudo-code for the Grid Search heuristic [28]. This
heuristic generates a neighborhood and determines at which points in the
neighborhood, if any, the objective function improves. If an improvement point
is found, it is made the current point, and the local search continues with
the new solution. We use the Algorithm 5 to create a solution in the $h$
neighborhood of the current solution $\bar{x}$. Consider the set of points
$v_{i}$ that are integer steps (of size $h$) away from $\bar{x}$. The
$h$-Neighborhood algorithm defines the projection of these points onto the
hyper-sphere of radius $h$, centered at $\bar{x}$ . We use the parameter
$\theta$ to control the intensity of this projection, maintaining a random key
fixed with a probability of $1-\theta$. An $h$-neighbor of the solution
$\bar{x}$ is defined by the current random keys, $h$, $v$, and the norm of
$v-\bar{x}$. The algorithm randomly selects solutions in the $h$-neighborhood,
one at a time.

Data: $x,h,n,\theta$

Result: The best neighbor found in the grid of solution $x$.

1 $\bar{x}\leftarrow x$;

2 $NumGridPoints\leftarrow n\cdot\lfloor 1.0/h\rfloor$;

3 $\mathit{NumPointsExamined}\leftarrow 0$;

4 while  _$\mathit{NumPointsExamined}\leq\mathit{NumGridPoints}$_ do

5 $\mathit{NumPointsExamined}\leftarrow\mathit{NumPointsExamined}+1$;

6 $y\leftarrow\mathit{h}\mathrm{Neighborhood}(\bar{x},h,n,\theta)$;

7 if  _$f(\mathcal{D}(y)) <f(\mathcal{D}(\bar{x}))$_ then

8 $\bar{x}\leftarrow y$;

9 $\mathit{NumPointsExamined}\leftarrow 0$;

10

11

return $\bar{x}$

Algorithm 4 GridSearch

Data: $\bar{x},h,n,\theta$

Result: A neighbor in the grid of solution $\bar{x}$.

1 for  _$i=1,...,n$_ do

2 if  _$\textrm{{UnifRand}}(0,1)\leq 0.5$_ then

3
$v_{i}\leftarrow\left\lceil\textrm{{UnifRand}}(1,(1.0-\bar{x}_{i})/h)\right\rceil$;

4

5 else

6
$v_{i}\leftarrow-1\cdot\left\lceil\textrm{{UnifRand}}(1,\bar{x}_{i}/h)\right\rceil$;

7

8

9for  _$i=1,...,n$_ do

10 if  _$\textrm{{UnifRand}}(0,1)\leq\theta$_ then

11 $x_{i}\leftarrow\bar{x}_{i}+h\cdot(v_{i}/\left\|v\right\|)$;

12

13

return $x$

Algorithm 5 $h$Neighborhood

The Nelder–Mead method, introduced by Nelder and Mead [38], is a numerical
technique for searching the minimum of an objective function within a
multidimensional space. This direct search approach, relying on function
comparisons, is commonly utilized in derivative-free nonlinear optimization.
This method begins with at least three solutions and can perform five moves:
reflection, expansion, inside contraction, outside contraction, and shrinking.
In this study, we always apply the Nelder-Mead Search with three solutions:
$x_{1}$, $x_{2}$, and $x_{3}$, where one is the current solution derived from
the constructive phase, while the others are randomly chosen from a pool of
elite solutions found during the search process. These solutions are ordered
by objective function value ($x_{1}$ is the best and $x_{3}$ is the worst).
Figure 1 shows an illustrative example of a simplex polyhedron and the five
moves.

Figure 1. Illustrative example of the simplex polyhedron and the five moves of
the Nelder-Mead Search. Source: based on the figure in [31]

Algorithm 6 presents pseudo-code for the Nelder-Mead Search adapted for
discrete optimization problems. We use an adapted Uniform Crossover (UX) [13]
to generate the moves. In the UX between solutions $a$ and $b$, for each
random key of the vector, we decide at random from which solution the random
key is copied to the corresponding random key of the new solution. We
incorporate a parameter $\mathit{factor}$ to determine whether, when using the
random key of solution $b$, to use the value of the random key itself
($\mathit{facto}r=1$) or an inverted value of the random key
($\mathit{factor}=-1$).

The algorithm begins with an initial simplex of three solutions
($x_{1},x_{2},x_{3}$). The simplex is sorted based on the objective function
values, and the simplex’s centroid ($x_{0}$) is computed between $x_{1}$ and
$x_{2}$ $(x_{0}=$UX$(x_{1},x_{2},1))$. The main loop iterates until a
termination condition is met. The algorithm performs a series of moves on the
simplex during each iteration to explore the search space. A reflection
solution $(x_{r}=$UX$(x_{0},x_{3},-1))$ is computed. If the objective function
value at $x_{r}$ is better than the current best solution ($x_{1}$), the
algorithm computes an expansion solution $(x_{e}=$UX$(x_{r},x_{0},-1))$. If
the objective function value at $x_{e}$ is better than at $x_{r}$, $x_{3}$ is
replaced by $x_{e}$; otherwise, $x_{3}$ is replaced by $x_{r}$. If neither the
reflection nor the expansion improves the solution, the algorithm contracts
towards the solution $x_{r}$ or $x_{3}$. For an outside contraction (when
$x_{r}$ is better than $x_{3}$), the contraction solution is
$x_{c}=$UX$(x_{r},x_{0},1)$. For an inside contraction (when $x_{r}$ is not
better than $x_{3}$), the contraction solution is $x_{c}=$UX$(x_{0},x_{3},1)$.
If the contraction step does not improve, the entire simplex is shrunk towards
the best solution $x_{1}$ $(x_{i}=$UX$(x_{1},x_{i},1),i=2,3)$. The algorithm
terminates when a specified termination condition is met. In this paper, we
use the same criterion of the Grid Search, the maximum number of iterations
equal $n\cdot\lfloor 1.0/h\rfloor$.

Data: $x_{1}$, $x_{2}$, $x_{3}$, $n$, $h$

Result: The best solution found in simplex $X$.

1 Initialize simplex: $X\leftarrow\\{x_{1},x_{2},x_{3}\\}$;

2 Sort simplex $X$ by objective function value;

3 Compute the simplex centroid $x_{0}\leftarrow$ UX$(x_{1},x_{2},1)$;

4 $\mathit{iter}\leftarrow 0$;

5 $numIter\leftarrow n\cdot\lfloor 1.0/h\rfloor$;

6 while  _$\mathit{iter} <\textit{numIter}$_ do

7 $shrink\leftarrow 0$;

8 $\mathit{iter}\leftarrow\mathit{iter}+1$;

9 Compute reflection solution $x_{r}\leftarrow$ UX$(x_{0},x_{3},-1)$;

10 if  _$f(\mathcal{D}(x_{r})) <f(\mathcal{D}(x_{1}))$_ then

11 Compute expansion solution $x_{e}\leftarrow$ UX$(x_{r},x_{0},-1)$;

12 if  _$f(\mathcal{D}(x_{e})) <f(\mathcal{D}(x_{r}))$_ then

13 $x_{3}\leftarrow x_{e}$;

14

15 else

16 $x_{3}\leftarrow x_{r}$;

17

18

19 else

20 if  _$f(\mathcal{D}(x_{r})) <f(\mathcal{D}(x_{2}))$_ then

21 $x_{3}\leftarrow x_{r}$;

22

23 else

24 if  _$f(\mathcal{D}(x_{r})) <f(\mathcal{D}(x_{3}))$_ then

25 Compute contraction solution $x_{c}\leftarrow$ UX$(x_{r},x_{0},1)$;

26 if  _$f(\mathcal{D}(x_{c})) <f(\mathcal{D}(x_{r}))$_ then

27 $x_{3}\leftarrow x_{c}$;

28 else

29 $shrink\leftarrow 1$;

30

31

32 else

33 Compute contraction solution $x_{c}\leftarrow$ UX$(x_{0},x_{3},1)$;

34 if  _$f(\mathcal{D}(x_{c})) <f(\mathcal{D}(x_{3}))$_ then

35 $x_{3}\leftarrow x_{c}$;

36 else

37 $shrink\leftarrow 1$;

38

39

40

41

42 if  _$shrink=1$_ then

43 Replace all solutions except the best $x_{1}$ with $x_{i}\leftarrow$
UX$(x_{1},x_{i},1),i=2,3$;

44

45 Sort simplex $X$ by objective function value;

46 Compute the simplex centroid $x_{0}\leftarrow$ UX$(x_{1},x_{2},1)$;

47

return $x_{1}$

Algorithm 6 NelderMeadSearch

Since a local optimal for a given local search heuristic may not be a local
optimal for another local search heuristic, using different heuristics usually
leads to a better improvement phase. However, the order in which local search
heuristics are explored must be determined. For this reason, we use a
variation of the Variable Neighborhood Descent (VND) [37] search that does not
have a predefined order for the heuristics to be applied. This search is
called Random Variable Neighborhood Descent (RVND) [45]. RVND randomly selects
the neighborhood heuristic order to be applied in each iteration. Algorithm 7
shows pseudo-code for the RVND.

The RVND algorithm starts by initializing a list of neighborhoods ($NL$)
representing different local search heuristics to explore the solution space.
It iteratively selects a neighborhood $\mathcal{N}^{i}$ randomly and searches
for the best neighboring solution ($x^{\prime}$) from that neighborhood. If
the new solution is better than the current solution ($x$), it updates $x$
with $x^{\prime}$ and restarts the exploration of neighborhoods. Otherwise,
the current neighborhood will be removed from the list. This process continues
until the list of neighborhoods is empty. Finally, the algorithm returns the
best solution found. RVND efficiently explores diverse solution spaces but can
also be applied to random-key spaces. Users can implement classic heuristics
for the problem on hand and encode the local optimal solution into the random-
key vector after the search process or implement random-key neighborhoods
independent of the problem on hand using the decoder to converge towards
better solutions.

Data: $x$

Result: The best solution in the neighborhoods.

1 Initialize the Neighborhood List ($NL$);

2 while  _$NL\neq 0$_ do

3 Choose a neighborhood $\mathcal{N}^{i}\in NL$ at random;

4 Find the best neighbor $x^{\prime}$ of $x\in\mathcal{N}^{i}$;

5 if  _$f(\mathcal{D}(x^{\prime})) <f(\mathcal{D}(x))$_ then

6 $x\leftarrow x^{\prime}$;

7 Restart $NL$;

8

9 else

10 Remove $\mathcal{N}^{i}$ from the $NL$ ;

11

12

return $x$

Algorithm 7 RandomVND

In this paper, we developed three additional problem-independent local search
heuristics that operate on the random-key vector. These heuristics encompass
distinct neighborhood structures: SwapRK, InvertRK, and FareyRK. The SwapRK
structure focuses on interchanging two positions within the random-key vector.
Conversely, the InvertRK structure perturbs the random-key values, inverting
the current value ($1-x_{i}$). Lastly, the FareyRK structure randomly changes
the current value of each random key with random values generated between each
two consecutive terms of the Farey sequence [39]. The Farey sequence of order
$\eta$ is the sequence of completely reduced fractions between 0 and 1, with
denominators less than or equal to $\eta$, arranged in increasing order. We
use the Farey sequence of order 7,

$F_{7}=\left\\{\frac{0}{1},\frac{1}{7},\frac{1}{6},\frac{1}{5},\frac{1}{4},\frac{2}{7},\frac{1}{3},\frac{2}{5},\frac{3}{7},\frac{1}{2},\frac{4}{7},\frac{3}{5},\frac{2}{3},\frac{5}{7},\frac{3}{4},\frac{4}{5},\frac{5}{6},\frac{6}{7},\frac{1}{1}\right\\}$

for our purposes. The random keys are explored in random order in each
iteration of these heuristics. The Grid Search, the Nelder-Mead Search, the
SwapRK, the InvertRK, and FareyRK heuristics are the neighborhoods of the RVND
algorithm, which systematically explores all neighbors of the current solution
using the first-improvement strategy, i.e., any improving move is immediately
applied.

##  4\. Applications

This section presents the applications of RK-GRASP to five combinatorial
optimization problems:

  1. (1)

Traveling salesman problem

  2. (2)

Hub location problem in graphs

  3. (3)

Steiner triple set covering problem

  4. (4)

Node capacitated graph partitioning problem

  5. (5)

Job sequencing and tool switching problem

We present a concise overview of each problem, including descriptions of the
solution encoding and a corresponding decoder. We aim to illustrate the wide
applicability of RK-GRASP as a problem-independent approach.

Initially, we directed our efforts toward applying the RK-GRASP to the
Traveling Salesman Problem (TSP). Our rationale behind this choice was
establishing an application in a widely recognized and classic problem domain.
By doing so, we aimed to facilitate the adoption of the RK-GRASP among new
users, leveraging the familiarity and accessibility of the TSP as a benchmark.
The framework of the RK-GRASP applied to solve the TSP is available in
www.github.com/antoniochaves19/RK-GRASP. Other problems reuse this framework
with the specific problem-dependent component (read data and decoder).

###  4.1. Traveling Salesman Problem.

####  4.1.1. Problem definition.

The Traveling Salesman Problem (TSP) [12] is a classic optimization problem in
which the objective is to find the shortest possible route that visits each
city exactly once and returns to the original city. Mathematically, given a
set of cities ($N$) and the distances between them, the goal is to determine
the optimal permutation of cities that minimizes the total distance traveled
by the salesman. The TSP is a well-studied combinatorial optimization problem
with logistics, transportation, and network design applications. We recommend
reading [48, 25] to review the literature on the TSP.

####  4.1.2. Solution encoding.

Each solution is encoded as a vector $x$ of random keys of length $n=|N|$,
where the i-th random key corresponds to the i-th city of $N$.

####  4.1.3. Solution decoder.

We use the best decoder developed for the TSP [6] to map the random key and
TSP solutions. The decoder starts by ordering the random-key vector in
ascending order. Then, using this sequence, the cheapest insertion
constructive heuristic is performed. A sub-route is created with the first
three cities, and from the fourth city in the sequence, the cheapest position
is calculated to insert the current city. The cost of inserting a city $k$
after city $j$ is calculated by $c_{k,j}=d_{j,k}+d_{k,j+1}-d_{j,j+1}$, where
$d_{i,j}$ is the distance or costs between the cities $i$ and $j$. In the end,
the cost of the TSP route is computed. Figure 2 shows an example of this
decoder with five cities.

Figure 2. Example of the TSP decoder with five cities.

###  4.2. Tree Hub Location Problem.

####  4.2.1. Problem definition.

The Tree Hub Location Problem (THLP) [9] consists of defining $p$ hubs in a
network that are connected by a tree with non-directional arcs. From this
tree, each non-hub node must be connected to a hub so that all the flow
between the nodes must use the hubs to travel. Each arc has a transportation
cost per unit, and we have a demand between each pair of nodes in the network.
The objective function involves two components: firstly, computing the
transportation costs for fulfilling all demands from non-hub nodes to a
designated hub, and secondly, factoring in the costs associated with the flow
within the hub tree. A discount factor is applied to achieve this, allowing
for a consideration of the overall cost dynamics involved in the
transportation process.

The literature on THLP encompasses both exact approaches, such as Mixed
Integer Linear Programming (MILP) models [8, 9] and Benders decomposition
[14], as well as heuristic methods like BRKGA [46]. Additionally, some
variants of THLP have been proposed [15, 30], expanding the scope of the
problem and its solution space. Notably, BRKGA has emerged as the state-of-
the-art method for THLP, presenting the best results reported in the
literature.

####  4.2.2. Solution encoding.

In this study, we adopt the encoder proposed by Pessoa et al. [46] for our
analysis. Consequently, the THLP solution is denoted by a random-key vector
composed of three distinct segments. As depicted in Figure 3, the first
segment of the vector, with a dimensionality of $|N|$, signifies the network
nodes. The second segment, sized at $|N|-p$, assigns non-hub nodes to their
corresponding hubs. Finally, the third segment, sized at $p(p-1)/2$,
represents all possible pairs of edges connecting the hubs.

####  4.2.3. Solution decoder.

Our decoder is also based on the decoder presented in [46]. It begins by
sorting the vector’s first segment to delineate hubs and non-hubs: the initial
$p$ positions denote hubs, while the subsequent $|N|-p$ positions signify non-
hubs. Subsequently, the second segment undertakes the task of assigning non-
hub nodes to their respective hubs. Achieving this involves partitioning the
interval [0, 1] into $p$ equal intervals, each associated with a specific hub
per the preceding ordering. Lastly, the decoder arranges the random keys
corresponding to inter-hub arcs in ascending order. The resulting tree
structure is then constructed using the Kruskal algorithm applied to the
sorted arcs. Figure 3 shows an example of the THLP decoder with $|N|=10$ and
$p=3$.

Figure 3. Example of the THLP decoder with ten points and three hubs.

###  4.3. Job Sequencing and Tool Switching Problem.

####  4.3.1. Problem definition.

The Job Sequencing and Tool Switching problem (SSP) [54] considers scheduling
a set of jobs ($T$) to be executed on a single machine. Each job $t$ requires
a particular set of tools ($F_{t}$) for processing, with the job being
executable only when this tool subset is available on the machine. However,
tool switches become critical when the machine is limited to holding a maximum
of C tools simultaneously, and C is less than the total number of tools
required for all jobs.

In the SSP, the objective is to sequence jobs to minimize the total number of
tool switches. A tool switch occurs when a tool is replaced with another in
the machine’s magazine. To put it simply, consider you have a toolbox with a
limited number of compartments. Each compartment can hold a different tool.
Now, if you need to use a tool that is not currently in the toolbox, you must
remove a tool from one of the compartments to make space for the new tool.

The SSP is a NP-hard problem ([10, 54]), in which we have two interconnected
decision processes: determine the optimal job processing sequence and identify
the necessary tools for each job’s processing to minimize tool switches, given
a fixed job processing sequence. For the last decision, Tang et al. [54]
introduce an optimal policy called KTNS (Keep Tool Needed Sooner). According
to this policy, when inserting a new tool and removing another from the
magazine, priority is given to retaining the tools required for upcoming job
tasks.

The literature review on the SSP highlights the research focus on heuristic
and metaheuristic strategies due to limited exact solution approaches.
Mathematical models are proposed by [54, 4, 10, 47, 32, 23, 11] with several
advances in terms of efficiency for solving SSP instances. The literature also
presents novel heuristics, such as Tabu Search [1] and Beam Search [58].
Furthermore, advancements like memetic algorithms by [2] and graph-based
representations [44] offered promising avenues for SSP solutions. [7] proposed
hybrid methods combining clustering search and BRKGA, while [35] presented a
hybrid genetic search, marking ongoing progress in SSP solution methodologies.
We recommend reading [11] and [35] for more in-depth literature reviews about
SSP.

####  4.3.2. Solution encoding.

We encode a solution of the SSP as a vector of $n$ random keys, where $n=|T|$
is the number of jobs and the $i$-th random key represents the $i$-th job.

####  4.3.3. Solution decoder.

Decoding the $n$ random keys of each solution into a sequence of jobs is
accomplished by sorting the jobs in ascending order of their corresponding
random key values. Figure 4 shows an example of the decoding process for the
SSP. In this example, there are five jobs. The sorted random keys correspond
to the jobs sequence (1, 4, 3, 2, 5). The objective function value of a
solution with the job sequence is obtained with the KTNS algorithm [54]. It
returns the number of tool switches needed to process the jobs in this
sequence.

Figure 4. Example of the SSP decoder with five jobs and a magazine with four
slots.

###  4.4. Steiner Triple Covering Problem.

####  4.4.1. Problem definition.

Fulkerson et al. [22] introduced the Steiner Triple Covering Problem as a
computationally challenging set covering problem category. A Steiner triple
system $B$ on a set $X$ with $N$ elements is a collection of triples (subsets
of $X$ with three elements), such that any pair of elements $x$ and $y$ from
$X$ is contained in exactly one triple of $B$. For example, the Fano plane
(Figure 5) exemplifies a Steiner triple system on a set $X$ with seven
elements, where the triples represent the seven lines, each containing three
points, ensuring that every pair of points is part of a unique line.

A Steiner triple system can be represented by a binary incidence matrix $A$,
which has a column for every element in $X$ and a row for each triple in $B$.
In matrix $A$, the entry $A(i,j)$ is equal to 1 if and only if the $j$-th
element is part of the $i$-th triple. Consequently, every row $i$ in matrix
$A$ will have precisely three entries where $A(i,j)$ equals 1. Furthermore,
for every pair of columns $j$ and $k$ there is exactly one row $i$ for which
$A_{i,j}=A_{i,k}=1$.

(a) Fano Plane with $X=\\{1,\ldots,7\\}$, and
$B=\\{\\{1,2,3\\},\\{2,4,6\\},\\{1,4,5\\},\\{1,6,7\\},$
$\\{3,4,7\\},\\{2,5,7\\}\\}$.

$\begin{bmatrix}1&1&1&0&0&0&0\\\ 0&1&0&1&0&1&0\\\ 1&0&0&1&1&0&0\\\
1&0&0&0&0&1&1\\\ 0&0&1&1&0&0&1\\\ 0&1&0&0&1&0&1\\\ \end{bmatrix}$

(b) Incidence matrix of Fano Plane.

Figure 5. Example of STCP input data.

Several heuristic methods for addressing the Steiner triple covering problem
have been suggested. Feo and Resende [18] introduced a GRASP heuristic,
through which they identified a cover of size 61 for $A_{81}$, a result later
confirmed as optimal by Mannino and Sassano [34]. Karmarkar et al. [29]
developed an interior point method that achieved a cover of size 105 for
$A_{135}$. In the same study, they also applied a GRASP heuristic, finding an
improved cover of size 104, a result matched by Mannino and Sassano [34]. In
1998, Odijk and van Maaren [40] produced a cover of size 103. This result was
recently validated as optimal by Ostrowski et al. [42, 43] and Östergård and
Vaskelainen [59]. Moreover, Östergård and Vaskelainen [59] identified all
possible solutions for $A_{135}$, demonstrating their isomorphism.

####  4.4.2. Solution encoding.

We encoded a solution of the STCP as a random-key vector of length $n=|X|$,
where the $i$-th random key corresponds to the $i$-th element of $X$ (a column
of $B$).

####  4.4.3. Solution decoder.

The decoder process for the STCP is designed to create an initial solution and
improve it by eliminating redundancy. The random-key vector is sorted and a
loop iterates through the column sequence, covering the necessary rows and
tracking the number of used columns and coverage redundancy. If a column is
used to cover at least one row, it is marked as used. After ensuring all rows
are covered, the function attempts to reduce redundancy by checking each used
column to see if it is necessary for maintaining coverage. Redundant columns
are removed, and the coverage vector is updated accordingly. Finally, we
calculates the objective function value, representing the number of columns
used. Figure 6 shows an example of the STCP decoder considering the Fano
Plane.

Figure 6. Example of the STCP decoder with seven elements and the incident
matrix of Fano Plane.

###  4.5. Node Capacitated Graph Partitioning Problem.

####  4.5.1. Problem definition.

The node capacitated graph partitioning problem is also known as the handover
minimization problem described as follows. Let $B$ denotes the set of base
stations, with $T_{b}$ representing the total traffic handled by base station
$b\in B$ and its connected transceivers. Let $N$ be the set of RNCs (Radio
Network Controllers), where each $\mathrm{RNC}$ $r\in N$ can handle a maximum
traffic capacity of $C_{r}$. Furthermore, define $H_{b_{1},b_{2}}$ as the
total number of handovers between base stations $b_{1}$ and $b_{2}$
($b_{1},b_{2}\in B,b_{1}\neq b_{2}$). Note that $H_{b_{1},b_{2}}$ and
$H_{b_{2},b_{1}}$ may differ.

The objective of the handover minimization problem is to allocate each base
station in $B$ to a specific RNC in $N$. Let $\rho_{b}$ denote the index of
the RNC to which base station $b\in B$ is assigned, and let $\Psi_{r}$ denote
the indices of base stations assigned to RNC $r\in N$. The assignments must
satisfy the capacity constraints of each RNC, meaning that for all $r\in N$:

$\sum_{b\in\Psi_{r}}T_{b}\leq C_{r}.$

Among all feasible assignments, we seek to minimize the total handovers
between base stations assigned to different RNCs:

$\sum_{b_{1},b_{2}\in B\mid\rho_{b_{1}}\neq\rho_{b_{2}}}H_{b_{1},b_{2}}.$

Ferreira et al. [21] were first to study this problem, proposing strong valid
inequalities for a branch and cut algorithm. They tested their algorithm on
three applications of partitioning: compiler design, finite element
computations associated with meshes, and design of electronic circuits. The
largest instances they were able to solve to optimality had 61 nodes and 187
edges for compiler design, 48 nodes and 81 edges for design of electronic
circuits, and 274 nodes and 469 edges for finite element computations on a
mesh. Mehrotra and Trick [36] proposed a branch and price algorithm for the
node capacitated graph partitioning problem. They tested their algorithm on
instances having 30-61 nodes and 47-187 edges. The largest instances they were
able to solve to optimality had 61 nodes and 187 edges. Deng and Bard [16]
proposed a reactive GRASP with path-relinking post-processing for the node
capacitated graph partitioning problem. They tested their heuristic on
instances varying in size from 30 to 82 nodes and 65 to 540 edges. They used
CPLEX 11 to solve to optimality instances having 30 nodes and most of the
instances having 40 nodes. They showed that their heuristic was able to find
solutions having the same objective function value than those found by CPLEX
on most instances, but in less time than CPLEX. On larger instances, CPLEX
failed to prove optimality and found worse solutions than those found by the
reactive GRASP with path-relinking for all but one instance where they were
tied. They also compared their reactive GRASP heuristic with the combinatorial
algorithm of Mehrotra and Trick [36]. Deng, Y., Bard, J.F. [16] showing that
their heuristic matched the optimal solutions found by the combinatorial
algorithm on all but one instance where it found a suboptimal solution.

####  4.5.2. Solution encoding.

Solutions of the NCGPP are encoded as a vector of $n=|B|+1$ random keys. The
first $|B|$ positions represent each base station, while the last random key
indicates the number of base stations assigned initially. Once sorted, the
$|B|$ random keys will dictate the order in which base stations are assigned
to the RNCs.

####  4.5.3. Solution decoder.

To decode a random-key vector and produce an assignment of base stations to
RNCs, the following steps are computed: sort the random-key vector and then
group base stations into RNCs, considering capacity constraints and costs. The
assignment process starts by ensuring the first $\\#N=\left\lceil
x_{n}\cdot|N|\right\rceil$ base stations are assigned to separate RNCs. For
each subsequent base station ($|B|-\\#N$), it searches for the best RNC that
can accommodate the base station without exceeding capacity, evaluating the
insertion cost based on the sum of handovers between the current base station
and base stations that are already in the RNCs. If a suitable RNC is found,
the base station is added, and the RNC’s capacity is updated. A penalty is
added to the solution if no RNC can accommodate the base station. Finally, the
decoder calculates the objective function value by summing the handover
between base stations of different RNCs and the penalties for ungrouped base
stations.

Figure 7 shows an example of the NCGPP decoder with six base stations
($|B|=6$) and two RNCs ($|N|=2$). In this example, $\\#N=\left\lceil 0.7\cdot
2\right\rceil=2$ and the ordered random-key vector provides the sequence of
base stations: $2,4,5,3,6,1$. From these $\\#N$ and sequence, base stations
$2$ and $4$ are allocated to separate RNCs. Base station $5$ is then allocated
to RNC 2, as base station $5$ has 191 handovers with base station $4$ and only
116 with base station $2$. Base station $3$ is allocated to RNC 1, as it has
no handover with base stations $4$ and $5$, only with base station $2$. Base
station $6$ is allocated to RNC 2, which has 307 handovers (157 with base
station $4$ and 150 with base station $5$) against 13 for RNC 1. Finally, base
station $1$ is allocated to RNC 1 because it also has handovers only with base
stations $2$ and $3$.

Figure 7. Example of the NCGPP decoder with six base stations and two RNCs.

##  5\. Computational experiments

The RK-GRASP framework was coded in C++ and compiled with GCC. The computer
used in all experiments was a Dual Xenon Silver 4114 20c/40t 2.2Ghz processor
with 96GB of DDR4 RAM and running CentOS 8.0 x64. We tested $20$ instances of
the TSP, $36$ instances of the THLP, five instances of the STCP, $20$
instances of the NCGPP, and $40$ instances of the SSP. Each algorithm was run
$5$ times for each instance with different seeds and a CPU time limit
($run_{time}$). The CPU time limit was configured in seconds, taking into
account the size of each tested instance. For the TSP, the CPU time
corresponds to the number of cities in the instance. For the THLP, it is based
on the number of nodes in the network; for the STCP, it corresponds to the
number of columns; for the SSP, it is determined by the number of tasks; and
for the NCGPP, it is based on the number of base stations. This stopping
criterion was selected to allow a fairer comparison between different methods
using the same architecture. However, it can be replaced by other criteria,
such as the maximum number of objective function evaluations. The RK-GRASP has
three parameters to be configured: $h_{s}$, $h_{e}$, and $\theta$. These
parameters are tuned based on previous computational tests as $h_{s}=0.125$,
$h_{e}=0.00098$, and $\theta=$ UnifRand$[0.2,0.4]$. The simulated annealing
acceptance criterion is $e^{-\Delta/(run_{time}-time)}$, where $\Delta$ is the
difference between $f(\mathcal{D}(x^{\prime\prime}))$ and $f(\mathcal{D}(x))$,
and $time$ is the current CPU time. The framework, instances and detailed
results are available at https://github.com/antoniochaves19/RK-GRASP.

Table 1 summarizes the results of the computational experiments of the three
versions of the RK-GRASP applied to solve five combinatorial optimization
problems. The entries in this table are the problem on hand, the instance set,
the number of tested instances, the best-known solution (BKS), the average
solution over five runs (Avrg), the gap between the best solution found and
the BKS (Gap), the average relative percentage deviation (ARPD) considering
solutions found in each run and the BKS
($RPD(\%)=\frac{(f(\mathcal{D}(x^{*}))-BKS)}{BKS}\times 100$), and the average
running time to find the best solution. The values in boldface show the best
average objective function value for each instance set.

Table 1. RK-GRASP: Computational results for different combinatorial
optimization problems.

|  |  |  | RK-GRASP Grid Search | RK-GRASP NelderMead Search | RK-GRASP RVND Search
---|---|---|---|---|---|---
Problem Class | Instance set | # inst | BKS | Avrg | Gap | ARPD | Time (s) | Avrg | Gap | ARPD | Time (s) | Avrg | Gap | ARPD | Time (s)
TSP | $|N|<100$ | 5 | 27639.2 | 27639.2 | 0.00 | 0.00 | 5.09 | 27639.2 | 0.00 | 0.00 | 0.74 | 27639.2 | 0.00 | 0.00 | 0.34
| $|N|<200$ | 5 | 57468.8 | 57469.6 | 0.00 | 0.01 | 25.18 | 57468.8 | 0.00 | 0.01 | 16.43 | 57468.8 | 0.00 | 0.00 | 10.30
| $|N|<400$ | 5 | 60110.8 | 60148.6 | 0.06 | 0.34 | 104.56 | 60151.2 | 0.07 | 0.19 | 90.75 | 60121.4 | 0.02 | 0.10 | 64.48
| $|N|<1000$ | 5 | 42340.3 | 42649.5 | 0.73 | 1.03 | 311.01 | 42572.8 | 0.55 | 0.75 | 356.19 | 42398.8 | 0.14 | 0.31 | 319.95
THLP | cab20 | 9 | 4397.9 | 4519.7 | 2.77 | 7.39 | 9.49 | 4436.2 | 0.87 | 2.84 | 8.76 | 4396.4 | -0.03 | 0.00 | 1.98
| cab25 | 9 | 7187.6 | 7540.8 | 4.91 | 9.64 | 11.26 | 7268.3 | 1.12 | 4.25 | 12.31 | 7177.7 | -0.14 | -0.06 | 4.20
| ap20 | 9 | 60397.5 | 62322.5 | 3.19 | 5.87 | 10.51 | 60928.9 | 0.88 | 2.85 | 7.02 | 60383.7 | -0.02 | -0.03 | 2.28
| ap25 | 9 | 61495.0 | 64705.5 | 5.22 | 9.31 | 13.53 | 63113.4 | 2.63 | 5.53 | 9.36 | 61443.8 | -0.08 | 0.01 | 5.20
STCP | stn81 | 1 | 61 | 61.0 | 0.00 | 0.00 | 0.746 | 61.0 | 0.00 | 0.00 | 0.072 | 61.0 | 0.00 | 0.00 | 0.281
| stn135 | 1 | 103 | 104.0 | 0.97 | 1.17 | 20.813 | 104.0 | 0.97 | 0.97 | 51.841 | 103.0 | 0.00 | 0.78 | 1.876
| stn243 | 1 | 198 | 198.0 | 0.00 | 1.92 | 140.494 | 198.0 | 0.00 | 0.91 | 57.495 | 198.0 | 0.00 | 0.00 | 20.125
| stn405 | 1 | 335 | 344.0 | 2.69 | 2.93 | 95.396 | 343.0 | 2.39 | 2.63 | 299.486 | 339.0 | 1.19 | 1.55 | 184.676
| stn729 | 1 | 617 | 648.0 | 5.02 | 5.25 | 467.249 | 648.0 | 5.02 | 5.19 | 369.002 | 640.0 | 3.73 | 3.95 | 637.719
NCGPP | 20/5 | 5 | 381.6 | 381.6 | 0.00 | 0.00 | 0.11 | 381.6 | 0.00 | 0.00 | 0.04 | 381.6 | 0.00 | 0.00 | 0.03
| 40/15 | 5 | 5940.0 | 5940.0 | 0.00 | 0.00 | 6.27 | 5940.0 | 0.00 | 0.01 | 1.62 | 5940.0 | 0.00 | 0.00 | 0.12
| 100/25 | 5 | 35981.2 | 36409.6 | 1.19 | 1.85 | 45.83 | 35950.8 | -0.08 | 0.19 | 45.82 | 35943.6 | -0.10 | -0.09 | 7.94
| 200/50 | 5 | 220969.6 | 220472.0 | -0.23 | 0.42 | 104.10 | 217267.6 | -1.68 | -1.34 | 165.30 | 215788.8 | -2.34 | -2.29 | 110.90
SSP | c10-10-4 | 10 | 9.1 | 9.1 | 0.00 | 0.00 | 0.07 | 9.1 | 0.00 | 0.00 | 0.05 | 9.1 | 0.00 | 0.00 | 0.01
| c15-20-6 | 10 | 20.6 | 20.7 | 0.49 | 0.53 | 1.11 | 20.7 | 0.49 | 1.01 | 1.57 | 20.6 | 0.00 | 0.00 | 0.16
| c30-40-15 | 10 | 92.7 | 95.2 | 2.70 | 4.38 | 12.80 | 93.7 | 1.08 | 2.63 | 10.84 | 91.7 | -1.08 | -0.41 | 7.35
| c40-60-20 | 10 | 182.7 | 189.8 | 3.89 | 5.68 | 19.49 | 183.9 | 0.66 | 2.35 | 20.74 | 179.8 | -1.59 | -0.76 | 15.41

The three versions of the RK-GRASP generated effective solutions for the
problems on hand. The RK-GRASP variant with Grid Search found high-quality
solutions for the TSP, STCP, and NCGPP, improving the Best Known Solution
(BKS) of the literature for the large-scale instances of the NCGPP. The RK-
GRASP with Nelder-Mead Search found better results than RK-GRASP with Grid
Search, being more competitive also for the THLP and SSP in terms of quality
solution. However, the RK-GRASP variant with RVND Search emerged as the most
effective, delivering better solutions for the five combinatorial optimization
problems. The results of this version are competitive for the TSP and STCP and
for the THLP, NCGPP, and SSP this variant outperformed state-of-the-art
methods, presenting several negative gaps. We can observe that the RK-GRASP
converges to high-quality solutions in an efficient CPU time. In summary, all
versions of RK-GRASP demonstrated efficiency in both solution quality and
computational time, showing their effectiveness across diverse problems and
instances.

To assess the efficacy of incorporating the random-key concept into GRASP, we
conducted experiments employing the same decoders with the BRKGA
metaheuristic. Additionally, we employed a Multi-Start algorithm, in which
random-key vectors are randomly generated and decoded, to highlight the
significance of the RK-GRASP and BRKGA methodologies in achieving high-quality
solutions. The results are presented in Tables 2 and 3, with detailed
information in the GitHub repository. Table 3 presents the overall results for
the five combinatorial optimization problems. For each problem and algorithm,
the median, the mean, and the average gap values across all runs and instances
are shown. Additionally, the table includes the number of instances tested for
each problem, the number of instances where the algorithm matched the best-
known solution from the literature (#BKS), and the number of instances where
the algorithm found a new best-known solution (#newBKS). We also conducted a
Wilcoxon signed-rank test (WSR) [51], a non-parametric statistical hypothesis
test, to compare the rank sums between each pair of methods (Table 2). In
interpretation, a $p$-value $<0.05$ with alternative = “less” implies that the
first group tends to produce smaller values than the second. In contrast, a
$p$-value $<0.05$ with alternative = “greater” indicates that the first group
tends to produce larger values. The values in boldface show the cases with
statistical significance ($p$-value $<0.05$).

Table 2. Wilcoxon test comparing the different algorithms in terms of
$p$-value.

Alternate hypothesis: solution set A is less than solution set B (A x B)
---
Problem Class | TSP | THLP | STCP | NCGPP | SSP
RK-GRASP Grid Search x RK-GRASP NelderMead Search | 9.99E-01 | 1.00E+00 | 9.96E-01 | 1.00E+00 | 1.00E+00
RK-GRASP Grid Search x RK-GRASP RVND Search | 1.00E+00 | 1.00E+00 | 1.00E+00 | 1.00E+00 | 1.00E+00
RK-GRASP Grid Search x BRKGA | 1.00E+00 | 1.00E+00 | 9.98E-01 | 1.00E+00 | 1.00E+00
RK-GRASP Grid Search x MultiStart | 8.46E-14 | 2.20E-16 | 1.22E-01 | 8.46E-14 | 2.20E-16
RK-GRASP NelderMead Search x RK-GRASP RVND Search | 1.00E+00 | 1.00E+00 | 9.99E-01 | 1.00E+00 | 1.00E+00
RK-GRASP NelderMead Search x BRKGA | 1.00E+00 | 1.00E+00 | 9.86E-01 | 1.00E+00 | 1.00E+00
RK-GRASP NelderMead Search x MultiStart | 8.46E-14 | 2.20E-16 | 3.60E-03 | 8.46E-14 | 2.20E-16
RK-GRASP RVND Search x BRKGA | 3.03E-02 | 2.20E-16 | 1.42E-02 | 2.48E-16 | 3.91E-09
RK-GRASP RVND Search x MultiStart | 8.46E-14 | 2.20E-16 | 3.15E-04 | 8.46E-14 | 2.20E-16
Alternate hypothesis: solution set A is greater than solution set B (A x B)
Problem Class | TSP | THLP | STCP | NCGPP | SSP
RK-GRASP Grid Search x RK-GRASP NelderMead Search | 9.24E-04 | 2.20E-16 | 6.44E-03 | 3.90E-10 | 8.75E-16
RK-GRASP Grid Search x RK-GRASP RVND Search | 4.66E-08 | 2.20E-16 | 2.22E-04 | 3.90E-10 | 2.20E-16
RK-GRASP Grid Search x BRKGA | 1.01E-07 | 2.20E-16 | 2.21E-03 | 3.90E-10 | 2.20E-16
RK-GRASP Grid Search x MultiStart | 1.00E+00 | 1.00E+00 | 8.98E-01 | 1.00E+00 | 1.00E+00
RK-GRASP NelderMead Search x RK-GRASP RVND Search | 5.34E-06 | 2.20E-16 | 7.62E-04 | 1.85E-08 | 2.20E-16
RK-GRASP NelderMead Search x BRKGA | 9.96E-06 | 4.16E-04 | 1.78E-02 | 3.34E-06 | 2.91E-13
RK-GRASP NelderMead Search x MultiStart | 1.00E+00 | 1.00E+00 | 9.97E-01 | 1.00E+00 | 1.00E+00
RK-GRASP RVND Search x BRKGA | 9.71E-01 | 1.00E+00 | 9.90E-01 | 1.00E+00 | 1.00E+00
RK-GRASP RVND Search x MultiStart | 1.00E+00 | 1.00E+00 | 1.00E+00 | 1.00E+00 | 1.00E+00

Table 3. Comparison of computational results for the three variants of RK-
GRASP, BRKGA, and Multi-Start.

|  | RK-GRASP | RK-GRASP | RK-GRASP | BRKGA | Multi-Start
---|---|---|---|---|---|---
|  | Grid Search | NelderMead Search | RVND Search |  |
TSP | median | 38974.5 | 38872.5 | 38694.0 | 38727.0 | 39545.0
| mean | 45248.0 | 45214.8 | 45147.6 | 45166.5 | 45810.9
| avrg gap | 0.25 | 0.24 | 0.07 | 0.12 | 1.66
| #inst | 20 | 20 | 20 | 20 | 20
| #BKS | 13 | 14 | 14 | 14 | 6
| #newBKS | 0 | 0 | 0 | 0 | 0
THLP | median | 24811.1 | 22965.5 | 22610.4 | 22962.4 | 31175.4
| mean | 35575.4 | 34583.1 | 33363.1 | 34255.5 | 42441.3
| avrg gap | 4.94 | 1.51 | -0.06 | 1.39 | 28.72
| #inst | 36 | 36 | 36 | 36 | 36
| #BKS | 7 | 12 | 36 | 10 | 0
| #newBKS | 0 | 0 | 5 | 1 | 0
STCP | median | 203.0 | 198.0 | 198.0 | 198.0 | 202.0
| mean | 272.2 | 271.5 | 268.8 | 270.7 | 272.4
| avrg gap | 1.74 | 1.68 | 0.98 | 1.40 | 1.83
| #inst | 5 | 5 | 5 | 5 | 5
| #BKS | 2 | 2 | 3 | 3 | 2
| #newBKS | 0 | 0 | 0 | 0 | 0
NCGPP | median | 20699.0 | 20666.0 | 20666.0 | 20666.0 | 22756.0
| mean | 66228.0 | 65109.5 | 64552.6 | 64848.4 | 70850.2
| avrg gap | 0.23 | -0.44 | -0.62 | -0.53 | 4.62
| #inst | 20 | 20 | 20 | 20 | 20
| #BKS | 14 | 19 | 20 | 20 | 7
| #newBKS | 0 | 1 | 6 | 1 | 0
SSP | median | 52.0 | 52.5 | 51.0 | 51.0 | 64.0
| mean | 79.9 | 78.0 | 75.8 | 76.7 | 96.4
| avrg gap | 1.75 | 0.57 | -0.67 | -0.33 | 14.95
| #inst | 40 | 40 | 40 | 40 | 40
| #BKS | 19 | 24 | 39 | 37 | 11
| #newBKS | 0 | 1 | 14 | 7 | 0

We observe that RK-GRASP with RVND emerged as the most efficient method for
the five combinatorial optimization problems, with statistical differences
evident among the other two variants of RK-GRASP, BRKGA, and Multi-Start.
Additionally, BRKGA demonstrated statistical superiority over RK-GRASP with
Grid Search and Nelder-Mead, becoming the second best choice. The RK-GRASP
with Nelder-Mead performed better than RK-GRASP with Grid Search. Overall, all
metaheuristics outperformed the Multi-Start approach. The RK-GRASP with RVND
improved the BKS in 5 instances of THLP (14%), 6 instances of NCGPP (30%), and
14 instances of SSP (35%). This RK-GRASP version also achieved a solution
equal to or better than the BKS in 92.5% of the tested instances.

The performance profile [17] was used to evaluate and compare the performance
of three variants of RK-GRASP, BRKGA, and Multi-Start methods based on CPU
time. The performance ratio is calculated by:

$r_{i,a}=\frac{t_{i,a}}{min\\{t_{i,a}:\forall a\\}}$

where $t_{i,a}$ is the average CPU time required to find the best solution to
instance $i$ using algorithm $a$. We define the level of accuracy as 1%. Each
instance where an algorithm had a gap greater than 1% was assigned an infinite
CPU time ($\mathit{INF}$) to indicate failure in meeting the convergence test.

The cumulative distribution function $\rho_{a}(\tau)$ represents the
probability that algorithm $a$ achieves a performance ratio $r_{i,a}$ within a
factor $\tau$ of the best possible ratio:

$\rho_{a}(\tau)=\frac{\left|\\{i\in I:r_{i,a}\leq\tau\\}\right|}{|I|}.$

Algorithms with higher $\rho_{a}(\tau)$ values are preferred.

Figure 8 presents a $log_{2}$-scaled view of the performance profiles of the
three variants of RK-GRASP, BRKGA, and Multi-Start methods. Again, we observe
the dominance of the RK-GRASP with RVND in terms of CPU time. This variant of
the RK-GRASP was the most efficient method in 97% of THLP instances, 50% of
NCGPP instances, and 83% of SSP instances, while BRKGA was the most efficient
method in 60% of TSP instances. These two methods always dominate the other
methods, except in the case of STCP, where no visible dominance is evident.
Considering a performance rate of $\tau=4$ (x-axis = 2), we note that BRKGA
can find the target solution for 95% of TSP, 15% of THLP, 40% of STCP, 80% of
NCGPP, and 73% of SSP instances, while RK-GRASP RVND solved 95% of TSP, 100%
of THLP, 40% of STCP, 75% of NCGPP, and 100% of SSP instances.

(a) TSP.

(b) THLP.

(c) STCP.

(d) NCGPP.

(e) SSP.

Figure 8. Performance profile of runtime for the RK-GRASP versions to find the
best solution.

##  6\. Concluding remarks

This paper introduces RK-GRASP, a problem-independent greedy randomized
adaptive search procedure (GRASP) for combinatorial optimization, leveraging
the random-key optimizer (RKO) paradigm. RK-GRASP distinguishes itself from
standard GRASP by operating in the continuous unit hypercube and employing a
decoder to map solutions to the discrete problem space. This design allows for
a problem-independent implementation, requiring only the development of a
specific decoder for each problem.

Three local search methods, namely Grid Search, Nelder-Mead Search, and Random
Variable Neighborhood Descent (RVND), are integrated into RK-GRASP. The paper
demonstrates the applicability of RK-GRASP across five NP-hard combinatorial
optimization problems: Traveling Salesman Problem (TSP), Tree of Hubs Location
Problem (THLP), Steiner Triple Covering Problem (STCP), Node Capacitated Graph
Partitioning Problem (NCGPP), and Job Sequencing and Tool Switching Problem
(SSP).

Computational experiments show that all three versions of RK-GRASP generate
effective solutions. RK-GRASP with Nelder-Mead Search performs well for TSP
and NCGPP instances, while RK-GRASP with RVND Search emerges as the most
effective for the five combinatorial optimization problems. Notably, RK-GRASP
with RVND Search outperforms state-of-the-art methods for THLP, SSP, and
NCGPP, achieving new best-known solutions.

The paper also compares RK-GRASP with BRKGA and a Multi-Start approach using
the same decoders. Again, RK-GRASP with RVND demonstrates significant
statistical superiority and is more effective for the problems on hand.

In conclusion, this paper presents RK-GRASP as a versatile and efficient
problem-independent metaheuristic for combinatorial optimization. The RVND was
the best choice of local search technique within RK-GRASP, considering the
specific problems and tested instances. The results highlight the
effectiveness of the random-key concept within the GRASP framework, offering a
promising approach for tackling complex optimization challenges across diverse
problem domains.

As we did for BRKGA in [3, 41, 55], we envision a future where an API for RK-
GRASP will facilitate the use of this metaheuristic to solve a wide range of
combinatorial optimization problems.

Acknowledgments

Antônio A. Chaves was supported by São Paulo Research Foundation (FAPESP)
under grants 2018/15417-8 and 2022/05803-3, and Conselho Nacional de
Desenvolvimento Científico e Tecnológico (CNPq) under grants 312747/2021-7 and
405702/2021-3. Ricardo M. A. Silva was partially supported by the Foundation
for the Support of Development of the Federal University of Pernambuco, Brazil
(FADE), the Office for Research and Graduate Studies of the Federal University
of Pernambuco (PROPESQI),and the Foundation for Support of Science and
Technology of the State of Pernambuco (FACEPE).

## References

  * Al-Fawzan and Al-Sultan [2003] M. A. Al-Fawzan and K. S. Al-Sultan.  A tabu search based algorithm for minimizing the number of tool switches on a flexible machine.  _Computers & Industrial Engineering_, 44(1):35–47, 2003.
  * Amaya et al. [2020] J. E. Amaya, C. Cotta, A. J. Fernández-Leiva, and P. García-Sánchez.  Deep memetic models for combinatorial optimization problems: application to the tool switching problem.  _Memetic Computing_ , 12(1):3–22, 2020.
  * Andrade et al. [2021] C. E. Andrade, R. F. Toso, J. F. Gonçalves, and M. G. C. Resende.  The multi-parent biased random-key genetic algorithm with implicit path-relinking and its real-world applications.  _European Journal of Operational Research_ , 289(1):17–30, 2021.
  * Bard [1988] J. F. Bard.  A heuristic for minimizing the number of tool switches on a flexible machine.  _IIE Transactions_ , 20(4):382–391, 1988.
  * Bean [1994] J. C. Bean.  Genetic algorithms and random keys for sequencing and optimization.  _ORSA Journal on Computing_ , 6(2):154–160, 1994\.
  * Chaves and Lorena [2021] A. A. Chaves and L. H. N. Lorena.  An adaptive and near parameter-free BRKGA using $Q$-learning method.  In _2021 IEEE Congress on Evolutionary Computation (CEC)_ , pages 2331–2338, 2021.
  * Chaves et al. [2016] A. A. Chaves, L. A. N. Lorena, E. L. F. Senne, and M. G. C. Resende.  Hybrid method with cs and brkga applied to the minimization of tool switches problem.  _Computers & Operations Research_, 67:174–183, 2016.
  * Contreras et al. [2009] I. Contreras, E. Fernández, and A. Marín.  Tight bounds from a path based formulation for the tree of hub location problem.  _Computers & Operations Research_, 36(12):3117–3127, 2009.
  * Contreras et al. [2010] I. Contreras, E. Fernández, and A. Marín.  The tree of hubs location problem.  _European Journal of Operational Research_ , 202(2):390–400, 2010.
  * Crama et al. [1994] Y. Crama, A. G. Oerlemans, and F. C. R. Spieksma.  Minimizing the number of tool switches on a flexible machine.  _International Journal of Flexible Manufacturing Systems_ , 6:33–54, 1994.
  * da Silva et al. [2021] T. T. da Silva, A. A. Chaves, and H. H. Yanasse.  A new multicommodity flow model for the job sequencing and tool switching problem.  _International Journal of Production Research_ , 59(12):3617–3632, 2021.
  * Dantzig et al. [1954] G. B. Dantzig, R. Fulkerson, and S. Johnson.  Solution of a large-scale traveling-salesman problem.  _Math. Programm. Study_ , 2:393–410, 1954.
  * Davis [1991] L. D. Davis, editor.  _Handbook of Genetic Algorithms_.  Van Nostrand Reinhold, 1991.
  * de Sá et al. [2013] E. M. de Sá, R. S. de Camargo, and G. de Miranda.  An improved Benders decomposition algorithm for the tree of hubs location problem.  _European Journal of Operational Research_ , 226(2):185–202, 2013.
  * de Sá et al. [2015] E. M. de Sá, I. Contreras, and J.-F. Cordeau.  Exact and heuristic algorithms for the design of hub networks with multiple lines.  _European Journal of Operational Research_ , 246(1):186–198, 2015.
  * Deng and Bard [2011] Y. Deng and J.F. Bard.  A reactive grasp with path relinking for capacitated clustering.  _Journal of Heuristics_ , 17:119–152, 2011.
  * Dolan and Moré [2002] E. D. Dolan and J. J. Moré.  Benchmarking optimization software with performance profiles.  _Mathematical programming_ , 91:201–213, 2002.
  * Feo and Resende [1989a] T. A. Feo and M. G. C. Resende.  A probabilistic heuristic for a computationally difficult set covering problem.  _Operations Research Letters_ , 8(2):67–71, 1989a.
  * Feo and Resende [1989b] T. A. Feo and M. G. C. Resende.  A probabilistic heuristic for a computationally difficult set covering problem.  _Operations Research Letters_ , 8:67–71, 1989b.
  * Feo and Resende [1995] T. A. Feo and M. G. C. Resende.  Greedy randomized adaptive search procedures.  _Journal of Global Optimization_ , 6(1):109–133, 1995.
  * Ferreira et al. [1998] C.E. Ferreira, A. Martin, C.C. de Souza, R. Weismantel, and L.A. Wolsey.  The node capacitated graph partitioning problem: a computational study.  _Math. Program._ , 81:229–256, 1998.
  * Fulkerson et al. [1974] D. R. Fulkerson, G. L. Nemhauser, and L. E. Trotter Jr.  Two computationally difficult set covering problems that arise in computing the 1-width of incidence matrices of steiner triple systems.  _Journal of the Operations Research Society of America_ , 2:72–81, 1974.
  * Ghiani et al. [2010] G. Ghiani, A. Grieco, and E. Guerriero.  Solving the job sequencing and tool switching problem as a nonlinear least cost Hamiltonian cycle problem.  _Networks: An International Journal_ , 55(4):379–385, 2010.
  * Gonçalves and Resende [2011] J. F. Gonçalves and M. G. C. Resende.  Biased random-key genetic algorithms for combinatorial optimization.  _Journal of Heuristics_ , 17(1):487–525, 2011\.
  * Gutin and Punnen [2006] G. Gutin and A. P. Punnen.  _The traveling salesman problem and its variations_ , volume 12.  Springer Science & Business Media, 2006.
  * Hirsch et al. [2007] M. J. Hirsch, C. N. Meneses, P. M. Pardalos, and M. G. C. Resende.  Global optimization by continuous GRASP.  _Optimization Letters_ , 1:201–212, 2007.
  * Hirsch et al. [2010] M. J. Hirsch, P. M. Pardalos, and M. G. C. Resende.  Speeding up continuous GRASP.  _European J. of Operational Research_ , 205:507–521, 2010\.
  * Hirsch et al. [2011] M. J. Hirsch, P. M. Pardalos, and M. G . C. Resende.  Correspondence of projected 3-d points and lines using a continuous grasp.  _International Transactions in Operational Research_ , 18(4):493–511, 2011.
  * Karmarkar et al. [1991] N. Karmarkar, M.G.C. Resende, and K.G. Ramakrishnan.  An interior point algorithm to solve computationally difficult set covering problems.  _Mathematical Programming_ , 52:597–618, 1991.
  * Kayışoğlu and Akgün [2021] B. Kayışoğlu and İ. Akgün.  Multiple allocation tree of hubs location problem for non-complete networks.  _Computers & Operations Research_, 136:105478, 2021.
  * Kolda et al. [2003] T. G. Kolda, R. M. Lewis, and V. Torczon.  Optimization by direct search: New perspectives on some classical and modern methods.  _SIAM review_ , 45(3):385–482, 2003.
  * Laporte et al. [2004] G. Laporte, J. J. Salazar-Gonzalez, and F. Semet.  Exact algorithms for the job sequencing and tool switching problem.  _IIE Transactions_ , 36(1):37–45, 2004.
  * Mangussi et al. [2023] A. D. Mangussi, H. Pola, H. G. Macedo, L. A. Julião, M. P. T. Proença, P. R. L. Gianfelice, B. V. Salezze, and A. A. Chaves.  Meta-heurísticas via chaves aleatórias aplicadas ao problema de localização de hubs em árvore.  In _Anais do Simpósio Brasileiro de Pesquisa Operacional_ , São José dos Campos, 2023. Galoá.
  * Mannino and Sassano [1995] C. Mannino and A. Sassano.  Solving hard set covering problems.  _Operations Research Letters_ , 18(1):1–5, 1995\.
  * Mecler et al. [2021] J. Mecler, A. Subramanian, and T. Vidal.  A simple and effective hybrid genetic search for the job sequencing and tool switching problem.  _Computers & Operations Research_, 127:105153, 2021.
  * Mehrotra and Trick [1998] A. Mehrotra and M. A. Trick.  Cliques and clustering: A combinatorial approach.  _Operations Research Letters_ , 22(1):1–12, 1998\.
  * Mladenović and Hansen [1997] N. Mladenović and P. Hansen.  Variable neighborhood search.  _Computers & Operations Research_, 24(11):1097–1100, 1997.
  * Nelder and Mead [1965] J. A. Nelder and R. Mead.  A simplex method for function minimization.  _The Computer Journal_ , 7(4):308–313, 1965.
  * Niven et al. [1991] I. Niven, H. S. Zuckerman, and H. L. Montgomery.  _An introduction to the theory of numbers_.  John Wiley & Sons, 1991.
  * Odijk and van Maaren [1998] M. A. Odijk and H. van Maaren.  Improved solutions to the steiner triple covering problem.  _Inf. Process. Lett._ , 65(2):67–69, 1998.
  * Oliveira et al. [2022] B. B. Oliveira, M. A. Carravilla, J. F. Oliveira, and M. G. C. Resende.  A C++ application programming interface for co-evolutionary biased random-key genetic algorithms for solution and scenario generation.  _Optimization Methods and Software_ , 37(3):1065–1086, 2022.
  * Ostrowski et al. [2009] J. Ostrowski, J. Linderoth, F. Rossi, and S. Smriglio.  Solving large steiner triple covering problems.  _Technical Report 1663, Computer Sciences Department, University of Wisconsin, Madison_ , 1663, 2009.
  * Ostrowski et al. [2010] J. Ostrowski, J. Linderoth, F. Rossi, and S. Smriglio.  Solving steiner triple covering problems.  _Optima_ , 1663, 2010.
  * Paiva and Carvalho [2017] G. S. Paiva and M. A. M. Carvalho.  Improved heuristic algorithms for the job sequencing and tool switching problem.  _Computers & Operations Research_, 88:208–219, 2017.
  * Penna et al. [2013] P. H. V. Penna, A. Subramanian, and L. S. Ochi.  An iterated local search heuristic for the heterogeneous fleet vehicle routing problem.  _Journal of Heuristics_ , 19:201–232, 2013.
  * Pessoa et al. [2017] L. S. Pessoa, A. C. Santos, and M. G. C. Resende.  A biased random-key genetic algorithm for the tree of hubs location problem.  _Optimization Letters_ , 11:1371–1384, 2017.
  * Privault and Finke [2000] C. Privault and G. Finke.  $k$-server problems with bulk requests: An application to tool switching in manufacturing.  _Annals of Operations Research_ , 96(1):255–269, 2000.
  * Reinelt [2003] G. Reinelt.  _The traveling salesman: Computational solutions for TSP applications_ , volume 840.  Springer, 2003.
  * Resende and Ribeiro [2016] M. G. C. Resende and C. C. Ribeiro.  _Optimization by GRASP: Greedy randomized adaptive search procedures_.  Springer Nature, 2016.
  * Schuetz et al. [2022] M. J. A. Schuetz, J. K. Brubaker, H. Montagu, Y. van Dijk, J. Klepsch, P. Ross, A. Luckow, M. G. C. Resende, and H. G. Katzgraber.  Optimization of robot trajectory planning with nature-inspired and hybrid quantum algorithms.  _Physical Review Applied_ , 18(5), 2022.
  * Siegal [1956] S. Siegal.  _Nonparametric statistics for the behavioral sciences_.  McGraw-hill, 1956.
  * Silva et al. [2013] R. M. A. Silva, M. G. C. Resende, P. M. Pardalos, and M. J. Hirsch.  A Python/C library for bound-constrained global optimization with continuous GRASP.  _Optimization Letters_ , 7:967–984, 2013.
  * Spears and De Jong [1991] W. M. Spears and K. A. De Jong.  On the virtues of parameterised uniform crossover.  In _Proceedings of the 4th International Conference on Genetic Algorithms_ , pages 230–236, San Diego, California, 1991.
  * Tang and Denardo [1988] C. S. Tang and E. V. Denardo.  Models arising from a flexible manufacturing machine, Part I: Minimization of the number of tool switches.  _Operations Research_ , 36(5):767–777, 1988.
  * Toso and Resende [2015] R. F. Toso and M. G. C. Resende.  A C++ application programming interface for biased random-key genetic algorithms.  _Optimization Methods and Software_ , 30(1):81–93, 2015.
  * Tsallis and Stariolo [1996] C. Tsallis and D. A. Stariolo.  Generalized simulated annealing.  _Physica A: Statistical Mechanics and its Applications_ , 233(1):395–406, 1996.
  * Xiang et al. [1997] Y. Xiang, D. Y. Sun, W. Fan, and X. G. Gong.  Generalized simulated annealing algorithm and its application to the Thomson model.  _Physics Letters A_ , 233(3):216–220, 1997.
  * Zhou et al. [2005] B. Zhou, L. Xi, and Y. Cao.  A beam-search-based algorithm for the tool switching problem on a flexible machine.  _The International Journal of Advanced Manufacturing Technology_ , 25:876–882, 2005.
  * Östergård and Vaskelainen [2011] P.R.J. Östergård and V.P. Vaskelainen.  Russian doll search for the steiner triple covering problem.  _Optimization Letter_ , 5:631–638, 2011.