††footnotetext: Please cite as: M. López-Ibáñez, F. Chicano, R. Gil-Merino.
The Asteroid Routing Problem: A Benchmark for Expensive Black-Box Permutation
Optimization. In J. L. Jiménez Laredo et al., editors, _Applications of
Evolutionary Computation_ , volume 13224 of Lecture Notes in Computer Science.
Springer Nature, Switzerland, 2022.
https://doi.org/10.1007/978-3-031-02462-7˙911institutetext: ITIS Software,
Universidad de Málaga, Spain
11email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# The Asteroid Routing Problem: A Benchmark for Expensive Black-Box
Permutation Optimization

Manuel López-Ibáñez 0000-0001-9974-1295 Francisco Chicano 0000-0003-1259-2990
Rodrigo Gil-Merino

###### Abstract

Inspired by the recent 11th Global Trajectory Optimisation Competition, this
paper presents the asteroid routing problem (ARP) as a realistic benchmark of
algorithms for expensive bound-constrained black-box optimization in
permutation space. Given a set of asteroids’ orbits and a departure epoch, the
goal of the ARP is to find the optimal sequence for visiting the asteroids,
starting from Earth’s orbit, in order to minimize both the cost, measured as
the sum of the magnitude of velocity changes required to complete the trip,
and the time, measured as the time elapsed from the departure epoch until
visiting the last asteroid. We provide open-source code for generating
instances of arbitrary sizes and evaluating solutions to the problem. As a
preliminary analysis, we compare the results of two methods for expensive
black-box optimization in permutation spaces, namely, Combinatorial Efficient
Global Optimization (CEGO), a Bayesian optimizer based on Gaussian processes,
and Unbalanced Mallows Model (UMM), an estimation-of-distribution algorithm
based on probabilistic Mallows models. We investigate the best permutation
representation for each algorithm, either rank-based or order-based. Moreover,
we analyze the effect of providing a good initial solution, generated by a
greedy nearest neighbor heuristic, on the performance of the algorithms. The
results suggest directions for improvements in the algorithms being compared.

###### Keywords:

Spacecraft Trajectory OptimizationUnbalanced Mallows Model Combinatorial
Efficient Global Optimization Estimation of Distribution Algorithms Bayesian
Optimization

##  1 Introduction

Several space programs from different countries focus on small bodies orbiting
around the Sun and the Earth. Apart from the search for answers about the
origin and evolution of the Solar System, there are other practical reasons
that space agencies consider. One of these reasons is related to the present
decrease of mineral resources on Earth. These resources, heavily used in
technological devices like mobile phones and computers, include silicon,
quartz, boronite and others, and are hard to mine for: the density of these
minerals is usually high, and they fell to inner strata during Earth
formation, which translates into an extra difficulty when mining for these
materials. As an alternative to Earth mining for techno-minerals, space
agencies have proposed asteroids and near Earth objects (NEOs) [3, 4] as near-
term mining targets. These objects contain varying amounts of rare earth
metals, minerals and water at surface levels, consequently much easier to mine
for than on Earth. A clear drawback to these plans is the transport from and
to the asteroids, which will imply high consumption of energy. A solution to
reduce the energy consumption is to build human settlements on large asteroids
and/or on orbit, together with space solar power stations that could fuel
spacecrafts and those settlements. This idea firstly appeared on a science-
fiction novel and later it was popularized as a “gedankenexperiment” by Dyson
in the 60’s. Nowadays, the idea of a swarm of space solar power stations is
called “the Dyson sphere”. In its last edition, the _11 ${}^{\text{th}}$
Global Trajectory Optimisation Competition_ (GTOC11,
https://gtoc11.nudt.edu.cn) proposed to computationally implement a Dyson
ring, which is a set of space stations (12 in this case) orbiting the Sun in a
circular orbit, to mine close to Earth passing-by asteroids. The Dyson ring
construction was proposed as an optimisation problem in the context of orbital
mechanics. The first part of the optimization problem consists in finding a
plan (trajectory) to visit a selected subset of asteroids with a spacecraft
launched from Earth. This spacecraft does not return to Earth. The second part
of the problem consists in activating a device, called _asteroid transfer
device_ or ATD, located in the visited asteroids to guide them to the space
stations in the Dyson ring. Once the Dyson ring is formed, regular visits to
the space stations in the Dyson ring are supposed to bring the minerals in the
asteroids to Earth. We are interested in this paper in solving the first
optimization problem, for which we provide a formulation as a challenging
expensive black-box permutation problem that we call _Asteroid Routing
Problem_ (ARP) that may be used to compare and benchmark optimization
approaches.

In recent years, there has been growing interest in tackling expensive black-
box permutation problems [11, 9, 1, 6]. In permutation problems, the decision
space is the space of permutations of a given length $n$, usually denoted with
$S_{n}$. Black-box optimization, as considered in this paper, tackles an
optimization problem without requiring an explicit analytical model of its
objective function and constraints. Finally, in expensive optimization, each
function evaluation is often the result of a costly simulation or physical
experiment, thus the total budget of evaluations allowed before returning the
best solution found is severely restricted, typically no more than 1000 [8,
9]. The combination of these three characteristics gives rise to a challenging
family of optimization problems. Typical approaches include estimation-of-
distribution algorithms [6] and combinatorial Bayesian optimizers [11] or
combinations of both [9].

A potential criticism of recent works is that the benchmark problems
considered —e.g., the travelling salesman problem, the linear ordering
problem, etc.— are taken from the classical combinatorial optimization
literature and they are neither truly expensive nor black-box. The ARP
proposed here is much more expensive to evaluate than classical combinatorial
optimization problems, even for short permutations. Moreover, there is no
explicit analytical model of the problem from the routing perspective, thus
classical techniques such as fast neighborhood evaluation or gradient
calculation are not possible. In our formal definition of the ARP, we aim to
preserve the most challenging features of the problem, without requiring
expert knowledge of the astrophysical details. Benchmark generation and
evaluation code is made publicly available. We also provide preliminary
results using state-of-the-art optimization methods. Our hope is that this
work will raise interest in the field of expensive black-box combinatorial
optimisation and its applications to trajectory optimisation.

This paper is structured as follows. Section 2 gives a brief introduction to
basic astrodynamical concepts required to understand the problem as an
optimization benchmark. We do not explain in detail the underlying
calculations because (1) these are standard formulas available in many physics
textbooks, (2) they are implemented by multiple open-source software libraries
and (3) they are not directly available to any of the optimization algorithms
considered here. The optimization model of the proposed Asteroid Routing
Problem is given in Section 3. The optimization algorithms evaluated in this
work are described in Section 4. We provide an experimental analysis of these
algorithms on the ARP in Section 5 and conclude in Section 6 with a summary of
findings and suggestions for further research motivated by these findings as
well as possible extensions of the ARP benchmark.

##  2 Background

In this section we will briefly revise the basic theory behind the
astrodynamics of the problem. We will start describing the two-body problem
and its solution leading to the Keplerian orbits. Then, we will describe one
type of maneuver used by the spacecrafts to change their trajectory in space
and we will end this section presenting the Lambert problem to determine the
impulses required to reach a point in space and time from any other point.

###  2.1 Two-body Problem

The dynamics of two bodies in space subject to the gravitational force is
determined by Newton’s second law of dynamics and gravitation law. Combined,
they lead to the following vector equation:

$\frac{d^{2}\mathbf{r}}{dt^{2}}=-\frac{\mu}{r^{3}}\mathbf{r},$ (1)

where $\mathbf{r}$ is the position vector of the second object respect to the
first one, $r=|\mathbf{r}|$ is the module of the position vector, $t$ is time,
and $\mu$ is usually called _gravitational parameter_ , which for our solar
system is typically assumed to be a constant due to the mass of the Sun being
much larger than the mass of any other object.

Despite the dynamics of the two bodies happening in 3D space, the conservation
of the angular momentum forces the movement to happen in a plane and the
analytical solutions of the problem are conic curves (also called Keplerian
orbits when they refer to the two-body problem): ellipse (and circumference as
a particular case), hyperbola and parabola111The trajectory of the second
object could also be linear, but this is not of practical interest, because
this would end with the object destroyed in the Sun.. Both, the hyperbola and
parabola are open curves, which means that the object following one of them
will come from outside the solar system and after approaching the Sun will
disappear outside the solar system without coming back again. Solar system
objects, like planets and asteroids, follow an elliptic (or circular) orbit
around the Sun.

Figure 1: The five parameters to describe an elliptic orbit.

When we fix a reference frame in space, elliptic orbits can be characterized
by five _orbital elements_ : ${a}$ (semi-major axis), $e$ (eccentricity), $i$
(inclination), $\omega$ (argument of periapsis222The periapsis is the nearest
point of an elliptic orbit to the object in the focus. When the object is the
Sun, the periapsis is also called perihelion.) and $\Omega$ (longitude of the
ascending node) (see Figure 1). The position of the object in the orbit is
characterized by a sixth parameter which changes with time: the true anomaly,
$\nu$. The term _epoch_ is used in this context to refer to a specific _date_
and _time_. The movement of an object, i.e., its position and velocity
vectors, in an elliptical orbit is completely determined by the six orbital
elements described above in a given epoch.

When there are more than two bodies, as it happens in our solar system, the
movement equations cannot be solved analytically, except for a few particular
cases. In that non-analytical case, the solution of the two-body problem is
still useful. The presence of additional bodies can be considered in most of
the cases a minor perturbation in the Hamiltonian describing the two-body
problem and, as a consequence, the orbital elements are slowly changing with
time [5, Chapter 11]. In our formulation of the problem, however, we do not
consider these perturbations, and we assume that all objects follow a
Keplerian orbit.

###  2.2 Maneuvers in Space

Artificial spacecrafts also follow a Keplerian orbit around the Sun when they
are far from the gravitational influence of the Earth or any other solar
system object. They usually have an engine that allows them to change the
orbit using operations called _maneuvers_. There are two main kinds of
maneuvers: _continuous_ and _impulsive_. In a continuous maneuver the engine
of the spacecraft is on during a long period of time (e.g., days) and
generates a force whose direction and magnitude can be controlled. We will
only consider impulsive maneuvers here.

Impulsive maneuvers cover a wide range of different space maneuvers, from
small orbital adjustments when two spacecrafts encounter each other, to large
orbital corrections between a low-Earth orbit and an interplanetary transit.
To carry out an impulsive maneuver, usually small engine firings take place,
changing only the direction and magnitude of the velocity vector
instantaneously. This means that the position vector does not change during
the impulse, which, although an idealisation, it works in most situations. The
advantage of treating the impulsive maneuver in this way is that it avoids
solving the equations of motion including the thrusts. Under these
assumptions, the trajectory of the spacecraft is still a Keplerian orbit where
the orbital elements change instantaneously at the point where the impulsive
maneuver happens.

###  2.3 Lambert Problem

At the heart of Astrodynamics, a key problem is to determine the Keplerian
orbit connecting two positions in a given period of time $\Delta t$. This is
the Lambert’s problem, also known as the Gauss’ problem [7].

We can use the solution of Lambert’s problem to determine the impulsive
maneuvers required for a spacecraft to do a rendezvous with an asteroid. The
transit orbit and impulsive maneuvers can be completely characterized by two
times: the epoch at which the transit starts $\tau$ and the transit time $t$.
These two times determine the two points $P_{1}$ and $P_{2}$ in space we need
to join. $P_{1}$ is the point at which the spacecraft is at $\tau$. $P_{2}$ is
the point at which the asteroid is at $\tau+t$. Lambert’s problem is solved to
go from $P_{1}$ to $P_{2}$ in time $t$. Two impulsive maneuvers are used to
move from $P_{1}$ to $P_{2}$. A first impulsive maneuver
$\Delta\mathbf{v}_{1}$ at $P_{1}$ inserts the spacecraft in the transit orbit
determined by solving Lambert’s problem. A second impulsive maneuver
$\Delta\mathbf{v}_{2}$ at $P_{2}$ inserts the spacecraft in the same orbit as
the asteroid (rendezvous). Let’s denote with $s$ the orbit of the spacecraft
(six orbital parameters in an epoch) and with $a$ the orbit of the
asteroid333We previously used $a$ to mean the semi-major axis of an orbit,
from now on we will use it to denote the complete orbit of an asteroid.. Then,
$\Delta\mathbf{v}_{1}$ and $\Delta\mathbf{v}_{2}$ are determined by $s$, $a$,
$\tau$ and $t$:

$(\Delta\mathbf{v}_{1},\Delta\mathbf{v}_{2})=\text{Lambert}(s,a,\tau,t)$ (2)

##  3 Asteroid Routing Problem

The _Asteroid Routing Problem_ consists in finding a route followed by a
spacecraft launched from Earth to visit a given set of $n$ asteroids
$A=\\{a_{1},a_{2},\ldots,a_{n}\\}$ that minimizes both the sum of velocity
impulses required by the route (related to consumption) and the total time
required to visit all of them. A solution to this problem is a pair
$(\pi,\mathbf{t})$, where $\pi\in S_{n}$ is a permutation representing the
order in which the asteroids are visited and
$\mathbf{t}\in\mathbb{R}^{2n}_{\geq 0}$ is a vector of $2n$ non-negative real
numbers representing parking and transit times to reach each asteroid.

We will consider a start epoch $\tau_{0}$, determined by the instance, when
the spacecraft is on Earth. The spacecraft remains on Earth during a time
$t_{1}$ after $\tau_{0}$. Then, it is launched at epoch $\tau_{0}+t_{1}$ to
reach the first asteroid in the sequence, $a_{\pi(1)}$, using a transit orbit
computed by solving Lambert’s problem between the Earth and the asteroid. In
the launch, the spacecraft experiments an impulse $\Delta\mathbf{v}_{1}$ that
inserts the spacecraft in a transit orbit to reach asteroid $a_{\pi(1)}$ in
time $t_{2}$ after the launch epoch. Once the asteroid is reached, a new
impulse $\Delta\mathbf{v}_{2}$ is used to make the spacecraft to follow the
orbit of asteroid $a_{\pi(1)}$.

In order to visit the $i$-th asteroid, $a_{\pi(i)}$, for $i>1$, the spacecraft
remains during time $t_{2i-1}$ in the orbit of asteroid $a_{\pi(i-1)}$ and
receives an impulse $\Delta\mathbf{v}_{2i-1}$ to reach asteroid $a_{\pi(i)}$
in time $t_{2i}$. Then, a new impulse $\Delta\mathbf{v}_{2i}$ changes the
orbit of the spacecraft to coincide with the one of $a_{\pi(i)}$.

For the transit between asteroids $a_{\pi(i-1)}$ and $a_{\pi(i)}$, the
impulses $\Delta\mathbf{v}_{2i-1}$ and $\Delta\mathbf{v}_{2i}$ are computed
using Equation (2). All times at odd positions of the time vector, $t_{2i-1}$,
are parking times, because the spacecraft is waiting in an asteroid (or
Earth’s) orbit, while the times at even positions of the time vector, $t_{2i}$
are transit times, because the spacecraft is traveling through a transit orbit
to reach the next asteroid.

The two objectives of the problem to minimize are the sum of all the velocity
impulses, which is related to energy consumption, and the sum of all parking
and transit times, which gives the elapsed time from launch to arriving at the
last asteroid. More formally, the two objectives are:

$\Delta
v=\sum_{i=1}^{2n}\left\lvert\Delta\mathbf{v}_{i}\right\rvert\quad\text{and}\quad
T=\sum_{i=1}^{2n}t_{i}\enspace,$ (3)

where $\Delta\mathbf{v}_{i}$ is computed as follows:

$(\Delta\mathbf{v}_{2i-1},\Delta\mathbf{v}_{2i})=\text{Lambert}\left(a_{\pi(i-1)},a_{\pi(i)},\tau_{i-1},t_{2i}\right)\enspace,$
(4)

where $\tau_{i-1}=\tau_{0}+\sum_{j=1}^{2i-1}t_{j}$ is the launch epoch from
orbit $a_{\pi(i-1)}$ and we consider that $a_{\pi(0)}$ is Earth’s orbit.

An instance of the problem is completely determined by the orbital parameters
at a given epoch of Earth and all asteroids, the starting epoch $\tau_{0}$ and
the gravitational parameter $\mu$.

In this paper, we will focus on a scalarized variant of the above bi-objective
problem that simply aggregates the two objectives in Eq. (3) as follows:

$f(\pi,\mathbf{t})=\Delta
v+\frac{2~{}\text{km/s}~{}}{30~{}\text{days}~{}}\cdot T,$ (5)

where the constant in front of $T$ is based on preliminary experiments.

##  4 Optimization Algorithms

In order to use the ARP as a benchmark problem for expensive black-box
combinatorial optimization, we tackle the problem in a hierarchical manner.
The _inner_ optimization problem decides the values of the vector of times
given a permutation of the asteroids, while the _outer_ problem aims to find
an optimal permutation of the asteroids. Since our focus here is the outer
problem, we always tackle the inner problem using the deterministic SLSQP
(Sequential Least Squares Programming) algorithm, thus ensuring that we obtain
the same objective function value given the same permutation. For the outer
problem, we evaluate and compare two algorithms for expensive black-box
optimization problems, namely, Unbalanced Mallows Model (UMM) [6] and
Combinatorial Efficient Global Optimization (CEGO) [11]. In addition, we study
the effect that a good initial solution has on the behavior of the black-box
algorithms. This initial solution is found by means of a gray-box greedy
nearest neighbor heuristic. All these algorithms are described below.

###  4.1 Sequential Least Squares Programming (SLSQP)

The inner optimization of the vector of times is done using the implementation
of SLSQP provided by SciPy [10] without an explicit gradient, which is
estimated by the SLSQP algorithm using 2-point finite differences with an
absolute step size of 1.49e-08. Each transfer is optimized independently, that
is, SLSQP solves a sequence of $n$ problems with two numerical decision
variables $(x_{0},x_{1})$, which correspond respectively to the parking and
transit times to visit a particular asteroid. Parking times are bounded by
$[0,730]$, while transit times are bounded by $[1,730]$, both measured in
days.

###  4.2 Greedy Nearest Neighbor Heuristic

A reasonably good permutation of the asteroids can be generated by iteratively
visiting the nearest asteroid, in Euclidean distance, to the last-visited
asteroid, after calculating the positions of all unvisited asteroids at the
arrival time at the last-visited asteroid. This heuristic is presented in
Algorithm 1. The algorithm keeps track of the orbit of the spacecraft $s$ and
its arrival time (epoch) $\tau$ at the last-visited asteroid, which are
initialized to the Earth’s orbit $a_{\pi(0)}$ and the initial epoch
$\tau_{0}$, and a set of unvisited asteroids $U$. At each iteration $i$, an
asteroid $a_{\pi(i)}$ is chosen from $U$ according to the minimum value of
$d(s,a,\tau)$, which calculates the Euclidean distance between the Cartesian
positions of the orbits of the spacecraft $s$ and asteroid $a$ when both are
considered at time $\tau$. After the asteroid $a_{\pi(i)}$ is chosen, we solve
the inner problem—using SLSQP as described above—to find the parking time
$t_{2i-1}$ at asteroid $a_{\pi(i-1)}$ and transit time $t_{2i}$ to asteroid
$a_{\pi(i)}$. The algorithm finishes when all asteroids are visited, which
results in a permutation $\pi$ of the asteroids and the corresponding vector
$\mathbf{t}$ of parking and transit times.

Algorithm 1 Greedy Nearest Neighbor Heuristic.

1: $s:=a_{\pi(0)}$ // Earth’s orbit

2: $\tau:=\tau_{0}$ // Epoch of the spacecraft

3: $U:=\\{1,\dotsc,n\\}$

4: for $i:=1$ to $n-1$ do

5: $\pi(i):=\operatorname*{arg\,min}_{j\in U}d(s,a_{j},\tau)$ // Euclidean
distance between orbits at epoch $\tau$

6: $(t_{2i-1},t_{2i}):=\textsf{SLSQP}(s$, $a_{\pi(i)})$ // Solve inner problem

7: $\tau:=\tau+t_{2i-1}+t_{2i}$

8: $U:=U\setminus\\{\pi(i)\\}$

9: $s:=a_{\pi(i)}$

10: return $(\pi,\mathbf{t})$

###  4.3 Unbalanced Mallows Model (UMM)

The UMM algorithm [6] is an estimation-of-distribution algorithm based on the
Mallows model and unbalanced Borda learning. The Mallows model for permutation
spaces is defined by a reference permutation $\sigma_{0}$ and a dispersion
parameter $\theta$, which are analogous to the mean and variance of the
Gaussian distribution. At each step of UMM, the permutations evaluated so far
together with their objective values are used to calculate a _weighted mean_
permutation $\hat{\sigma}_{0}$. From this $\hat{\sigma}_{0}$ and a value of
the dispersion parameter $\theta$, a new permutation is sampled at each
iteration from the Mallows model. The dispersion parameter $\theta$ is set by
UMM in such a way that, at the first iteration, the expected Kendall’s-$\tau$
distance of newly sampled permutations from the reference permutation is at
half of the expected distance of an uniform sample. In subsequent iterations,
the parameter $\theta$ is adjusted such that the expected Kendall’s-$\tau$
distance linearly decreases until the expected distance is 1 at the last
iteration.

Empirical results of UMM on expensive black-box variants of the linear
ordering problem and the permutation flowshop scheduling problem showed
results [6] comparable to the Combinatorial Efficient Global Optimization
(CEGO), with the additional benefit of being computationally much less
expensive than CEGO.

###  4.4 Combinatorial Efficient Global Optimization (CEGO)

CEGO extends classical Bayesian optimizers based on Gaussian processes to
unconstrained black-box combinatorial optimization problems [11]. A Gaussian
process may be used as a surrogate of the continuous landscape of an expensive
optimization problem. In order to model the landscape of a permutation space,
CEGO uses a distance measure for permutation spaces, such as the
Kendall’s-$\tau$ distance, to _interpolate_ between solutions. CEGO uses a
genetic algorithm (GA) to explore the landscape of the surrogate model. The GA
optimizes the expected improvement criterion, which takes into account both
the expected mean and variance of a solution. The solution returned by the GA
is then evaluated on the true objective function and this information is used
to update the Gaussian process surrogate model.

##  5 Experimental Study

In this section we describe the experiments performed and the results
obtained. Our goal is to answer two research questions:

RQ1:

How do the algorithms perform on the problem in a black-box setting (without
any problem knowledge)?

RQ2:

How do they perform when some a priori information about the problem, in the
form of a good initial solution, is introduced?

In Subsection 5.1 the methodology used to answer both research questions is
presented. Then, the experimental results to answer RQ1 and RQ2 are presented
and described in Subsections 5.2 and 5.3, respectively.

###  5.1 Experimental Methodology

We prepared an instance generator that, given an instance size $n$ and seed
for the random number generator, randomly selects $n$ asteroids from the 83
453 asteroids provided by the GTOC11 competition to create an instance of the
ARP. In this manner, we generated two instances (seeds 42 and 73) of each size
$n=\\{10,15,20,25,30\\}$. In the remainder of the paper, instances are named
$n$_seed, that is, instance 10_73 is an instance of size 10 generated with
seed 73.

In our first set of experiments, designed to answer RQ1, both UMM and CEGO
start the search from a set of 10 initial solutions generated by a _max-min-
distance sequential design_ , that is, solutions are added iteratively to a
set of existing solutions by maximizing the minimum distance to solutions
already in the set. The first solution in the set is generated uniformly at
random. Generating the max-min-distance design only takes into account the
distance between permutations and, thus, it does not require any evaluation of
the objective function. These 10 solutions are then evaluated on the objective
function function (Eq. 5) and they become the initial population of either
CEGO or UMM. In order to answer RQ2, we initialize the max-min-distance design
with the heuristic solution returned by the greedy nearest neighbor heuristic
described above.

Both UMM and CEGO stop after evaluating 400 permutations using the objective
function (Eq. (5)). Each evaluation requires solving the inner problem by
optimizing the vector of times $\mathbf{t}$ using a sequence of $n$ runs of
SLQSP. Each run $i$ of SLSQP optimizes the pair of parking and transit times
$(t_{2i-1},t_{2i})$, starting from the initial solution $(0,30)$ and
performing a maximum of 1000 iterations.

For each algorithm, we consider two representations, namely, whether the
permutations are represented as the _order_ or as the _ranking_ of the
asteroids in the visiting sequence. The _order_ -based representation is the
one used in the description of the problem in Section 3, i.e., $\pi(i)=j$
denotes that asteroid $j$ is visited at step $i$ of the route. Its _ranking_
-based counterpart is given by its inverse $\pi^{-1}$, where
$\pi(i)=j\Leftrightarrow\pi^{-1}(j)=i$. This distinction is important when
calculating distances between permutations, learning a probability
distribution, as done by UMM, or a surrogate model, as done by CEGO [6].

The GA used by CEGO to explore the surrogate model uses a budget of $10^{4}$
evaluations of the surrogate model, population size of 20 individuals, cycle
crossover for permutations with crossover rate of 0.5, swap mutation (i.e.,
exchanging two randomly selected elements) with mutation rate of $1/n$, and
tournament selection of size 2 with probability of 0.9. The population used by
the GA is different from the solutions evaluated by CEGO because the GA
explores the surrogate model, which does not require evaluations of the
expensive problem.

We use here the original implementation of GECO v2.4.2
(https://cran.r-project.org/package=CEGO). UMM is implemented in Python based
on the original code [6]. We have implemented the proposed ARP benchmark in
Python using the software package poliastro (v0.16) [2] for astrophysical
calculations.

Each experiment was repeated 30 times with different random seeds. The
experiments were run in the Picasso supercomputing facility of the University
of Málaga with 126 SD530 servers with Intel Xeon Gold 6230R (26 cores each) at
2.10GHz, 200 GB of RAM and an InfiniBand HDR100 network.

###  5.2 Results of the Black-box Setting

Table 1: Objective value of the best solution found by each algorithm in a single run. Mean (and standard deviation) over 30 independent runs. RS denotes random search. | CEGO | UMM |  |
---|---|---|---|---
Instance | Order | Rank | Order | Rank | RS
10_42 | 379.3 | (17.5) | 390.1 | (16.4) | 413.6 | (15.2) | 388.7 | (20.8) | 419.3 | (14.1)
10_73 | 346.0 | (17.7) | 361.0 | (14.3) | 374.5 | (16.6) | 361.4 | (18.9) | 375.2 | (15.3)
15_42 | 575.1 | (29.2) | 610.1 | (27.4) | 626.9 | (26.1) | 594.6 | (37.7) | 636.2 | (18.5)
15_73 | 582.4 | (33.1) | 591.6 | (29.1) | 626.2 | (24.7) | 595.3 | (32.6) | 636.2 | (25.2)
20_42 | 806.2 | (45.2) | 833.2 | (38.7) | 868.9 | (26.9) | 841.6 | (43.7) | 877.8 | (26.4)
20_73 | 839.2 | (49.8) | 838.0 | (43.0) | 883.0 | (28.4) | 876.8 | (32.7) | 897.2 | (26.3)
25_42 | 1048.8 | (51.2) | 1074.5 | (56.5) | 1121.2 | (28.3) | 1076.9 | (48.5) | 1122.1 | (20.8)
25_73 | 1069.1 | (49.5) | 1096.8 | (50.1) | 1131.4 | (35.8) | 1115.5 | (45.5) | 1139.7 | (25.6)
30_42 | 1272.2 | (65.8) | 1334.4 | (66.6) | 1372.0 | (24.3) | 1347.3 | (35.6) | 1357.3 | (38.4)
30_73 | 1320.9 | (63.2) | 1325.1 | (66.4) | 1385.9 | (27.5) | 1360.6 | (51.1) | 1398.6 | (24.9)

In this section we compare UMM and CEGO (combined with SLSQP) in a black-box
setting to solve the 10 instances of the ARP. Table 1 shows the results of
CEGO and UMM using the two representations, rank-based or order-based. As a
baseline, we show the results obtained by a simple _random search_ that
iteratively evaluates 400 random permutations. We can also see in Figure 2 the
evolution of the search for two particular instances. Each plot shows as a
solid line the mean, over 30 runs, of the best objective value found up to a
given number of function evaluations. The 95% confidence interval around the
mean is shown as a shaded region. We applied the non-paired Wilcoxon sumrank
test to compare the different algorithm configurations in each instance. The
p-values are shown in Table 2. We assume that the differences are
statistically significant when the p-value is below $0.01$ (confidence level).

Results show that, whenever there are significant differences (see columns two
and three of Table 2), rank-based UMM outperforms order-based UMM in terms of
the objective function, whereas order-based CEGO outperforms rank-based CEGO.
This result matches what we expect given the assumptions made by each
algorithm. That is, UMM internally assumes that the permutation $\pi$ it
generates represents a ranking ($\pi(3)=6$ means that object 3 has position 6)
and, thus, it makes sense that it works better when the fitness function is
aligned with this assumption [6]. CEGO, on the other hand, assumes that the
permutation it generates represents an order ($\pi(3)=6$ means that the third
object is 6), and that explains the good results of CEGO when the objective
function interprets the permutation as an order.

If we compare both algorithms, we observe that CEGO is usually better than UMM
using both representations (with only two exceptions for the ranking
representation in instances 10_42 and 15_42). The best CEGO variant is clearly
better than the best UMM variant, although not always with statistical
significance (see column four in Table 2).

Table 2: p-values of a Wilcoxon sumrank test of all the hypotheses checked in during the experimental study for each instance of the problem. Columns two and three correspond to a test comparing the two representations (rank and order) for CEGO and UMM in a black-box setting. The fourth column compares the best black-box UMM (using ranking) with the best black-box CEGO (using order). The final two columns compare UMM-ranking and CEGO-order in a black-box versus informed settings (including the greedy initialization). | Rank vs. Order | CEGO-order | Greedy vs. BB
---|---|---|---
Instance | CEGO | UMM | vs. UMM-rank | UMM | CEGO
10_42 | 6.46e-03 | 2.00e-06 | 4.48e-02 | 2.65e-02 | 1.18e-01
10_73 | 2.37e-03 | 1.33e-02 | 7.07e-04 | 8.34e-07 | 1.65e-02
15_42 | 3.77e-05 | 6.12e-04 | 5.14e-02 | 2.95e-11 | 1.72e-12
15_73 | 3.35e-01 | 2.75e-04 | 1.50e-01 | 6.52e-09 | 3.69e-10
20_42 | 1.42e-02 | 8.01e-03 | 2.74e-03 | 2.41e-05 | 3.25e-07
20_73 | 7.75e-01 | 4.67e-01 | 2.22e-03 | 2.84e-11 | 2.31e-11
25_42 | 3.71e-02 | 7.74e-05 | 3.45e-02 | 3.02e-11 | 4.04e-11
25_73 | 4.46e-02 | 1.77e-01 | 5.43e-04 | 2.26e-11 | 2.84e-11
30_42 | 5.43e-04 | 4.11e-03 | 1.30e-05 | 2.95e-11 | 3.02e-11
30_73 | 6.12e-01 | 3.71e-02 | 1.83e-02 | 2.78e-11 | 3.00e-11

Given these results, we select rank-based UMM and order-based CEGO for the
next step in our analysis.

Figure 2: Convergence of the best-so-far objective value over the number of
evaluations. Each line is the mean value over 30 runs and the shaded area is
the 95% confidence interval around the mean. The instances shown are 15_42
(top) and 25_73 (bottom).

Figure 3: Convergence of the best-so-far objective value over the number of
evaluations. Each line is the mean value over 30 runs and the shaded area is
the 95% confidence interval around the mean. The instances shown are 10_73
(top) and 30_42 (bottom).

###  5.3 Results of the Informed Setting

In this section, we evaluate the effect of providing a good initial solution
to the UMM-rank and CEGO-order variants selected above. This initial solution
is generated using the greedy nearest neighbor heuristic as the initial point
of the max-min-distance sequential design, as explained above, that creates
the 10 initial solutions of both UMM and CEGO.

Table 3: Objective value of the best solution found by each algorithm in a single run. For CEGO and UMM, the values are the mean (and standard deviation) over 30 independent runs. |  |  | CEGO |  | UMM
---|---|---|---|---|---
Inst. | GreedyNN |  | Order | Ord+Greedy |  | Rank | Rank+Greedy
10_42 | 391.3 |  | 379.3 | (17.5) |  | 374.9 | (6.9) |  | 388.7 | (20.8) |  | 382.4 | (1.1)
10_73 | 398.3 |  | 346.0 | (17.7) |  | 355.9 | (18.3) |  | 361.4 | (18.9) |  | 385.0 | (6.4)
15_42 | 508.1 |  | 575.1 | (29.2) |  | 497.2 | (1.1) |  | 594.6 | (37.7) |  | 501.7 | (3.8)
15_73 | 576.4 |  | 582.4 | (33.1) |  | 525.6 | (4.7) |  | 595.3 | (32.6) |  | 545.0 | (9.3)
20_42 | 841.7 |  | 806.2 | (45.2) |  | 737.0 | (31.0) |  | 841.6 | (43.7) |  | 797.6 | (28.3)
20_73 | 691.5 |  | 839.2 | (49.8) |  | 661.8 | (7.8) |  | 876.8 | (32.7) |  | 684.7 | (6.1)
25_42 | 946.3 |  | 1048.8 | (51.2) |  | 881.5 | (23.3) |  | 1076.9 | (48.5) |  | 921.7 | (15.2)
25_73 | 918.3 |  | 1069.1 | (49.5) |  | 873.6 | (15.3) |  | 1115.5 | (45.5) |  | 911.9 | (10.0)
30_42 | 1131.7 |  | 1272.2 | (65.8) |  | 1084.6 | (13.4) |  | 1347.3 | (35.6) |  | 1115.4 | (11.0)
30_73 | 1024.7 |  | 1320.9 | (63.2) |  | 967.7 | (17.9) |  | 1360.6 | (51.1) |  | 1017.3 | (6.9)

Figure 3 shows the evolution of the algorithms. In addition to the results
provided by CEGO and UMM with and without the initial greedy solution, we plot
the objective function value of the greedy solution as a horizontal line.
Table 3 provides a summary of the results of the two best variants of CEGO and
UMM together with the result obtained by the Greedy Nearest Neighbor
Heuristic. A first observation is that UMM only reaches the fitness level of
the greedy solution after 10 evaluations. The reason is that we evaluate first
the 9 random solutions generated by the max-min-distance design to show how
much the greedy solution improves over them. Moreover, without knowledge of
the greedy solution, CEGO or UMM are only able to match the greedy solution
under 400 evaluations on a few small instances ($n\leq 20$). Unsurprisingly,
when starting from the greedy solution, both algorithms outperform their
_uninformed_ counterparts. However, UMM struggles to find any improving
solutions, whereas CEGO is typically able to further improve the greedy
solution, although not by much. In the case of UMM, this behavior is explained
by the fact that the sampling of new permutations is strongly biased by the
relative objective values of the permutations already evaluated. Since the
greedy solution is so much better than any other initial permutation, the
algorithm remains stuck at the initial permutation. In the case of CEGO,
starting from a good initial solution provides a good starting point around
which the surrogate model can be built. However, in several instances, after a
quick initial improvement of the greedy solution, CEGO appears stuck and
unable to further find any improvements. This behavior of CEGO was already
reported for other problems [6] and it is attributed to the inability of the
surrogate model to accurately estimate the underlying objective function
beyond a certain point, thus leading to a blind search.

Table 4: Mean CPU-time (minutes) of a single run of each algorithm. |  |  | CEGO |  | UMM
---|---|---|---|---|---
Instance | GreedyNN |  | Order | Ord+Greedy |  | Rank | Rank+Greedy
10_42 | 0.1 |  | 2365.4 | 2147.8 |  | 26.0 | 26.2
10_73 | 0.1 |  | 2372.9 | 2365.3 |  | 24.2 | 25.7
15_42 | 0.1 |  | 2445.9 | 2394.1 |  | 34.5 | 32.4
15_73 | 0.1 |  | 2375.6 | 2404.2 |  | 35.5 | 36.2
20_42 | 0.2 |  | 2411.4 | 2348.4 |  | 47.8 | 48.7
20_73 | 0.2 |  | 2223.9 | 2387.4 |  | 48.2 | 50.8
25_42 | 0.2 |  | 2338.1 | 2452.5 |  | 59.4 | 58.9
25_73 | 0.2 |  | 2397.3 | 2657.6 |  | 47.3 | 60.3
30_42 | 0.2 |  | 2375.9 | 2553.7 |  | 61.9 | 69.1
30_73 | 0.2 |  | 2138.3 | 2611.0 |  | 56.0 | 67.4

We also report the mean runtime of CEGO, UMM and the greedy heuristic in Table
4. The runtime of the greedy algorithm is relatively fast, taking around 6
seconds to generate and evaluate a single solution. UMM requires between 24
and 69 minutes per run on average, depending on the instance size, where each
run involves 400 solution evaluations. The runtime per evaluation is larger
than in GreedyNN due to the update of and sampling from the probabilistic
model. CEGO also evaluates 400 solutions per run, however, it is more than 100
times slower than UMM, requiring almost 2 CPU-days per run on average for the
largest instances tested here. Although part of this difference is due to
implementation choices, the main difference is the fact that building the
Gaussian Process model in CEGO is significantly more expensive than building
the probabilistic model in UMM.

Finally, we show in Fig. 4 a visualization of three solutions found by the
greedy heuristic, CEGO-order and UMM-rank, the two last ones starting from the
greedy solution. In each plot, the legend shows the order in which asteroids
are visited (from top to bottom), the epoch at which an impulse was applied to
change orbit, and the arrival epoch at each asteroid. Solid lines indicate
either parking or transfer orbits. When the impulse and arrival epochs
coincide, the parking time at the previous orbit was zero and no parking orbit
is shown. Earth’s orbit is shown as a dashed blue line and the Sun’s position
as a yellow point. Distances are shown at scale.

(a) Greedy Nearest Neighbor.

(b) CEGO-order $+$ Greedy.

(c) UMM-rank $+$ Greedy.

Figure 4: Example solutions for instance 10_73.

##  6 Conclusions

We have proposed in this paper a realistic benchmark for evaluating algorithms
designed to tackle expensive black-box optimization problem in permutation
space. The results highlight already some limitations of the evaluated
algorithms. First, although CEGO consistently obtains better results than UMM
as long as the right choice of representation is made, this comes at a cost of
100 times more computation time. Second, without any additional information,
very rarely either algorithm is able to match the quality of the solution
found by a greedy approach, suggesting that their results are far from
optimal. Finally, only CEGO makes some use of the information provided by a
good initial solution, whereas UMM is very rarely able to improve over it.

Future work should evaluate other black-box algorithms for permutation
problems in the ARP. Improvements to CEGO and UMM to make better use of good
initial solutions would also be welcome. In fact, it is possible to generate
multiple greedy solutions that may be used as starting points by the black-box
optimizers. Finally, we plan to change the optimizer of the inner problem from
the current deterministic SLSQP to a stochastic one, which changes the context
of the outer problem from a deterministic permutation problem to a stochastic
one. Under such context, a black-box optimizer for the outer problem must
decide whether it is worth evaluating again already seen permutation in the
hope that a new evaluation may further improve the objective function.

Reproducibility. Datasets and source code for reproducing the experiments
reported are available from 10.5281/zenodo.5725837 .

Acknowledgements. This work is partially funded by the Universidad de Málaga,
Consejería de Economía y Conocimiento de la Junta de Andalucía and FEDER
(grant UMA18-FEDERJA-003) and MCIN/AEI/10.13039/501100011033 (grant PID
2020-116727RB-I00). Thanks to the Supercomputing and Bioinnovation Center
(SCBI) of the Univerity of Málaga for their provision of computational
resources and support. M. López-Ibáñez is a “Beatriz Galindo” Senior
Distinguished Researcher (BEAGAL 18/00053) funded by the Spanish Ministry of
Science and Innovation (MICINN).

## References

  * Bartz-Beielstein and Zaefferer [2017] Bartz-Beielstein, T., Zaefferer, M.: Model-based methods for continuous and discrete global optimization. Applied Soft Computing 55, 154–167 (Jun 2017), 10.1016/j.asoc.2017.01.039
  * Cano Rodríguez et al. [2015] Cano Rodríguez, J.L., et al.: poliastro: Astrodynamics in Python. Zenodo (2015), 10.5281/zenodo.174
  * European Space Agency [2019] European Space Agency: Hera Mission. https://whttps://www.esa.int/Safety˙Security/Hera/Hera (2019), [Online; accessed 22-Nov-2021]
  * European Space Agency [2021] European Space Agency: Science & Exploration: Asteroids. https://www.esa.int/Science˙Exploration/Human˙and˙Robotic˙Exploration/Exploration/Asteroids (2021), [Online; accessed 22-Nov-2021]
  * Goldstein [1980] Goldstein, H.: Classical Mechanics. Addison-Wesley (1980)
  * Irurozki and López-Ibáñez [2021] Irurozki, E., López-Ibáñez, M.: Unbalanced mallows models for optimizing expensive black-box permutation problems. In: Proceedings of GECCO (2021), 10.1145/3449639.3459366
  * Izzo [2014] Izzo, D.: Revisiting Lambert’s problem. arXiv 1403.2705 [astro-ph.EP] (2014)
  * Knowles et al. [2009] Knowles, J.D., Corne, D., Reynolds, A.P.: Noisy multiobjective optimization on a budget of 250 evaluations. In: Proccedings of EMO (LNCS 5467), pp. 36–50 (2009)
  * Pérez Cáceres et al. [2015] Pérez Cáceres, L., López-Ibáñez, M., Stützle, T.: Ant colony optimization on a limited budget of evaluations. Swarm Intelligence 9(2-3), 103–124 (2015), 10.1007/s11721-015-0106-x
  * Virtanen et al. [2020] Virtanen, P., et al.: SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nature Methods 17, 261–272 (2020), 10.1038/s41592-019-0686-2
  * Zaefferer et al. [2014] Zaefferer, M., Stork, J., Friese, M., Fischbach, A., Naujoks, B., Bartz-Beielstein, T.: Efficient global optimization for combinatorial problems. In: Proceedings of GECCO, pp. 871–878 (2014), 10.1145/2576768.2598282


&=\frac{3}{\sqrt{\pi}}(I_{\texttt{+-}}(h)+I_{\texttt{-+}}(h))\\\
&=\log|\mathbf{\Gamma}_{\mathfrak{a},h}(\mathfrak{fb}^{-1})|^{2},\end{split}$

proving Theorem 3.1.

Remark. The proof of Lemma 4.2.5 implies the following property of the numbers
$\displaystyle\mathbf{\Gamma}_{\mathfrak{a},h}(L)$. Note that the function
$\displaystyle\Gamma_{a,b}$ defined in (2.4) has the property

$\begin{split}&\Gamma_{a,b}(-w,z;L)\Gamma_{a,b}(w,z;L)\\\
&=\frac{\prod_{\delta\in(C_{\texttt{+-}}(a,b)\cap
H(b))/\mathbf{Z}\gamma}\left(1-e^{-2\pi
i(\delta(z)-w)/\gamma(z)}\right)\left(1-e^{-2\pi
i(\delta(z)+w)/\gamma(z)}\right)}{\prod_{\delta\in(C_{\texttt{-+}}(a,b)\cap
H(a))/\mathbf{Z}\gamma}\left(1-e^{-2\pi
i(\delta(z)-w)/\gamma(z)}\right)\left(1-e^{-2\pi
i(\delta(z)+w)/\gamma(z)}\right)}.\end{split}$

Using the notation in the above proof, this implies that

(4.32)
$\mathbf{\Gamma}_{\mathfrak{a},h}(L)\mathbf{\Gamma}_{\mathfrak{a},h}(-h(\mathbf{x}),\mathbf{x};L)=\frac{\theta_{0}^{(N)}(z_{0},\sigma)}{\theta_{0}^{(N)}(z_{0},-\tau)}=\zeta,$

where $\displaystyle\zeta^{3}=1$; moreover, $\displaystyle\zeta=1$ if
$\displaystyle(6,\mathrm{N}(\mathfrak{a}))=1$.

####  4.2.6. A geometric interpretation of the elliptic gamma function

Felder–Henriques–Rossi–Zhu [24] gave a geometric interpretation of the
elliptic gamma function. A _triptic curve_ is a pointed holomorphic stack of
the the form $\displaystyle[\mathbf{C}/\mathbf{Z}^{3}]$ where the action is
induced by a homomorphism $\displaystyle\mathbf{Z}^{3}\to\mathbf{C}$ of real
rank $\displaystyle 2$. It plays a role analogous to that of the elliptic
curve for the theta function. The moduli space of oriented triptic curves is
the complex stack

(4.33)
$\left[\mathrm{SL}_{3}(\mathbf{Z})\backslash(\mathbf{C}\mathrm{P}^{2}-\mathbf{R}\mathrm{P}^{2})\right].$

The elliptic gamma function is a section of a hermitian abelian gerbe over the
universal oriented triptic curve. Our construction can be understood in terms
of a modified version of this gerbe, defined using the ideal
$\displaystyle\mathfrak{a}$. The stack (4.33) is however far from being
algebraic, thus proving our conjecture certainly requires new ideas. We
nevertheless want to point out that — up to replacing
$\displaystyle\mathrm{SL}_{3}(\mathbf{Z})$ by
$\displaystyle\mathrm{SL}_{2}(\mathbf{Z}[1/p])$ and
$\displaystyle\mathbf{C}\mathrm{P}^{2}-\mathbf{R}\mathrm{P}^{2}$ by the
product of $\displaystyle\mathscr{H}$ by the Drinfeld upper half-plane — the
above geometric picture seems quite parallel to that of the recent work of
Darmon and Vonk [15] on explicit class field theory for real quadratic fields.

##  5\. A few additional numerical examples

In this Section we present numerical evidence for Conjecture 2.3 in some
examples. The computation relies on the summation formula [25, eq. (15)] for
the logarithm of the elliptic gamma function, whose convergence is slow if
$\displaystyle\tau$ and $\displaystyle\sigma$ have small imaginary parts. This
is the main bottleneck towards computing with larger input data. It would be
interesting to find faster methods to evaluate
$\displaystyle\mathbf{\Gamma}_{\mathfrak{a},h}(\mathfrak{fb}^{-1})$, for
example by making more efficient use of the
$\displaystyle\mathrm{SL}_{3}(\mathbf{Z})$-modularity of
$\displaystyle\Gamma(z,\tau,\sigma)$ [25].

###  5.1. $\displaystyle K=\mathbf{Q}(\beta)$ where
$\displaystyle\beta^{3}-3=0$

The discriminant is $\displaystyle-3^{5}$, the class number is $\displaystyle
1$ and $\displaystyle\beta^{2}-2$ is a fundamental unit. Take
$\displaystyle\mathfrak{f}$ to be the unique ideal of norm $\displaystyle 3$.
Then the fundamental unit belongs to
$\displaystyle\mathscr{O}_{\mathfrak{f}}^{+,\times}$.

The narrow ray class field $\displaystyle K(\mathfrak{f})$ is a quadratic
extension of $\displaystyle K$. We begin with this simple example because we
carried out our first tests on it.

We first take the smoothing ideal $\displaystyle\mathfrak{a}$ to be the single
prime ideal with norm $\displaystyle 5$. Then the complex numbers
$\displaystyle\mathbf{\Gamma}_{\mathfrak{a},h}(\mathfrak{f})$ appear to be
independent of $\displaystyle h$ and equal to $\displaystyle
0.18108\ldots+i\cdot 0.21746\ldots$, a complex number that Pari-GP identifies
as a root of

$\displaystyle x^{6}-3x^{5}+6x^{4}+17x^{3}+6x^{2}-3x+1$

to 1000 digits. Here again the splitting field is the quadratic extension
$\displaystyle K(\mathfrak{f})$ of $\displaystyle K$. Since the class number
formula implies Stark’s conjecture in this setting [37, Thm. 5.4], our Theorem
3.1 shows that $\displaystyle|\mathbf{\Gamma}_{\mathfrak{a},h}(\mathfrak{f})|$
is algebraic.

Let us now take the smoothing ideal $\displaystyle\mathfrak{a}$ to be the
single prime ideal with norm $\displaystyle 2$. The algdep command of Pari-GP
allows to show that the complex numbers
$\displaystyle\mathbf{\Gamma}_{\mathfrak{a},h}(\mathfrak{f})$ are very close
to a root of a degree $\displaystyle 24$ irreducible polynomial over
$\displaystyle\mathbf{Q}$. As predicted by one of the remarks following our
conjecture, up to multiplication by $\displaystyle\pm 1$, the complex number
$\displaystyle\mathbf{\Gamma}_{\mathfrak{a},h}(\mathfrak{f})^{4}$ is
independent of $\displaystyle h$, and Pari-GP identifies it as
$\displaystyle\pm\mathbf{u}_{2}$ where $\displaystyle\mathbf{u}_{2}$ is a root
of

$\displaystyle x^{6}+24x^{5}+168x^{4}+98x^{3}+168x^{2}+24x+1$

to 1000 digits. The splitting field of this polynomial is the quadratic
extension $\displaystyle K(\mathfrak{f})$ of $\displaystyle K$.

###  5.2. $\displaystyle K=\mathbf{Q}(\beta)$ where
$\displaystyle\beta^{3}-7=0$

The discriminant is $\displaystyle-3^{3}\cdot 7^{2}$, the class number is
$\displaystyle 3$ and $\displaystyle 2-\beta$ is a fundamental unit. Take
$\displaystyle\mathfrak{f}$ to be the unique ideal of norm $\displaystyle 3$.
The corresponding narrow ray class field $\displaystyle K(\mathfrak{f})$ has
relative degree $\displaystyle 6$ over $\displaystyle K$ and the fundamental
unit is a generator of $\displaystyle\mathscr{O}_{\mathfrak{f}}^{+,\times}$.

In this case it turns out that one can find many choices of $\displaystyle h$
and $\displaystyle L$ such that that the set $\displaystyle
F/\mathbf{Z}\gamma$ appearing in (2.7) is a singleton. This explains why there
is only one occurence of the basic elliptic gamma function
$\displaystyle\Gamma$ in the example of the introduction. There the smoothing
ideal $\displaystyle\mathfrak{a}$ is the single ideal with norm $\displaystyle
5$, and the complex numbers
$\displaystyle\mathbf{\Gamma}_{\mathfrak{a},h}(\mathfrak{f})$ appear
numerically to be independent of $\displaystyle h$.

Let $\displaystyle\mathfrak{b}$ be the unique ideal of norm $\displaystyle 2$.
Its class is a generator of the ray class group. We identify the alleged units
$\displaystyle\mathbf{u}_{\mathfrak{f}\mathfrak{b}^{-k},\mathfrak{a}},\ (0\leq
k\leq 5)$, with their complex images given by the elliptic gamma values. As
predicted by our conjecture one verifies that these numbers coincide with the
roots of (1.3) to 1000 digits and that these roots satisfy parts (2) and (3)
of our conjecture. The number (1.2) is
$\displaystyle\mathbf{u}_{\mathfrak{fb}^{-5},\mathfrak{a}}$ (taking
$\displaystyle\beta\in\mathbf{C}$ to have negative imaginary part). It is very
close to the image, under $\displaystyle\sigma_{\mathfrak{b}}^{-1}$, of the
root of (1.3) that is

$\displaystyle\mathbf{u}_{\mathfrak{f},\mathfrak{a}}\approx-0.06261629382\ldots-i\cdot
0.2669850711\ldots.$

We have similarly tested our conjecture with other smoothing ideals, e.g.
$\displaystyle\mathfrak{a}_{11}$ the single prime ideal with norm
$\displaystyle 11$.

###  5.3. $\displaystyle K=\mathbf{Q}(\beta)$ where
$\displaystyle\beta^{3}-2=0$

The discriminant is $\displaystyle-108$, the class number is $\displaystyle 1$
and $\displaystyle\beta-1$ is a fundamental unit. Fix a complex embedding
$\displaystyle K\hookrightarrow\mathbf{C}$.

Following Ren and Sczech [31, §4.1] we first take
$\displaystyle\mathfrak{f}=(3)$. The corresponding narrow ray class field
$\displaystyle K(\mathfrak{f})$ has relative degree $\displaystyle 6$ over
$\displaystyle K$ and $\displaystyle(\beta-1)^{3}$ is a generator of
$\displaystyle\mathscr{O}_{\mathfrak{f}}^{+,\times}$. The field $\displaystyle
K(\mathfrak{f})$ contains $\displaystyle 18$ roots of unity. We fix an
embedding $\displaystyle K(\mathfrak{f})\hookrightarrow\mathbf{C}$ extending
the fixed complex embedding of $\displaystyle K$.

Let our smoothing ideal $\displaystyle\mathfrak{a}$ be the single prime ideal
with norm $\displaystyle 5$. The class of the ideal
$\displaystyle\mathfrak{b}=(\beta)$ generates the ray class group. Using our
conjecture we obtain complex numbers that coincide with the roots of

$\begin{split}x^{6}&-(3\beta^{2}-9\beta+6)x^{5}+(23\beta^{2}-17\beta-16)x^{4}\\\
&+(48\beta^{2}-43\beta-23)x^{3}+(23\beta^{2}-17\beta-16)x^{2}-(3\beta^{2}-9\beta+6)x+1\end{split}$

to $\displaystyle 1000$ decimal digits. This polynomial is irreducible over
$\displaystyle K$ with splitting field $\displaystyle K(\mathfrak{f}).$ Its
image under the real embedding of $\displaystyle K$ has all its roots lying on
the unit circle.

Ren and Sczech [31] have checked numerically that among the eighteen roots of

$x^{18}-7767x^{17}+51550065x^{16}-199524692622x^{15}+520755985257966x^{14}\\\
-1828056747902004x^{13}+24870880029533226x^{12}-80588629212013080x^{11}\\\
+116076408275027511x^{10}-118102314911180623x^{9}+116076408275027511x^{8}\\\
-80588629212013080x^{7}+24870880029533226x^{6}-1828056747902004x^{5}\\\
+520755985257966x^{4}-199524692622x^{3}+51550065x^{2}-7767x+1,$

six have absolute value $\displaystyle 1$ and the twelve others are the
complex images of the Galois transforms of “the” alleged Stark unit
$\displaystyle\mathbf{u}_{\rm Stark}$ (uniquely defined up to the
$\displaystyle 18$ roots of unity in $\displaystyle K(\mathfrak{f})$) that is
denoted by $\displaystyle\eta_{6}$ in [31]. We verify numerically that

$\displaystyle\mathbf{u}_{\rm
Stark}^{\mathrm{N}(\mathfrak{a})-\sigma_{\mathfrak{a}}}\approx\mathbf{\Gamma}_{\mathfrak{a},h}(\mathfrak{f})^{18}.$

We have also numerically checked other ramification ideals, e.g.
$\displaystyle\mathfrak{f}=(5)$.

###  5.4. $\displaystyle K=\mathbf{Q}(\beta)$ where
$\displaystyle\beta^{3}-\beta+1=0$

It is the complex cubic field with minimal discriminant, equal to
$\displaystyle-23$. It is of class number $\displaystyle 1$ and
$\displaystyle-\beta$ is a fundamental unit. Fix an embedding $\displaystyle
K\hookrightarrow\mathbf{C}$.

As in [31, §4.2] we consider $\displaystyle\mathfrak{f}=(5)$. The
corresponding narrow ray class field $\displaystyle K(\mathfrak{f})$ has
relative degree $\displaystyle 4$ over $\displaystyle K$ and
$\displaystyle\beta^{24}$ is a generator of
$\displaystyle\mathscr{O}_{\mathfrak{f}}^{+,\times}$. The field $\displaystyle
K(\mathfrak{f})$ contains $\displaystyle 10$ roots of unity. We fix an
embedding $\displaystyle K(\mathfrak{f})\hookrightarrow\mathbf{C}$ extending
the fixed complex embedding of $\displaystyle K$.

Let our smoothing ideal $\displaystyle\mathfrak{a}$ be the single prime ideal
with norm $\displaystyle 7$. Using our conjecture we obtain complex numbers
that coincide with the complex roots of

$\displaystyle
x^{4}+(70\beta^{2}+15\beta-104)x^{3}-(65\beta^{2}+155\beta+89)x^{2}+(70\beta^{2}+15\beta-104)x+1$

to $\displaystyle 1000$ decimal digits. This polynomial is irreducible over
$\displaystyle K$ with splitting field $\displaystyle K(\mathfrak{f}).$ Its
image under the real embedding of $\displaystyle K$ has all its roots lying on
the unit circle.

Ren and Sczech have checked numerically that, in this example, “the” alleged
Stark unit $\displaystyle\mathbf{u}_{\rm Stark}$ (uniquely defined up to the
$\displaystyle 10$ roots of unity in $\displaystyle K(\mathfrak{f})$) is a
root of

$\displaystyle
x^{4}-(395\beta^{2}-60\beta-776)x^{3}+(495\beta^{2}-385\beta-1374)x^{2}-(395\beta^{2}-60\beta-776)x+1.$

It is denoted by $\displaystyle\eta_{4}$ in [31]. We verify numerically that

$\displaystyle\mathbf{u}_{\rm
Stark}^{\mathrm{N}(\mathfrak{a})-\sigma_{\mathfrak{a}}}\approx\mathbf{\Gamma}_{\mathfrak{a},h}(\mathfrak{f})^{10}.$

## Appendix: Eisenstein’s Jugendtraum

This section provides some historical background underpinning our work, deeply
rooted in Gotthold Eisenstein’s largely ignored paper [22].

In 1844, three years before his masterpiece [23] laying alternative
foundations for a theory of elliptic functions, Eisenstein published in
Crelle’s journal a short prospective note [22] where he first sketches his own
construction of classical trigonometric and elliptic functions by considering
nested conditionally convergent infinite products.

In the last pages, he indicates how one could go beyond and investigate higher
dimensional cases. Aware of the fact that no meromorphic functions over
$\displaystyle\mathbf{C}$ can be triply periodic, he first states:

> Following the remark [of his §1] related to simple and double infinite
> products relative to circular and elliptic functions, we should consider as
> an analog of greater degree the quotients of quotients of infinite triple
> products555Eisenstein’s emphasis; this sentence has been ridiculed by Jacobi
> in [26], who misjudged Eisenstein’s original point of view; see also [1]. of
> the form
>
> (5.1)
> $\prod\left(1-\frac{x}{\lambda+\lambda^{\prime}A+\lambda^{\prime\prime}A^{\prime}}\right),$
>
> where $\displaystyle A$ and $\displaystyle A^{\prime}$ are constants and
> $\displaystyle\lambda,\lambda^{\prime},\lambda^{\prime\prime}$ represent
> indices with the ordinary meaning attached to them. […]
>
> Unlike the simple or double products, these multiple products would not have
> a determined analytical sense if one wishes to assign to the indices all
> possible values.”

Eisenstein goes further and suggests:

> Nevertheless nothing refrains from considering these multiple products as
> soon as the indices are constrained by certain restrictions that enable to
> discard precisely those circumstances contrary to convergence : for example,
> certain inequality conditions could be assigned to the indices
> $\displaystyle\lambda,\lambda^{\prime},\lambda^{\prime\prime}$ and adjoined
> to the product, to the effect that the multiplication sign is extended to
> the values of the indices satisfying these conditions, the other values
> being omitted. Generally speaking, one cannot say here anything further
> about the choice of these conditional equations. There exists however a
> whole class of similar functions which have very close connections to
> certain results in Number Theory.

According to this paragraph, it looks plausible to us that Eisenstein had
anticipated the construction of (quotients of) the elliptic gamma function by
means of the simple manipulations for multiple infinite products he pioneered.
As we shall now exemplify, Eisenstein’s own convergence method quickly allows
for an explicit relation between a triple product of the shape (5.1), limited
by the inequality conditions
$\displaystyle\\{(\lambda,\lambda^{\prime},\lambda^{\prime\prime})\in\mathbf{Z}^{3}\;\big{|}\;\lambda^{\prime},\lambda^{\prime\prime}>0\\},$
and the function $\displaystyle\Gamma(z,A,A^{\prime}).$

Recall (see for instance [38, 10] for a survey) that in dimension
$\displaystyle 1$ his method provides the classical identity for trigonometric
functions

$\begin{split}\prod_{\lambda\in\mathbf{Z}}^{(\textrm{eis})}\left(1-\frac{z}{\lambda+u}\right)&=\lim_{M\rightarrow+\infty}\prod_{\lambda=-M}^{M}\left(1-\frac{z}{\lambda+u}\right)\\\
&=\frac{\sin\pi(u-z)}{\sin(\pi u)}\\\ &=e^{i\pi z}\left(\frac{1-e^{2\pi
i(u-z)}}{1-e^{2\pi iu}}\right).\end{split}$

In order to improve the convergence of the next infinite products, it is
convenient to get rid of the exponential factor by introducing a smoothed
version with an extra integral parameter $\displaystyle N>1:$

$\prod_{\lambda\in\mathbf{Z}}^{(\textrm{eis})}\frac{\left(1-\frac{z}{\lambda+u}\right)^{N}}{\left(1-\frac{Nz}{\lambda+Nu}\right)}=\frac{(1-e^{2\pi
i(u-z)})^{N}(1-e^{2\pi iNu})}{(1-e^{2\pi iu})^{N}(1-e^{2\pi i(Nu-Nz)})}\cdot$

Now we substitute $\displaystyle
u=A\lambda^{\prime}+A^{\prime}\lambda^{\prime\prime}+v,$ with $\displaystyle
A,A^{\prime}\in\mathscr{H},$ and proceed to the absolutely convergent double
product over
$\displaystyle\lambda^{\prime},\lambda^{\prime\prime}\in\mathbf{Z}$
constrained by the inequality
$\displaystyle\lambda^{\prime},\lambda^{\prime\prime}>0.$ This immediately
leads to

$\displaystyle\displaystyle Q_{\texttt{++}}=$
$\displaystyle\displaystyle\,\prod_{\lambda\in\mathbf{Z},\lambda^{\prime},\lambda^{\prime\prime}>0}^{(\textrm{eis})}\frac{\left(1-\frac{z}{\lambda+A\lambda^{\prime}+A^{\prime}\lambda^{\prime\prime}+v}\right)^{N}}{\left(1-\frac{Nz}{\lambda+NA\lambda^{\prime}+NA^{\prime}\lambda^{\prime\prime}+Nv}\right)}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\,\prod_{\lambda^{\prime},\lambda^{\prime\prime}>0}\frac{(1-e^{2\pi
i(A\lambda^{\prime}+A^{\prime}\lambda^{\prime\prime}+v-z)})^{N}(1-e^{2\pi
i(NA\lambda^{\prime}+NA^{\prime}\lambda^{\prime\prime}+Nv)})}{(1-e^{2\pi
i(A\lambda^{\prime}+A^{\prime}\lambda^{\prime\prime}+v)})^{N}(1-e^{2\pi
i(NA\lambda^{\prime}+NA^{\prime}\lambda^{\prime\prime}+Nv-Nz)})}\cdot$

We write $\displaystyle Q_{\texttt{--}}$ for the similar ratio, now indexed by
$\displaystyle\lambda^{\prime},\lambda^{\prime\prime}\in\mathbf{Z}$ satisfying
the inequalities $\displaystyle\lambda^{\prime},\lambda^{\prime\prime}\leq 0.$
It is then straightforward to express the quotient of quotients $\displaystyle
Q_{\texttt{++}}/Q_{\texttt{--}}$ solely in terms of elliptic gamma functions
as

(5.2)
$\frac{Q_{\texttt{++}}}{Q_{\texttt{--}}}=\frac{\Gamma(z-v,A,A^{\prime})^{N}/\Gamma(-v,A,A^{\prime})^{N}}{\Gamma(Nz-
Nv,NA,NA^{\prime})/\Gamma(-Nv,NA,NA^{\prime})}\cdot$

Additionally, Eisenstein highlights as a case of special interest for Number
Theory when the inequalities correspond to grouping together values of the
indices “in geometric progression.” We are inclined following [31] to infer
that Eisenstein could have foreseen applications to the arithmetic of cubic
complex fields, for which the unit group is of rank one.

He concludes:

> Pursuing this path, the functions one is led to seem to possess very
> remarkable properties; they open up a field in which the most prolific
> investigations offer themselves and appear to be ground where the most
> difficult parts of Analysis and Number Theory join.

This cryptic prediction,“as good as lost” as many original ideas of Eisenstein
([38], p.4), predates the utterance of Kronecker’s Jugendtraum by 36 years. As
the very name indicates, Kronecker must have carried his own aspiration over a
long period of time. We hope our paper provides some evidence that their two
programs, maybe conceived while they were together students in Berlin, have
tight connections. It is however to be observed that only Eisenstein’s text
points to specific functions (“infinite triple products”) to pursue his dream
beyond the $\displaystyle\mathrm{GL}_{2}$ setting.

Considering that Eisenstein was only 21 when he outlined this rough sketch of
a program, and that the arithmetic applications we have unveiled are, in the
case of complex cubic fields, quite close to Kronecker’s Jugendtraum, or
Hilbert’s twelfth problem, we felt entitled to name it “Eisenstein’s
Jugendtraum.”

As a closing note, let us recast our initial numerical data somewhat
differently to confirm that Eisenstein’s meromorphic quotients of quotients
(5.2) undoubtedly possess striking number theoretic properties. In parallel to
(1.2), it is legitimate to choose the specializations $\displaystyle z=1/3$,
$\displaystyle v=-1/3$, $\displaystyle N=5$ together with $\displaystyle
A=-(2\beta+\beta^{2})/15$, $\displaystyle A^{\prime}=-(2+\beta)/15$ where
$\displaystyle\beta=\sqrt[3]{7}e^{-\frac{2i\pi}{3}}.$ We deduce that, even if
neither $\displaystyle Q_{\texttt{++}}$ nor $\displaystyle Q_{\texttt{--}}$
seem to be algebraic individually, their ratio

(5.3) $\frac{Q_{\texttt{++}}}{Q_{\texttt{--}}}\approx-1735.727884\ldots+i\cdot
336.859173\ldots$

does. To $\displaystyle 1000$ digits of accuracy, this is the square of the
root (1.2) of the degree six polynomial with coefficients in
$\displaystyle\mathbf{Z}[\beta]$ displayed in the Introduction.

## References

  * [1] Allan Adler.  Eisenstein and the Jacobian varieties of Fermat curves.  Rocky Mountain J. Math., 27(1):1–60, 1997.
  * [2] Avner Ash and Lee Rudolph.  The modular symbol and continued fractions in higher dimensions.  Invent. Math., 55(3):241–250, 1979.
  * [3] Hohto Bekki.  The relative Hecke integral formula for an arbitrary extension of number fields.  J. Number Theory, 197:185–217, 2019.
  * [4] Nicolas Bergeron, Pierre Charollois, and Luis E. García.  Cocycles de groupe pour GLn et arrangements d’hyperplans, volume 39 of CRM Monograph Series.  AMS and Centre de Recherches Mathématiques, first edition, 2023.
  * [5] Nicolas Bergeron, Pierre Charollois, and Luis E. García.  Eisenstein cohomology classes for GLN over imaginary quadratic fields.  J. Crelle, 797:1–40, 2023.
  * [6] Jean-Michel Bismut, Henri Gillet, and Christophe Soulé.  Analytic torsion and holomorphic determinant bundles. III. Quillen metrics on holomorphic determinants.  Comm. Math. Phys., 115(2):301–351, 1988.
  * [7] Armand Borel.  Stable real cohomology of arithmetic groups.  Ann. Sci. École Norm. Sup. (4), 7:235–272 (1975), 1974.
  * [8] Armand Borel and Lizhen Ji.  Compactifications of symmetric and locally symmetric spaces.  In Lie theory, volume 229 of Progr. Math., pages 69–137. Birkhäuser Boston, Boston, MA, 2005.
  * [9] Pierre Charollois and Henri Darmon.  Arguments des unités de Stark et périodes de séries d’Eisenstein.  Algebra Number Theory, 2(6):655–688, 2008.
  * [10] Pierre Charollois and Robert Sczech.  Elliptic functions according to Eisenstein and Kronecker: an update.  Eur. Math. Soc. Newsl., (101):8–14, 2016.
  * [11] John Coates and Andrew Wiles.  On the conjecture of Birch and Swinnerton-Dyer.  Invent. Math., 39(3):223–251, 1977.
  * [12] Henri Darmon and Samit Dasgupta.  Elliptic units for real quadratic fields.  Ann. of Math. (2), 163(1):301–346, 2006.
  * [13] Henri Darmon, Lennart Gehrmann, and Michael Lipnowski.  Rigid meromorphic cocycles for orthogonal groups, 2023.  ArXiv:2308.14433.
  * [14] Henri Darmon, Alice Pozzi, and Jan Vonk.  The values of the Dedekind–Rademacher cocycle at real multiplication points, 2021.  ArXiv:2103.02490.
  * [15] Henri Darmon and Jan Vonk.  Singular moduli for real quadratic fields: a rigid analytic approach.  Duke Math. J., 170(1):23–93, 2021.
  * [16] Samit Dasgupta.  Stark’s Conjectures, 1999.  Senior honors thesis, Harvard University, https://services.math.duke.edu/~dasgupta/papers/Dasguptaseniorthesis.pdf.
  * [17] Samit Dasgupta and Mahesh Kakde.  On the Brumer-Stark conjecture.  Ann. of Math. (2), 197(1):289–388, 2023.
  * [18] Ehud de Shalit.  Iwasawa theory of elliptic curves with complex multiplication, volume 3 of Perspectives in Mathematics.  Academic Press, Inc., Boston, MA, 1987.  $\displaystyle p$-adic $\displaystyle L$ functions.
  * [19] David S. Dummit, Brett A. Tangedal, and Paul B. van Wamelen.  Stark’s conjecture over complex cubic number fields.  Math. Comp., 73(247):1525–1546, 2004.
  * [20] Isaac Efrat.  Determinants of Laplacians and a second limit formula in $\displaystyle{\rm GL}(3)$.  Duke Math. J., 55(2):349–360, 1987.
  * [21] Isaac Efrat.  On a $\displaystyle{\rm GL}(3)$ analog of $\displaystyle|\eta(z)|$.  J. Number Theory, 40(2):174–186, 1992.
  * [22] Gotthold Eisenstein.  Bemerkungen zu den elliptischen und Abelschen Transcendenten.  J. Crelle, 27:185–191, 1844. French translation = Remarques sur les transcendantes elliptiques et abéliennes. J. Liouville, 10:445-450, 1845.
  * [23] Gotthold Eisenstein.  Beiträge zu Theorie der elliptischen Functionen.  J. Crelle, 35:137–274, 1847.
  * [24] Giovanni Felder, André Henriques, Carlo A. Rossi, and Chenchang Zhu.  A gerbe for the elliptic gamma function.  Duke Math. J., 141(1):1–74, 2008.
  * [25] Giovanni Felder and Alexander Varchenko.  The elliptic gamma function and $\displaystyle{\rm SL}(3,{\bf Z})\ltimes{\bf Z}^{3}$.  Adv. Math., 156(1):44–76, 2000.
  * [26] Carl G. Jacobi.  Note sur les fonctions abéliennes.  J. Crelle, 30:183–184, 1846.
  * [27] Kazuya Kato.  $\displaystyle p$-adic Hodge theory and values of zeta functions of modular forms.  Number 295, pages ix, 117–290. 2004.  Cohomologies $\displaystyle p$-adiques et applications arithmétiques. III.
  * [28] Serge Lang.  Elliptic functions, volume 112 of Graduate Texts in Mathematics.  Springer-Verlag, New York, second edition, 1987.  With an appendix by J. Tate.
  * [29] Varghese Mathai and Daniel Quillen.  Superconnections, Thom classes, and equivariant differential forms.  Topology, 25(1):85–110, 1986.
  * [30] The PARI Group, Univ. Bordeaux.  PARI/GP version 2.15.2, 2022, available from http://pari.math.u-bordeaux.fr/.
  * [31] Tian Ren and Robert Sczech.  A refinement of Stark’s conjecture over complex cubic number fields.  J. Number Theory, 129(4):831–857, 2009.
  * [32] Gilles Robert.  Unités elliptiques et formules pour le nombre de classes des extensions abéliennes d’un corps quadratique imaginaire.  Bull. Soc. Math. France Mém., No. 36. Société Mathématique de France, Paris, 1973.  Supplément au Bull. Soc. Math. France, Tome 101.
  * [33] Karl Rubin.  Elliptic curves with complex multiplication and the conjecture of Birch and Swinnerton-Dyer.  In Arithmetic theory of elliptic curves (Cetraro, 1997), volume 1716 of Lecture Notes in Math., pages 167–234. Springer, Berlin, 1999.
  * [34] Simon Ruijsenaars.  First order analytic difference equations and integrable quantum systems.  J. Math. Phys., 38(2):1069–1146, 1997.
  * [35] Carl Ludwig Siegel.  Advanced analytic number theory, volume 9 of Tata Institute of Fundamental Research Studies in Mathematics.  Tata Institute of Fundamental Research, Bombay, second edition, 1980.
  * [36] Harold M. Stark.  $\displaystyle L$-functions at $\displaystyle s=1$. IV. First derivatives at $\displaystyle s=0$.  Adv. in Math., 35(3):197–235, 1980.
  * [37] John Tate.  Les conjectures de Stark sur les fonctions $\displaystyle L$ d’Artin en $\displaystyle s=0$, volume 47 of Progress in Mathematics.  Birkhäuser Boston, Inc., Boston, MA, 1984.  Lecture notes edited by Dominique Bernardi and Norbert Schappacher.
  * [38] André Weil.  Elliptic functions according to Eisenstein and Kronecker.  Classics in Mathematics. Springer-Verlag, Berlin, 1999.  Reprint of the 1976 original.