# Structuring Concept Space With The Musical Circle Of Fifths By Utilizing
Music Grammar Based Activations

1st Tofara Moyo Bulawayo , Zimbabwe
<EMAIL_ADDRESS>Mazusa AI

###### Abstract

In this paper, we explore the intriguing similarities between the structure of
a discrete neural network, such as a spiking network, and the composition of a
piano piece. While both involve nodes or notes that are activated sequentially
or in parallel, the latter benefits from the rich body of music theory to
guide meaningful combinations. We propose a novel approach that leverages
musical grammar to regulate activations in a spiking neural network, allowing
for the representation of symbols as attractors. By applying rules for chord
progressions from music theory, we demonstrate how certain activations
naturally follow others, akin to the concept of attraction. Furthermore, we
introduce the concept of modulating keys to navigate different basins of
attraction within the network. Ultimately, we show that the map of concepts in
our model is structured by the musical circle of fifths, highlighting the
potential for leveraging music theory principles in deep learning algorithms.

##  I Introduction

###  I-A Music Theory

####  I-A1 Background

Pythagoras, the ancient Greek mathematician and philosopher, made significant
discoveries in the field of music that have had a lasting impact on our
understanding of the relationship between mathematics and music. One of his
most famous discoveries is the Pythagorean tuning system, which relates the
frequencies of musical notes to simple ratios of whole numbers.

In the Pythagorean tuning system, the frequency of a note is determined by the
ratio of the lengths of two vibrating strings. For example, if one string is
twice as long as another, the note produced by the longer string will have a
frequency that is half that of the note produced by the shorter string. This
simple relationship between frequency and ratio forms the basis of the musical
scale and allows for the creation of harmonious chords and melodies.

Pythagoras also discovered the concept of the musical circle of fifths, which
is a diagram that shows the relationships between different musical keys. The
circle of fifths is based on the mathematical concept of group theory, which
is a branch of mathematics that studies the properties of groups of objects
and the relationships between them.

In the circle of fifths, each key is represented by a point on the circle, and
the keys are arranged in a way that shows the relationships between them based
on the intervals of a perfect fifth. Moving clockwise around the circle, each
key is a perfect fifth higher than the previous key, creating a harmonious
progression of chords and melodies.

####  I-A2 Music Groups

Group theory is a branch of mathematics that deals with the study of groups,
which are sets equipped with an operation that combines any two elements to
produce a third element in the set. Groups have various properties, such as
closure, associativity, identity element, and inverses, that make them a
fundamental concept in mathematics.

In the context of music theory, group theory can be applied to analyze the
structure and relationships between musical elements. One interesting
application of group theory in music is the study of music groups, which are
sets of musical elements (such as notes, chords, or scales) that form a group
under a certain operation.

One example of a music group is the circle of fifths, which is a diagram that
shows the relationships between the 12 tones of the Western musical scale. The
circle of fifths is generated by the operation of moving up a perfect fifth
from one note to the next, creating a cycle that eventually returns to the
starting note.

The generators of the circle of fifths are the notes that form the basis of
the cycle, such as the note C in the key of C major. By applying group theory
concepts, we can analyze the structure of the circle of fifths and understand
how the generators define the relationships between the different keys and
chords in music.

Overall, group theory provides a powerful framework for studying the structure
and relationships in music, and the analysis of music groups such as the
circle of fifths can offer valuable insights into the underlying principles of
music theory.

##  II Creating Musical Grammar Based Activations

###  II-A Introduction

We propose a novel approach that leverages musical grammar to regulate
activations in a spiking neural network, allowing for the representation of
symbols as attractors.

By applying rules for chord progressions from music theory, we show how
certain activations naturally follow others, akin to the concept of
attraction. Furthermore,we show the concept of modulating keys to navigate
different basins of attraction within the network. Ultimately, the map of
concepts in our model is structured by the musical circle of fifths.

Additionally, we propose a method to regularize a deep neural network by
assigning a one-to-one mapping of the notes on a musical keyboard to the nodes
on a neural network layer. This allows for the calculation of perceptual
consonance in the activations of that layer.

We could train the neural network to provide an additional reward
(proportional to the level of consonance in the layers)for agents. We further
adapt the learning rule so that the value of perceptual consonance between
nodes in different layers determines the learning rate of connections between
them. Connections associated with higher perceptual consonance will have a
higher learning rate than those with lower perceptual consonance.

The equations necessary are as follows:

Let $P_{ij}$ represent the perceptual consonance between nodes $i$ and $j$ in
a neural network layer. The learning rate $\zeta_{ij}$ for the connection
between nodes $i$ and $j$ is then determined by:

$\zeta_{ij}\propto P_{ij}$

When training a deep neural network we use an update rule for changing the
weights as we learn. This can be interpreted as determining how much each
weight $w_{jk}$ contributes to the error signal by weighting the error by the
magnitude of the output activation from the previous (hidden) layer. The
gradients with respect to each $w_{jk}$ are thus considered to be the
“contribution” of that parameter to the total error signal and should be
“negated” during learning. This gives the following gradient descent update
rule for the output layer weights:

$\displaystyle w_{jk}$ $\displaystyle\leftarrow w_{jk}-\eta\frac{\partial
E}{\partial w_{jk}}$ (1) $\displaystyle\leftarrow
w_{jk}-\eta({\delta_{k}}{a_{j}})$ (2)

Where $\eta$ is the learning rate E is the error of the network and
$\delta_{k}$ can be interpreted as the network output error after being “back
propagated” through the output activation function .In order to incorporate
the learning rate $\zeta$ we need to add it to $\eta$.

We also believe that notes that are an octave apart, such as middle C and the
C above it, share contexts of firing features and should be represented by
related notes in the neural network. This could be or something that captures
relationship.

###  II-B Concept Space

In this paper we take the view that populations of activations represent
concepts. In fact a concept represents an attractor that modulates populations
until it is stable. Just as an object falling into a basin. Here a population
is a chord and basins are scales.

As we can see, with this system we are able to order the features detected by
the neural network in a very organised and regularised manner.

Because D12, the music group , has 2 generators that span the group along
different directions this makes a very complex feature representation system.
Remember that within a layer, ALL the activations are related. And there are
multiple ways in which they relate. This may be by them being a similar chord
,the same chord but different inversions, the same chord but higher extensions
, the same chord but same scale, different chord but same chord type, chord in
same scale , chordf rom parallel minor….and so forth.

Also the relationships between layers moves in the direction of the circle of
fifths and fourths compounding relationships and regularizing activations
across layers.

In fact since music theory is an unsolved theory, we have more relating
activations than just the theory of music groups to ”file” concepts.

##  III Related Works

Regularization is a technique used in machine learning and deep learning to
prevent overfitting and improve the generalization of a model. The earliest
form of regularization is known as L1 regularization, also called Lasso
regularization. L1 regularization adds a penalty term to the loss function
that is proportional to the absolute value of the weights of the model. This
encourages the model to learn sparse representations by setting some weights
to zero.

The next form of regularization is L2 regularization, also known as Ridge
regularization. L2 regularization adds a penalty term to the loss function
that is proportional to the square of the weights of the model. This
encourages the model to learn smooth representations by preventing the weights
from becoming too large.

Another form of regularization is dropout, which randomly sets a fraction of
the input units to zero during training. This helps prevent the model from
relying too heavily on any one feature and encourages it to learn more robust
representations.

Batch normalization is another regularization technique that normalizes the
input to each layer of the neural network. This helps stabilize the training
process and can improve the generalization of the model.

Orthogonal deep neural networks are a type of neural network architecture that
enforces orthogonality constraints on the weight matrices of the network. This
can help prevent the weights from becoming too correlated and can improve the
stability and generalization of the model.

Overall, regularization techniques such as L1 and L2 regularization, dropout,
batch normalization, and orthogonal deep neural networks can all help improve
the generalization of deep learning models and prevent overfitting.

The techniques in this paper describe a method of regularization. In
particular when surveying the literature our work is most closely related to
that of orthogonal deep neural networks. Our work is similar but instead of
separating features by a 90 degree angle we separate them by a 5th and a 4th
according to the music group generators.

.

The advantages we gain are that the music group has 2 non trivial generators
and we can have more complex dynamics as we pack representability into our
networks.This should lead to greater generalisability and stability in the
network.

##  IV Conclusion

In conclusion, this paper, basing off of a population theory of concepts as
attractors, shows how we can use a music grammar to represent these dynamics.
chords then form populations, and scales handle the movement between chords or
attraction. moving from one basin to the other as done in thought involves
modulating keys and moving from one scale into another. This is more involved
than has been described because music theory is as yet an unsolved theory.
Later work will function on incorporating rhythm. This may be done with
oscillatory spiking networks.More research needs to be done.

## References

  * [1] Interpretable statistical representations of neural population dynamics and geometry Adam Gosztolai, Robert L. Peach, Alexis Arnaudon, Mauricio Barahona, Pierre Vandergheynst
  * [2] Representation learning for neural population activity with Neural Data Transformers Joel Ye, Chethan Pandarinath
  * [3] The Geometry of Information Coding in Correlated Neural Populations Rava Azeredo da Silveira, Fred Rieke
  * [4] Grandmother cells and localist representations: a review of current thinking Jeffrey S. Bowers[5] A. Punnappurat
  * [5] Mathematics and group theory in music Athanase Papadopoulos (IRMA)
  * [6] Unsolved Problems in Group Theory. The Kourovka Notebook E. I. Khukhro, V. D. Mazurov
  * [7] A Group-theoretic view of Music-Dimitris ChatzisDimitris Chatzis
  * [8] [1710.10686] Regularization for Deep Learning: A Taxonomy - Authors: Kukačka, Jan, Golkov, Vladimir, Cremers, Daniel
  * [9] [2208.04413] An automatic L1-based regularization method for the analysis of FFC dispersion profiles with quadrupolar peaks Authors: Landi, Germana, Spinelli, Giovanni V., Zama, Fabia…
  * [10] [2108.05623] Existence, Stability and Scalability of Orthogonal Convolutional Neural Networks
  * [11] Orthogonal-Nets: A Large Ensemble of 2D Neural Networks for 3D Brain Tumor Segmentation Authors:Kamlesh Pawar, Shenjun Zhong, Dilshan Sasanka Goonatillake, Gary Egan, Zhaolin Chen
  * [12] Heaps’ Law and Vocabulary Richness in the History of Classical Music Harmony-Authors: Serra-Peralta, Marc, Serrà, Joan, Corral, Álvaro


# Learning Instance-Specific Parameters of Black-Box Models Using
Differentiable Surrogates

Arnisha Khondaker
University of Alberta
Edmonton, AB, Canada
<EMAIL_ADDRESS>Nilanjan Ray
University of Alberta
Edmonton, AB, Canada
<EMAIL_ADDRESS>
###### Abstract

Tuning parameters of a non-differentiable or black-box compute is challenging.
Existing methods rely mostly on random sampling or grid sampling from the
parameter space. Further, with all the current methods, it is not possible to
supply any input specific parameters to the black-box. To the best of our
knowledge, for the first time, we are able to learn input-specific parameters
for a black box in this work. As a test application we choose a popular image
denoising method BM3D as our black-box compute. Then, we use a differentiable
surrogate model (a neural network) to approximate the black-box behaviour.
Next, another neural network is used in an end-to-end fashion to learn input
instance-specific parameters for the black-box. Drawing inspiration from the
work of Tseng et al. [37], we applied our method to the Smartphone Image
Denoising Dataset (SIDD) for image denoising. The results are compelling,
demonstrating a significant increase in PSNR and a notable improvement in SSIM
nearing 0.93. Experimental results underscore the effectiveness of our
approach in achieving substantial improvements in both model performance and
optimization efficiency. For code and implementation details, please refer to
our GitHub repository here.

##  1 Introduction

Black box models are systems or processes whose internal workings are unknown
or inaccessible. These models are widely utilized across diverse domains,
including image processing units (ISPs) in consumer-grade cameras, specialized
imaging hardware, and simulations employed for modeling intricate real-world
scenarios or generating synthetic data for machine learning tasks [31, 28, 30,
2]. 111Preprint. Under Review. The configuration parameters of these black box
models often exhibit complex interactions with the output, playing a crucial
role in determining efficacy and enhancing performance. Traditionally, manual
tuning of these parameters has been labor-intensive, requiring extensive
domain expertise. While effective, a manual method is neither scalable nor
consistent. Automating the parameter tuning process can significantly reduce
time and effort, leading to more efficient model configurations. However, a
challenge arises because many black boxes are opaque and contains complex,
non-differentiable components.

Several methods exist for optimizing such black box models, depending on the
availability of gradients [17]. Efficient non-differentiable black box
optimization remains an open research problem with various solutions. For low-
dimensional problems or scenarios where fast function evaluation is feasible,
grid search remains a practical strategy. Early methods like Powell’s
bisection search [25] and the Nelder-Mead method [21] have also been employed.
Random search methods have been proposed for parameter optimization problems
involving mixed categorical and continuous parameters [3], building
probabilistic priors over function evaluations to approximate the objective
landscape. Stochastic gradient estimators [20] such as REINFORCE [39] have
been utilized to estimate gradients for non-differentiable functions [6, 35]
and subsequently perform gradient-based optimization.

More recent approaches leverage surrogate or proxy models to approximate and
optimize black box performance, mitigating the high variance associated with
score function gradient estimators and efficiently solving high-dimensional
parameter optimization problems [37, 33].

Bayesian optimization methods [4, 32, 34, 36] differ in their surrogate
modeling approaches but operate similarly in practice, though scaling to
higher-dimensional parameter spaces remains challenging. Evolutionary
algorithms [24, 11, 18] can scale to higher-dimensional spaces but require
efficient objective evaluations, which limits their use with large training
datasets.

The Local Generative Surrogate Optimization (L-GSO) method is designed for
gradient-based optimization of black-box simulators using differentiable local
surrogate models [33]. This approach is particularly beneficial in fields such
as physics and engineering, where processes are often modeled with non-
differentiable simulators that have intractable likelihoods. L-GSO employs
deep generative models to iteratively approximate the simulator within local
neighborhoods of the parameter space. This facilitates the approximation of
the simulator’s gradient, enabling gradient-based optimization of the
simulator parameters. L-GSO has proven efficient in high-dimensional parameter
spaces, especially when the parameters reside on a low-dimensional
submanifold. This approach demonstrates faster convergence to minima compared
to baseline methods such as Bayesian optimization [19], numerical optimization
[22], and score function gradient estimators [40].

In this context, Tseng et al. have pioneered the use of differentiable proxies
for parameter optimization in black-box image processing [37]. Their method
automates the parameter tuning process in both hardware and software image
processing pipelines. By employing convolutional neural networks (CNNs) to
generate differentiable proxies, they establish a mappable relationship
between configuration parameters and evaluation metrics. This innovative
technique enables efficient high-dimensional parameter search using stochastic
first-order optimizers, eliminating the need to explicitly model the intricate
lower-level transformations in image processing. The authors demonstrate that
by optimizing with this approach, traditional algorithms can outperform recent
deep learning methods on specific benchmarks [9, 10, 7].

While existing methods demonstrate the capability to tune black-box
parameters, they lack the ability to deliver input-specific parameter tuning.
Inspired by Tseng et al.’s work [37], this work automates input-specific
black-box model parameter optimization using a differentiable surrogate.
Unlike [37], we adopt UNETR architecture [12] as our surrogate CNN to model
the black box and extend the method with the following contributions:

  1. 1.

End-to-end training with the black-box in the loop to dynamically update its
parameters to generate training instances, thereby improving the surrogate
approximator’s adaptivity and performance.

  2. 2.

Using a predictive model to learn black-box parameters specific to individual
inputs.

  3. 3.

Demonstrating improved image denoising performance on the SIDD dataset [1]
using the BM3D denoiser [7] compared to other proxy-based optimization
methods.

Section 2 discusses a brief background of the foundational concepts drawn from
existing literature and research relevant to our work, while Section 3
presents our approach. Section 4 details the implementation setup and results,
and Section 5 concludes with future work.

##  2 Background

###  2.1 Black-Box Approximation

The work of Randika et al. aims to optimize the parameters of an Optical
Character Recognition (OCR) preprocessor by approximating the gradient of the
OCR engine, treating it as a black-box due to the inaccessibility of its
internal mechanisms and the potential presence of non-differentiable
components [27]. It stands in contrast to previous OCR-agnostic preprocessing
techniques, such as Otsu’s binarization [23], skeletonization [16],
geometrical correction methods [5], and deep learning-based methods like
DeepOtsu [13].

The approach involves an alternating optimization process, where noise is
added to the input in an inner loop to perturb the OCR engine and accumulate
errors between the OCR output and an approximating neural network. This neural
network is trained to approximate the OCR’s gradient, enabling the
optimization of the preprocessor parameters. The algorithm alternates between
minimizing the error of the preprocessor and the approximating neural network,
thereby enhancing OCR performance through effective gradient approximation and
parameter tuning for the preprocessor. This method demonstrates improved OCR
accuracy on challenging datasets [14, 15].

Our present setup uses the same alternating differentiable approximating
mechanism. While Randika et al. [27] uses a CRNN model for approximation, we
use a UNETR architecture.

###  2.2 Black-Box Parameter Optimization

In optimizing black-box imaging systems’ parameters using differentiable
proxies, Tseng et al.[37] propose Algorithm 1. First, a Differentiable Proxy
Function (modelled with a CNN) is trained to emulate the behavior of a given
Image Signal Processor (ISP) on input images with its parameters. This
training phase minimizes loss between proxy’s output and the actual ISP output
across multiple images and parameter samples. Then, the optimized weights
obtained from the first stage are applied to instantiate a relaxed parameter
optimization problem which involves minimizing a task-specific loss function
over the ISP parameters, utilizing the proxy’s partial derivatives with
respect to the parameters. This framework allows for gradient-based
optimization methods, enhancing scalability across larger parameter spaces and
nonlinearities in black-box ISP systems. As we build upon these insights,
concepts derived from [37] for our method are explained in sections 2.3 and
2.4.

###  2.3 Modeling the Black-Box System

The black-box model is conceptualized as abstract function $f_{\text{bb}}$
that map an input $x$ to an output $y$. These models are parameterized by a
set of parameters $\theta$, such that $y=f_{bb}(x;\theta)$. Our objective is
to address the parameter optimization problem, which can be formulated as
follows:

$\theta^{*}=\underset{\theta}{\operatorname{arg\,min}}\sum_{i=1}^{N}\mathcal{L}\left(f_{\text{bb}}(x_{i};\theta),y_{i}\right)$
(1)

where $\theta^{*}$ represents the optimal set of parameters, $x_{i}$ the
$i$-th input instance, $y_{i}$ the corresponding target output,
$f_{\text{bb}}(x_{i};\theta)$ the predicted output,
$\mathcal{L}\left(\cdot,\cdot\right)$ the task-specific loss function, and $N$
the number of target instances.

Further, we want to learn an instance-specific parameter learner
$\theta_{i}=f_{\text{pl}}(x_{i};\psi)$ parameterized by $\psi$:

$\psi^{*}=\underset{\psi}{\operatorname{arg\,min}}\sum_{i=1}^{N}\mathcal{L}\left(f_{\text{bb}}(x_{i};f_{\text{pl}}(x_{i};\psi)),y_{i}\right)$
(2)

with the hope that the instance-specific learning of black-box parameters
would be more effective.

###  2.4 Training an Approximator Model to Learn the Black-Box Function

Typically, the function $f_{\text{bb}}$ of the black box model is unknown or
non-differentiable, complicating the solution of Eq. (1) and (2). This
challenge can be addressed using a differentiable approximation
$f_{\text{am}}$ for $f_{\text{bb}}$. The approximation model, similar to
$f_{\text{bb}}$, maps inputs to outputs and also accepts parameters $\theta$
as inputs. $f_{\text{am}}$ is parameterized by its weights $\omega$ and
ensures it is differentiable with respect to $\theta$ and $\omega$. The
optimization problem for the approximator model is:

$\omega^{*}=\underset{\omega}{\operatorname{arg\,min}}\sum_{i=1}^{N}\underset{\theta}{\mathbb{E}}~{}\mathcal{L}(f_{\text{am}}(x_{i},\theta;\omega),f_{\text{bb}}(x_{i};\theta))$
(3)

where the loss $\mathcal{L}$ ensures that the approximator closely resembles
the output of the black-box model. Since $f_{\text{am}}$ is differentiable
with respect to $\omega$, Tseng et al. [37] used first-order techniques to
mimic black-box outputs and obtain optimized weights $\omega^{*}$ (see
Algorithm 1) . This optimization simplifies the black-box parameter
optimization problem to:

$\theta^{*}=\underset{\theta}{\operatorname{arg\,min}}\sum_{i=1}^{N}\mathcal{L}(f_{\text{am}}(x_{i},\theta;\omega^{*}),y_{i})$
(4)

Here, substituting $f_{\text{bb}}$ with $f_{\text{am}}$ from Eq. (1)
transforms the problem into a more tractable form, allowing evaluation of
$f_{\text{am}}$ and its gradient with respect to the parameters. This allows
the application of first-order optimization methods as shown in Algorithm (1).
Note that in Algorithm (1) Tseng et al. [37] needed to create an instance of
parameter vector $\theta_{i}$ for every input instance $x_{i}$. This was done
by randomly sampling $\theta$. We implemented Algorithm (1) using mini-batches
from the training set.

Algorithm 1 Black-Box Optimization

Dataset $D=\\{(x_{i},y_{i},\theta_{i})\\}_{i=1}^{N}$, initial surrogate model
weights $\omega$, number of iterations $\tau$, learning rate $\alpha>0$

Surrogate $f_{\text{am}}$ with optimized weights $\omega^{*}$ and optimum
black-box parameters $\theta^{*}$

for $t\leftarrow 0$ to $\tau-1$ do

$L_{a}\leftarrow\sum_{i}\mathcal{L}(f_{\text{am}}(x_{i},\theta;\omega),f_{\text{bb}}(x_{i};\theta))$

$\omega\leftarrow\omega-\alpha\nabla_{\omega}L_{a}$

end for

Initialize $\theta$

for $t\leftarrow 0$ to $\tau-1$ do

$L_{s}\leftarrow\sum_{i}\mathcal{L}(f_{\text{am}}(x_{i},\theta;\omega),y_{i})$

$\theta\leftarrow\theta-\alpha\nabla_{\theta}L_{s}$

end for

return Learned Black-Box parameters $\theta$

##  3 Proposed Method

Figure 1: Illustration of instance-specific parameter learning in BM3D
denoising application (Algorithm 3). The surrogate model $f_{\text{am}}$ takes
the noisy image and parameters as input channels (Figure 2). The loss $\omega$
is updated using the model output and BM3D output. The parameter learner
$f_{\text{pl}}$ outputs a parameter layer, which in the case of BM3D is
$H\times W\times 5$. This layer is averaged across each channel and quantized
back into the parameter space to obtain the BM3D output with updated
parameters. The input to the parameter learner is noisy images, and $\psi$ is
updated using the approximator output and ground truth images.

###  3.1 Dynamic Optimization with Differentiable Surrogate

Algorithm 2 Dynamic Black-Box Optimization

Dataset $D=\\{(x_{i},y_{i},\theta_{i})\\}_{i=1}^{N}$, initial surrogate model
weights $\omega$, number of iterations $\tau$, learning rate $\alpha>0$

Surrogate $f_{\text{am}}$ with optimized weights $\omega^{*}$ and optimum
black-box parameters $\theta^{*}$

for $t\leftarrow 0$ to $\tau-1$ do

$L_{a}\leftarrow\sum_{i}\mathcal{L}(f_{\text{am}}(x_{i},\theta;\omega),f_{\text{bb}}(x_{i};\theta))$

$\omega\leftarrow\omega-\alpha\nabla_{\omega}L_{a}$

$L_{s}\leftarrow\sum_{i}\mathcal{L}(f_{\text{am}}(x_{i},\theta;\omega),y_{i})$

$\theta\leftarrow\theta-\alpha\nabla_{\theta}L_{s}$

end for

return Learned Black-Box parameters $\theta$

Algorithm (1) uses a differentiable surrogate function $f_{\text{am}}$ to
model $f_{\text{bb}}$ using randomly sampled parameters of the black-box. Here
we propose a dynamic variation of Algorithm (1) following Ayantha et al. [27]
(see Algorithm 2). The important difference between Algorithm (1) and (2) is
that in every iteration the surrogate function dynamically approximates the
black-box and simultaneously updates the black-box parameter vector. As in
Algorithm (1), we used mini-batching to implement gradient-based optimization,
e.g., the Adam optimizer. Gradient update with mini-batch for $\omega$ is
followed by gradeient mini-batch update of $\theta$. For the image denosing
application, we demonstrate that the dynamic version (Algorithm 2) yields
significantly better accuracy for BM3D with more effective parameter tuning.
Note also that we used proper quantization of continuous $\theta$ before
passing to $f_{\text{bb}}$ when needed, but we pass continuous $\theta$ to the
approximator $f_{\text{am}}$.

###  3.2 Learning Parameters of Black-Box

Algorithm 3 Learning Black-box Parameters

Dataset $D=\\{(x_{i},y_{i})\\}_{i=1}^{N}$, initial surrogate model weights
$\omega$, initial parameter learner weights $\psi$, number of iterations
$\tau$, learning rate $\alpha>0$

Surrogate $f_{\text{am}}$ with optimized weights $\omega^{*}$ and parameter
learner $f_{\text{pl}}$ with optimized weights $\psi$

for $t\leftarrow 0$ to $\tau-1$ do

$\theta_{i}\leftarrow f_{\text{pl}}(x_{i};\psi),~{}\forall i$

$L_{a}\leftarrow\sum_{i}\mathcal{L}(f_{\text{am}}(x_{i},\theta_{i};\omega),f_{\text{bb}}(x_{i};\theta_{i}))$

$\omega\leftarrow\omega-\alpha\nabla_{\omega}L_{a}$

$L_{s}\leftarrow\sum_{i}\mathcal{L}(f_{\text{am}}(x_{i},f_{\text{pl}}(x_{i};\psi);\omega),y_{i})$

$\psi\leftarrow\psi-\alpha\nabla_{\psi}L_{s}$

end for

return Parameter learner $f_{\text{pl}}$

Futher to Algorithm (2), we integrate an additional model, a parameter learner
$f_{\text{pl}}$ to learn input instance-specific black-box parameters (see Eq.
(2)). This model takes as input the current instance $x$ and learns to predict
instance-specific optimal parameters that maximize the performance of the
black-box model $f_{\text{bb}}$ for that specific instance. This is
illustrated in Algorithm (3). Our experiments show that learning of instance-
specific black-box parameter further enhances the accuracy of BM3D.

##  4 Experimental Setup and Results

Figure 2: Architecture of the Surrogate Model($f_{\text{am}}$), UNETR. Output
sizes are given for patch resolution = 16, filters = [32 , 64, 128, 256],
attention heads = 12 and embedding size = 78. Parameter layer is concatenated
with input as well as hidden layers 3, 6 and 9.

###  4.1 Dataset

We employed the Smartphone Image Denoising Dataset (SIDD) developed by
Abdelhamed et al. [1] in our experiments. Specifically, we focused on 23 scene
instances captured using an iPhone 7 under normal lighting conditions with ISO
settings greater than 200. From these instances, we extracted 150 random crops
of size 512 × 512 pixels, yielding a total of 6,900 images. This dataset
includes 3,450 pairs of ground truth (GT) and noisy images. The dataset
components are as follows( see Figure 3):

  * •

Unprocessed Input Images ($I$): Raw noisy images serving as inputs.

  * •

BM3D Output with Parameters ($y_{\text{BM3D}}+\theta$): Images denoised by
BM3D after applying parameters $\theta$.

  * •

Ground Truth Images ($I^{\text{gt}}$): Corresponding denoised ground truth
images for each raw input.

For training and validation, we used 5,514 images (2,757 GT-NOISY pairs). The
training set consists of 2,413 image pairs (approximately 87.5%), while the
validation set comprises of 344 image pairs (approximately 12.5%). For
testing, 1,386 images (693 GT-NOISY pairs) were utilized to evaluate the
denoising performance.

(a)

(b)

(c)

Figure 3: Each sample of data comprises three key elements: (a) the initial
unprocessed noisy image, (b) the image processed using the BM3D algorithm with
specified parameters, and (c) the corresponding ground truth image with low
noise.

###  4.2 Evaluation Metrics

We evaluated our proposed method using Peak Signal-to-Noise Ratio (PSNR)[26]
and Structural Similarity Index Measure (SSIM)[38], metrics for measuring the
quality of an image or signal. Higher values of PSNR and SSIM indicate better
image quality.

PSNR measures the ratio between the maximum possible power of a signal and the
power of corrupting noise and is defined as:

$10\cdot\log_{10}\left(\frac{1}{N}\sum_{i=1}^{N}\left(y_{\text{output}}[i]-y_{\text{gt}}[i]\right)^{2}\right)^{-1}$
(5)

In this equation, $y_{\text{output}}$ represents the output image,
$y_{\text{gt}}$ is the ground truth image, and $N$ is the number of pixels in
the images.

On the other hand, SSIM assesses the similarity between two images and is
defined as:

$\frac{(2\mu_{x}\mu_{y}+c_{1})(2\sigma_{xy}+c_{2})}{(\mu_{x}^{2}+\mu_{y}^{2}+c_{1})(\sigma_{x}^{2}+\sigma_{y}^{2}+c_{2})}$
(6)

where $x$ and $y$ are the images being compared, $\mu_{x}$ and $\mu_{y}$ are
the average pixel values of $x$ and $y$, $\sigma_{x}^{2}$ and $\sigma_{y}^{2}$
are the variances of $x$ and $y$, $\sigma_{xy}$ is the covariance of $x$ and
$y$, and $c_{1}$ and $c_{2}$ are constants to stabilize the division.

We compared our proposed methods with Algorithm (1) as well as with random and
exhaustive search to do a thorough performance test. The PSNR and SSIM values
for each method are summarized in Table (2) and (3).

Figure 4: Architecture of the Parameter Learner ($f_{\text{pl}}$), U-Net[29].
Output sizes are given for filters = [64, 128, 256, 512]. The model takes in
noisy image $(H\times W\times 3)$ and predicts parameter layer $(H\times
W\times 5)$.

###  4.3 Black Box: BM3D

In our experiments, we utilize the BM3D denoiser—a 3D block-matching algorithm
known for effectively reducing noise in images—as our black box model [7].
Consequently, Eq. (2) is formulated as:

$\psi^{*}=\underset{\psi}{\operatorname{arg\,min}}\sum_{i=1}^{N}\mathcal{L}\left(f_{\text{BM3D}}(I_{i};f_{\text{pl}}(I_{i};\psi)),I^{\text{gt}}_{i}\right)$
(7)

Here, $f_{\text{BM3D}}(I_{i};f_{\text{pl}}(I_{i};\psi))$ denotes the output of
the BM3D denoiser applied to input image $I_{i}$ with predicted parameters,
and $I_{i}^{\text{gt}}$ denotes the corresponding ground truth clean image.

The key tuning parameters (see Table 1) of BM3D include patch size ($n1$),
color space ($cspace$), Wiener transform type ($wtransform$), linear scale for
noise ($cff$), and neighborhood size ($neighborhood$). These parameters are
systematically optimized to enhance denoising performance, aiming for metrics
that closely match ground truth images.

Parameter Name | Search Range
---|---
cff | Continuous in [1, 20]
n1 | Discrete in { 4, 8}
cspace | Discrete in {0, 1}
wtransfor | Discrete in {0, 1}
neighborhood | Discrete in { 3, 4, 5, …, 15 }
Table 1: BM3D Parameters Search Ranges

###  4.4 Algoritm 1: Black-Box Optimization

For Algorithm (1), we utilize UNETR (U-Net Transformers) (Figure 2) to train a
differentiable surrogate model instead of U-Net [29] used in the original work
for fair comparison with our method. UNETR integrates the Vision Transformer
(ViT) [8] with a U-Net architecture [29], combining a transformer-based
encoder with a U-shaped architecture that excels at learning sequence
representations from input volumes and capturing essential global and multi-
scale information [12].

To minimize weights $\omega$ of the surrogate model, we utilize two components
from the dataset: unprocessed images $I$ and outputs from BM3D along with its
parameters $\theta$. During training, these inputs are fed into UNETR to
approximate BM3D’s output and functional behavior. The approximation stage
uses an eight-channel input composed of three image channels and five
parameters (treated as additional channels). This input is processed by the
ViT to extract global features using self-attention mechanisms. The model then
transitions into a U-Net-like structure with encoder and decoder layers.

Encoders’ filter sizes increase from 32 to 64, 128, and 256, while the
decoders reverse this operation. Transformer hidden states 3, 6 and 9 are
integrated at the first 3 stages of the encoder-decoder structure (see Figure
2). Before passing them onto the encoders, parameter channels are concatenated
with these hidden states. The model is trained with a learning rate of 0.005
using ADAM optimizer and a batch size of 1, achieving convergence around epoch
20. Note that, using U-Net as in the original work has a slower convergence
speed where it occurs over epoch 30.

Upon obtaining the approximator UNETR, the 5-channel parameter layer is
optimized using ADAM optimizer with a Mean Squared Error (MSE) loss function.
The initial learning rate is set to 0.02 and is reduced every 2 epochs
according to the schedule $\text{LR}_{i+1}=0.8\times\text{LR}_{i}$. The
parameters are adjusted until convergence, which typically occurs rapidly,
usually within 8 epochs.

(a)

(b)

Figure 5: Illustration of the losses during Algorithm 1 (a) shows the
surrogate approximator loss, whereas (b) shows the parameter optimization loss
over epochs.

###  4.5 Algoritm 2: Dynamic Black-Box Optimization

For Algorithm (2), we also employ UNETR to train a differentiable surrogate
model, similar to Algorithm (1). The surrogate loss $L_{a}$ is computed based
on the BM3D outputs given current optimized parameters, and the surrogate
model weights $\omega$ are updated accordingly. Once a sufficiently accurate
approximation of the black-box model is achieved, typically around 10-20
iterations of training, we proceed to optimize the parameters $\theta$. The
parameters are optimized with ADAM optimizer using MSE loss and a learning
rate of 0.02. After optimization, the tuned parameter layer is averaged across
each channel and quantized for input into BM3D. Subsequently, the BM3D
processed images with the updated $\theta$ are used for calculating the loss
of the approximator for next iteration, maintaining an end-to-end training
approach. The non-quantized parameter layer is concatenated with the input and
hidden states 3, 6 and 9 for subsequent iterations. The convergence plot for
this training is illustrated in Figure 6.

(a)

(b)

Figure 6: Illustration of the losses during Algorithm 2 (a) shows the
surrogate approximator loss, whereas (b) shows the parameter optimization loss
over epochs.

###  4.6 Algorithm 3: Learning Black-Box Parameters

In addition to using UNETR for approximation, in Algorithm 3 we use a U-Net
architecture (Figure 4) to determine the optimal set of parameters ${\psi}$
for each input image, $I_{i}$. The U-Net model is trained using Mean Squared
Error (MSE) loss and optimized with the Adam optimizer, employing a learning
rate set at 0.02. Our U-Net architecture initiates with convolutional layers
featuring 64 filters, progressively scaling up to 512 filters (see Figure 4),
following the original design of U-Net [29]. The input to the parameter
learner, $f_{\text{pl}}$ is the noisy image of dimension $H\times W\times 3$,
and it produces an output of dimensions $H\times W\times 5$. This learner is
integrated into the training loop in an end-to-end fashion, as depicted in
Figure 1. The parameter processing is analogous to Algorithm 2, and the
surrogate model is also trained similarly.

The approximation loss typically converges around epoch 60, while the
optimization loss stabilizes by epoch 4. Figure 7 shows detailed convergence
trajectories of both the approximation and optimization losses for Algorithm
3.

(a)

(b)

Figure 7: Illustration of the losses during Algorithm 3 (a) shows the
approximation loss, whereas (b) shows the parameter learner optimization loss
over epochs. Figure 8: Scatter plot of Algorithm 1, Algorithm 2, and Algorithm
3 performance on the test samples based on PSNR per image.

###  4.7 Results

Our experiments were conducted on a system equipped with an NVIDIA GeForce RTX
4090 GPU, offering 24 GiB of memory. We explored baseline methods such as
exhaustive and random search. For the exhaustive search, we evenly sampled the
continuous parameter ($cff$) and treated other parameters as discrete. Due to
the exponential complexity of the problem, particularly with slow BM3D calls,
completing this experiment would take several months given the large parameter
space and number of images. Random search proved feasible by randomly sampling
from the parameter space and iteratively optimizing parameters based on the
PSNR metric through BM3D. Each set of 10 iterations in random search took
approximately 72 hours. Algorithms (1), (2), and (3) were run for the same
duration. During inference, we use BM3D with optimized parameters.

Method | Train | Test
---|---|---
Proxy-Opt. (Algorithm 1) | 33.41 | 33.23
Random Grid Search | 36.19 | 36.0
Dynamic-Opt. (Algorithm 2) |  36.31$\Uparrow$ |  36.12$\Uparrow$
Instance-Specific Opt. (Algorithm 3) |  36.80 $\Uparrow$ |  36.51$\Uparrow$
Table 2: PSNR values obtained by different methods on the training and test
sets. Algorithm 1, Algorithm 2 and Algorithm 3 all use UNETR as the surrogate.

The performance metrics against other baseline are shown in Table 2 and Table
3. Our method demonstrates superior performance compared to Algorithm 1 [37],
achieving PSNR values of 36.31 and 36.12 on the training and test sets,
respectively. Random search also performs competitively, yielding a PSNR
around 36.19 on train and 36.0 on test set. Transitioning from Algorithm 2 to
Algorithm 3 shows a slight improvement, with PSNR values increasing to 36.80
and 36.51, respectively. When evaluating SSIM against other baselines, our
method achieves a value of 0.929, outperforming Algorithm 1 and random search.

Method | SSIM
---|---
Proxy-Opt. (Algorithm 1) | 0.87
Random Grid Search | 0.92
Dynamic-Opt. (Algorithm 2) |  0.93 $\Uparrow$
Instance-Specific Opt. (Algorithm 3) |  0.93 $\Uparrow$
Table 3: SSIM results of different methods on the test set. Algorithm 1,
Algorithm 2 and Algorithm 3 all use UNETR as the surrogate.

(a)

(b)

(c)

Figure 9: Visualization of optimization outputs. (a) shows the optimization
output after employing Algorithm 1 [37], (b) shows our Algorithm 2, and (c)
shows 3.

The scatter plot in Figure 8 shows the test set performance based on PSNR for
individual images, with Algorithm 3 demonstrating some increase in PSNR per
instance, marked in blue.

##  5 Conclusion

In this study, we have introduced a novel approach utilizing surrogate
architecture for automated optimization of input-specific parameters in black-
box models, focusing on image denoising with the BM3D denoiser on the SIDD
dataset. Our method integrates end-to-end training with the black-box model
and parameter learning, demonstrating adaptivity and performance through
dynamic parameter updates.

Moving forward, future research directions include applying our algorithm to
diverse applications beyond image denoising, while exploring alternative
surrogate and parameter learner model architectures to further improve
accuracy and efficiency. Additionally, efforts will focus on scaling the
approach to handle larger datasets and higher-dimensional parameter spaces,
aiming to advance automated parameter optimization in complex, real-world
systems.

## References

  * [1] Abdelrahman Abdelhamed, Stephen Lin, and Michael S. Brown.  A high-quality denoising dataset for smartphone cameras.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 1692–1700. IEEE, 2018.
  * [2] Harkirat Singh Behl, Atılım Güneş Baydin, Ran Gal, Philip H. S. Torr, and Vibhav Vineet.  Autosimulate: (quickly) learning synthetic data generation.  In 16th European Conference on Computer Vision (ECCV), 2020.
  * [3] James Bergstra and Yoshua Bengio.  Random search for hyper-parameter optimization.  Journal of Machine Learning Research, 13(Feb):281–305, 2012.
  * [4] James Bergstra, Dan Yamins, and David D Cox.  Hyperopt: A python library for optimizing the hyperparameters of machine learning algorithms.  In 12th Python in Science Conference, pages 13–20. Citeseer, 2013.
  * [5] W. Bieniecki, S. Grabowski, and W. Rozenberg.  Image preprocessing for improving ocr accuracy.  In 2007 International Conference on Perspective Technologies and Methods in MEMS Design, pages 75–80. IEEE, 2007.
  * [6] Yutian Chen, Matthew W. Hoffman, Sergio Gomez Colmenarejo, Misha Denil, Timothy P. Lillicrap, and Nando de Freitas.  Learning to learn for global optimization of black box functions.  arXiv, 2016.
  * [7] Kostadin Dabov, Alessandro Foi, Vladimir Katkovnik, and Karen Egiazarian.  Image denoising by sparse 3-d transform-domain collaborative filtering.  IEEE Transactions on Image Processing, 16(8):2080–2095, 2007.
  * [8] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  CoRR, abs/2010.11929, 2020.
  * [9] Michael Elad and Michal Aharon.  Image denoising via sparse and redundant representations over learned dictionaries.  IEEE Transactions on Image Processing, 15(12):3736–3745, 2006.
  * [10] Shi Guo, Zifei Yan, Kai Zhang, Wangmeng Zuo, and Lei Zhang.  Toward convolutional blind denoising of real photographs.  arXiv preprint arXiv:1807.04686, 2018.
  * [11] Nikolaus Hansen, Sibylle D. Müller, and Petros Koumoutsakos.  Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (cma-es).  Evolutionary Computation, 11(1):1–18, 2003.
  * [12] Ali Hatamizadeh, Yucheng Tang, Vishwesh Nath, Dong Yang, Andriy Myronenko, Bennett Landman, Holger Roth, and Daguang Xu.  Unetr: Transformers for 3d medical image segmentation, 2021.
  * [13] S. He and L. Schomaker.  Deepotsu: Document enhancement and binarization using iterative deep learning.  Pattern Recognition, 91:379–390, 7 2019.
  * [14] Z. Huang, K. Chen, J. He, X. Bai, D. Karatzas, S. Lu, and C. V. Jawahar.  Icdar2019 competition on scanned receipt ocr and information extraction.  In 2019 International Conference on Document Analysis and Recognition (ICDAR), pages 1516–1520, 2019.
  * [15] Max Jaderberg, Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman.  Synthetic data and artificial neural networks for natural scene text recognition.  In Workshop on Deep Learning, NIPS, 2014.
  * [16] L. Lam and C. Y. Suen.  An evaluation of parallel thinning algorithms for character recognition.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 17:914–919, 1995.
  * [17] Gang Lei, Jianguo Zhu, Youguang Guo, Chengcheng Liu, and Bo Ma.  A review of design optimization methods for electrical machines.  Energies, 10:1962, 11 2017.
  * [18] Ilya Loshchilov, Tobias Glasmachers, and Hans-Georg Beyer.  Limited-memory matrix adaptation for large scale black-box optimization.  CoRR, abs/1705.06693, 2017.
  * [19] Jonas Mockus, Vytautas Tiesis, and Antanas Zilinskas.  Bayesian approach to global optimization.  Springer, 2(1):2–3, 1975.
  * [20] Shakir Mohamed, Mihaela Rosca, Michael Figurnov, and Andriy Mnih.  Monte carlo gradient estimation in machine learning.  arXiv preprint arXiv:1906.10652, 2019.
  * [21] John A. Nelder and Roger Mead.  A simplex method for function minimization.  The Computer Journal, 7(4):308–313, 1965.
  * [22] Jorge Nocedal and Stephen J Wright.  Numerical optimization.  Springer Science & Business Media, 2006.
  * [23] Nobuyuki Otsu.  A threshold selection method from gray-level histograms.  IEEE Transactions on Systems, Man, and Cybernetics, 9:62–66, 1979.
  * [24] Panos Y. Papalambros and Douglass J. Wilde.  Principles of Optimal Design: Modeling and Computation.  Cambridge University Press, 2000.
  * [25] MJD Powell.  A method for minimizing a sum of squares of non-linear functions without calculating derivatives.  Computational Journal, 7(4):303–307, 1965.
  * [26] William K. Pratt.  Digital Image Processing.  Wiley-Interscience, 1st edition, 1974.
  * [27] Ayantha Randika, Nilanjan Ray, Xiao Xiao, and Allegra Latimer.  Unknown-box approximation to improve optical character recognition performance.  CoRR, abs/2105.07983, 2021.
  * [28] Stephan R. Richter, Vibhav Vineet, Stefan Roth, and Vladlen Koltun.  Playing for data: Ground truth from computer games.  In European Conference on Computer Vision, pages 102–118. Springer, 2016.
  * [29] O. Ronneberger, P. Fischer, and T. Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Medical Image Computing and Computer-Assisted Intervention (MICCAI), volume 9351 of LNCS, pages 234–241. Springer, 2015.  (available on arXiv:1505.04597 [cs.CV]).
  * [30] German Ros, Laura Sellart, Joanna Materzynska, David Vazquez, and Antonio M. Lopez.  The synthia dataset: A large collection of synthetic images for semantic segmentation of urban scenes.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3234–3243, 2016.
  * [31] Nataniel Ruiz, Samuel Schulter, and Manmohan Chandraker.  Learning to simulate.  In International Conference on Learning Representations, 2019.
  * [32] Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P Adams, and Nando De Freitas.  Taking the human out of the loop: A review of bayesian optimization.  IEEE, 104(1):148–175, 2016.
  * [33] Sergey Shirobokov, Vladislav Belavin, Michael Kagan, Andrey Ustyuzhanin, and Atilim Günes Baydin.  Differentiating the black-box: Optimization with local generative surrogates.  CoRR, abs/2002.04632, 2020.
  * [34] Jasper Snoek, Hugo Larochelle, and Ryan P Adams.  Practical bayesian optimization of machine learning algorithms.  In Advances in Neural Information Processing Systems, pages 2004–2012, 2012.
  * [35] Freek Stulp and Olivier Sigaud.  Policy improvement: Between black-box optimization and episodic reinforcement learning.  Technical report, Technical Report hal-00738463, 2013.
  * [36] Kevin Swersky, Jasper Snoek, and Ryan P Adams.  Multi-task bayesian optimization.  In Advances in Neural Information Processing Systems, 2013.
  * [37] Ethan Tseng, Felix Yu, Yuting Yang, Fahim Mannan, Karl St. Arnaud, Derek Nowrouzezahrai, Jean-Francois Lalonde, and Felix Heide.  Hyperparameter optimization in black-box image processing using differentiable proxies.  ACM Transactions on Graphics (TOG), 38(4), 7 2019.
  * [38] Zhou Wang, Alan C. Bovik, Hamid R. Sheikh, and Eero P. Simoncelli.  Image quality assessment: From error visibility to structural similarity.  IEEE Transactions on Image Processing, 13(4):600–612, 2004.
  * [39] Ronald J. Williams.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  Machine Learning, 8(3-4):229–256, May 1992.
  * [40] Ronald J Williams.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  In Reinforcement Learning, pages 5–32. Springer, 1992.


# An evolutionary algorithm for online, resource constrained, multi-vehicle
sensing mission planning

Nikolaos Tsiogkas1 and David M. Lane1 1 Ocean Systems Lab, Heriot Watt
University, Edinburgh UK {nt95<EMAIL_ADDRESS>
###### Abstract

Mobile robotic platforms are an indispensable tool for various scientific and
industrial applications. Robots are used to undertake missions whose execution
is constrained by various factors, such as the allocated time or their
remaining energy. Existing solutions for resource constrained multi-robot
sensing mission planning provide optimal plans at a prohibitive computational
complexity for online application [1, 2, 3]. A heuristic approach exists for
an online, resource constrained sensing mission planning for a single vehicle
[4]. This work proposes a _Genetic Algorithm_ (GA) based heuristic for the
_Correlated Team Orienteering Problem_ (CTOP) that is used for planning
sensing and monitoring missions for robotic teams that operate under resource
constraints. The heuristic is compared against optimal _Mixed Integer
Quadratic Programming_ (MIQP) solutions. Results show that the quality of the
heuristic solution is at the worst case equal to the 5% optimal solution. The
heuristic solution proves to be at least 300 times more time efficient in the
worst tested case. The GA heuristic execution required in the worst case less
than a second making it suitable for online execution.

##  I Introduction

Robotic technology has significantly advanced during the past years. The
maturing process of the scientific breakthroughs, combined with engineering
effort, has led to reliable and affordable systems, that are becoming an
indispensable assistant to science and industry. Some examples include space
exploration [5], underwater archaeology [6] and search and rescue [7].

In many cases the robot operation is constrained by the nature of the problem
the robot is solving. For example, in emergency situations time is of essence.
It is required to gather the maximum amount of information in a limited amount
of time so that a response plan can be devised. In other cases the dynamic
nature of the environment affects the energy requirements, leaving the vehicle
with less energy to complete the required mission. For example, in a highly
tidal environment the developed currents are strong and have an effect on the
mission execution energy and time [8]. In such cases it is required to
maximise the sensing outcome with limited energy by being able to replan a
maximising trajectory online.

This paper focuses on sensing missions where specific resource constraints
need to be respected. The mission requires a team of vehicles to sample a user
defined area for scalar field estimation. A scalar field is used to represent
the physical quantities in space. Its importance is presented in [9], where
phenomena such as algal blooms, oil spills and chemical pollution are studied.
These phenomena are dynamic and they need to be monitored fast and precisely
in order to be tackled. An example can be seen in [10] with the use of
Autonomous Underwater vehicles for monitoring in the Deepwater Horizon oil
spill. For these cases [1, 2] presented a variant of the _Orienteering
Problem_ (OP) for persistent monitoring tasks using _Unmanned Aerial Vehicles_
(UAVs). The presented approach takes into account the correlation of
information among the sampling points, while respecting some budget
constraints. The problem was named the _Correlated Orienteering Problem_ (COP)
and _Mixed Integer Quadratic Programming_ (MIQP) formulations were provided
for a single (COP) and multiple vehicles (CTOP). The method was able to
generate paths for one or more vehicles that were optimal for the monitoring
task.

In [3] the COP method was adapted for an underwater sensing scenario. Results
showed its potential over the standard OP for field estimation, but its
computational complexity made it prohibitive for online use. To overcome that,
in [4], a GA based heuristic was presented. The presented heuristic was able
to generate paths for a single vehicle that were close to optimal, with low
computational cost, making it suitable to be used online.

The work presented in this paper is focused on heuristic solutions for the
multi-vehicle case of the COP. It provides a GA based heuristic method that
aims to generate close to optimal paths for the whole team, while keeping the
time complexity low. This would allow plans to be regenerated online in case
of a change in the mission parameters. For example, if some constraint
changes, or in case of a vehicle failure, a new plan can be generated for the
team, guaranteeing a close to optimal execution any time.

The rest of the paper is organised as follows. Section II presents a brief
review on previous work in the field of resource constrained optimisation. In
section III the CTOP and the proposed heuristic solution are detailed. Section
IV presents the results obtained through simulated experiments. Finally, in
section V the paper concludes discussing potential future paths for research.

##  II Previous work

In the literature various approaches to maximise the utility gathered from an
agent that operates under some resource constraint can be found. A formal
definition of the problem was first presented in [11] as the selective
travelling salesman problem. The work of [12] coined the term _orienteering
problem_ (OP), as the same problem describes the sport of orienteering.
Applications of the problem to transportation and logistics systems are
presented in [13, 14].

Extending the classic OP to multiple agents is mentioned in the literature as
the _team orienteering problem_ (TOP). Over the years multiple heuristic
solutions have been presented. In [15] a tabu search heuristic is used. Local
search methods are presented in [16] and [17]. The work of [18] presents an
ACO approach, while [19] uses simulated annealing to solve the problem.
Variable neighbourhood search is used in [20]. In [21] _variable neighbourhood
search_ (VNS) is used to tackle the Capacitated TOP with time windows. A path
relinking approach for the TOP is presented in [22]. In [23] the TOP with time
windows is solved using an evolutionary local search algorithm.

The tasks of sampling or information gathering have been modelled as instances
of the OP in various works. In [24] it is applied to lake and river
monitoring. Approximation methods of a submodular OP are used to maximise the
information gathering using autonomous surface vehicles. In [25] a recursive
greedy algorithm is used to solve the additional information gathering problem
in a sensor network. The work of [26] found the submodular OP inapplicable in
real-time tasks. It proposes the usage of a linear approximation of the
submodular OP for area coverage using a micro-UAV. A branch-and-bound solution
to the OP is used in [27] for scalar field estimation using Gaussian
processes. In [28] a GA based heuristic is used to solve an OP for marine
sampling. In [29] an agricultural monitoring application is shown using
multiple robots for sensing. The work of [30] presents a vehicle routing
problem with time windows variant that can’t serve all the customers so it is
forced to select customers based on revenue. That problem formulation is
directly related with the TOP. The problem is solved using a branch-and-price
approach. In [31] an evolutionary strategy is used for online terrain
classification using unmanned aerial vehicles. Informative path planning under
disturbances is studied in the work of [32]. They provide a dynamic
programming solution to the addressed path planning problem. In [33, 34] they
present a VNS approach for Dubins vehicles to the OP and COP respectively.
Their approach successfully computes optimised paths, but the computational
cost is prohibitive for online application. The work of [35] provides a
randomised algorithm for the informative path planning for a single vehicle.
The computational time is allowing it to be applied online. To the best of our
knowledge there is no prior work to tackle the online correlated informative
path planning under budget constraints for multiple vehicles.

##  III Methods

In this section the proposed solution and the methods that it consists of are
presented in detail. Subsection III-A briefly presents a method for scalar
field estimation using Gaussian processes. In III-B the MIQP formulation of
the CTOP is presented. Subsection III-C presents the details of the proposed
evolutionary heuristic.

###  III-A Scalar field estimation

As mentioned in section I the mission requires a scalar field to be estimated.
A method to model such a field based on a _Gaussian process_ (GP) is presented
in [27]. To achieve that, a grid of sampling points is used to approximate the
field. The field values are modelled as a multivariate Gaussian distribution,
where one dimension is used for each sampling point. The covariance is
calculated using a kernel function, such as the square exponential kernel
shown in (1).

$\displaystyle k(x_{i},x_{j})=exp(-\frac{\|x_{i}-x_{j}\|}{2l^{2}})$ (1)

The kernel is used to model the drop in correlation of two points based on
their distance. The parameter $l$ governs that behaviour and can be estimated
from previously gathered data. In general, the distance in which two points
can be correlated is larger, when $l$ grows larger. This is the correlation
taken into account in the aforementioned CTOP formulation.

###  III-B Correlated Team Orienteering Problem

The CTOP, which was first defined in [1, 2], maximises the following objective
function:

$\displaystyle\sum_{i\in V}(r_{i}x_{i}+\sum_{v_{j}\in
N_{i}}r_{j}w_{ij}x_{i}(x_{i}-x_{j}))$ (2)

For any given solution the objective function represents the sum of the reward
of each of the visited vertices in that solution. All the vertices are
included in set $V$. The variable $r_{i}$ represents the reward for visiting
vertex $i$ in a solution. To encode that a vertex $i$ is visited in a given
solution, the binary variable $x_{i}$ is used. In addition to the reward
received for visiting a vertex $i$, extra reward is obtained from its
neighbourhood $N_{i}$. This reward is calculated by summing the reward of each
unvisited neighbour $j$ multiplied by a weight $w_{ij}$. To ensure that only
reward from unvisited vertices will be added, the quadratic term
$x_{i}(x_{i}-x_{j})$ is used.

To enforce that $m$ vehicles start and finish in user defined points,
constraints (3) and (4) are used.

$\displaystyle\sum_{i\in V\setminus\\{s\\}}x_{isk}=\sum_{i\in
V\setminus\\{f\\}}x_{fik}=0,$ $\displaystyle\forall 1\leq k\leq m$ (3)
$\displaystyle\sum_{i\in V\setminus\\{s\\}}x_{sik}=\sum_{i\in
V\setminus\\{f\\}}x_{ifk}=m,$ $\displaystyle\forall 1\leq k\leq m$ (4)

In this formulation a binary variable $x_{ijk}$ is showing that a path exists
from vertex $i$ to vertex $j$ in the path of vehicle $k$. Specifically,
$x_{sik}$ shows that any vertex $i$ is visited right after the starting vertex
$s$ for the path of vehicle $k$. Likewise, variable $x_{ifk}$ defines that the
finishing vertex is after vertex $i$. For all the other vertices in the
solution constraints (5)-(9) are applied.

$\displaystyle\sum_{j\in V\setminus\\{sf\\}}x_{ijk}\leq 1,$
$\displaystyle\forall i\in V\setminus{s,f},1\leq k\leq m$ (5)
$\displaystyle\sum_{j\in V\setminus\\{sf\\}}x_{jik}\leq 1,$
$\displaystyle\forall i\in V\setminus{s,f},1\leq k\leq m$ (6)

In the solution, each vertex is visited at most once, as it is ensured from
constraints (5) and (6). This behaviour is enforced by allowing at most one
path entering and one leaving in each of the vertices in V, other than the
start and finish. This happens for all the $m$ paths.

$\displaystyle\sum_{k=1}^{m}\sum_{j\in V\setminus\\{sf\\}}x_{ijk}=x_{i}\leq
1,$ $\displaystyle\forall i\in V\setminus{f}$ (7)
$\displaystyle\sum_{k=1}^{m}\sum_{j\in V\setminus\\{sf\\}}x_{jik}=x_{i}\leq
1,$ $\displaystyle\forall i\in V\setminus{s}$ (8)
$\displaystyle\sum_{\mathclap{\j\in
V\setminus\\{sf\\}}}x_{ijk}=\sum_{\mathclap{j\in V\setminus\\{sf\\}}}x_{jik},$
$\displaystyle\forall i\in V,1\leq k\leq m$ (9)

Constraints (7) and (8) make sure that if a vertex $j$ is visited in a
solution, a path that connects $j$ with a next vertex exists. Constraint (9)
prevents a vertex having a path from one vehicle entering and the path of a
different vehicle exiting it. The maximum amount of resources used is enforced
by constraint (10).

$\displaystyle\sum_{i\in V}\sum_{j\in V}(c_{ij}+c_{i})x_{ijk}\leq C_{k},$
$\displaystyle\forall 1\leq k\leq m$ (10)

Term $c_{i}$ represents a fixed cost of performing sensing in vertex $i$. To
that it is added the cost of travelling from vertex $i$ to the next vertex
$c_{ij}$. Then the cost is multiplied by the binary variable $x_{ijk}$, used
to denote the existence of a path between $i$ and $j$ in path of robot $k$.
The sum of all these costs must be less or equal to the maximum allowed
resource usage for robot $k$. It must be noted that the fixed cost of
operations of the user defined start and finish vertices is zero.

$\displaystyle\begin{split}&u_{ik}-u_{jk}+1\leq(\left|{V}\right|-1)(1-x_{ijk})\\\
&\forall i,j\in V,i\neq j,\forall 1\leq k\leq m\end{split}$ (11)
$\displaystyle\begin{split}&0\leq u_{ik}\leq\left|{V}\right|\\\ &\forall i\in
V,\forall 1\leq k\leq m\end{split}$ (12)

Finally, constraints (11) and (12) are subtour elimination constraints,
allowing only a single tour to be generated for each robot.

###  III-C Genetic Algorithm heuristic

Genetic algorithms are one of the main meta-heuristics used for optimisation.
They try to mimic the evolution process occurring in nature by using a
population of chromosomes that undergo processes of natural selection,
reproduction and mutation as described in [36]. The genetic algorithm
presented by this paper uses the chromosome encoding as presented in [37]. It
is composed by a set vectors, or genes, each one of which is representing a
path for a single robot. In each vector the indices of the vertices to be
visited by that robot is stored. The whole chromosome, therefore, represents a
candidate CTOP solution. An outline of the procedure used during the
optimisation can be seen in algorithm 1. It can be seen that the optimisation
procedure is separated in two distinct phases. The first phase is responsible
for the population generation, while the second phase repeatedly applies the
operations of selection, crossover and mutation to optimise the population.

Algorithm 1 Genetic Algorithm for the Correlated Team Orienteering Problem

1:$TravelCosts,Rewards,MaxCosts,Start,Finish$

2:$FittestChromosome$

3:InitialisePopulation

4:while $gen\leq MaxGen$ or $FittestChromosome$ stable for 10 generations do

5: SelectNewPopulation

6: Crossover

7: Mutate

8:select $FittestChromosome$

9:return $FittestChromosome$

In the scope of this paper two different chromosome generation methods are
presented and compared. The difference is based on the way the genes,
composing the chromosome, are generated. The first method is a random method,
similar to the one presented in [4]. It can be seen in algorithm 2.

Algorithm 2 Random Gene Generation

1:$TravelCosts,Rewards,MaxCost,$ $Start,Finish,AvailVertices$

2:$Gene,AvailVertices$

3:$\textsc{Insert}(Gene,Start)$

4:$Done$ $\leftarrow$ $False$

5:while $!Done$ and $AvailVertices\neq\varnothing$ do

6: $v\leftarrow$ SelectRandom($AvailVertices$)

7: if $\textsc{CostInsert}(Gene,v)\leq MaxCost$ then

8: $\textsc{Insert}(Gene,v)$

9: $AvailVertices\leftarrow AvailVertices\setminus v$

10: else

11: $Done\leftarrow True$

12:$\textsc{Insert}(Gene,Finish)$

13:return $Gene,AvailVertices$

In the random generation method genes are created by randomly inserting
vertices until the maximum cost is reached or until there are no more free
vertices to insert. To generate a chromosome this process is repeated a number
of times, equal to the number of robots. In the chromosome generation process,
the set of available vertices is shared among the genes, so that if a vertex
is picked up by a robot, it cannot be assigned to any other robot.

The second gene generation method is named _Nearest Neighbour Randomised
Adaptive Search Procedure_ (NN-RASP). It is inspired by the _Greedy Randomised
Adaptive Search Procedure_ (GRASP) that is used in [22] to provide a heuristic
solution to the TOP. This procedure is taking into consideration the
peculiarities of the CTOP problem, regarding the correlation of information.
It and can be seen in algorithm 3.

Algorithm 3 Nearest Neighbour Randomised Adaptive Search Procedure

1:$TravelCosts,Rewards,MaxCost,Start,Finish$

2:$Gene$

3:$\textsc{Insert}(Gene,Start)$

4:$v\leftarrow\textsc{GetLast}(Gene)$

5:$N\leftarrow\textsc{GetNeighbours}(v)$

6:while $\textsc{GetSize}(N)>0$ do

7: for each $n\in N$ do

8: $\textsc{GetDistanceWeight}(n)$

9: $\textsc{GetNeighbourWeight}(n)$

10: $\textsc{GetCombinedWeight}(n)$

11: $nv\leftarrow\textsc{GetWeightedRandom}(N)$

12: if $\textsc{CostInsert}(Gene,nv)<=MaxCost$ then

13: $\textsc{Insert}(Gene,nv)$

14: $v\leftarrow\textsc{GetLast}(Gene)$

15: $N\leftarrow\textsc{GetNeighbours}(v)$

16: else

17: $N\leftarrow\\{\varnothing\\}$

18:$\textsc{Insert}(Gene,Finish)$

19:return $Gene$

This method is based on progressively constructing a path by choosing one of
the neighbours of the last inserted vertex in the path. The neighbour to be
inserted is chosen based on a categorical distribution where each one has a
probability to be chosen. This probability is based on the distance of the
neighbouring vertex from all other visited vertices in the solution and the
number of free neighbours it has. The distance from other vertices is chosen
as this can reduce overlapping of paths. The number of free neighbours is used
as it gives a better reward from correlated information, as well as, it will
allow the path to have more options to continue its construction. When a
vertex is chosen it is attempted to be inserted in the gene. If the gene that
is resulting is feasible the vertex is inserted and the search continues, else
the search stops and the gene is returned. The chromosome generation is again
performed by iteratively calling the gene generation methods until paths are
generated for all the robots.

To create the population both chromosome generation procedures are repeated
until the specified amount of chromosomes is generated. Following their
generation, the chromosomes are evaluated according to algorithm 4.

Algorithm 4 Chromosome Evaluation

1:$TravelCosts,Rewards,MaxCost$

2:$Fitness$

3:for each $g\in Genes$ do

4: $\textsc{TwoOpt}(g)$

5: $\textsc{RemoveDuplicates}(g)$

6: $GeneFit\leftarrow\textsc{GetReward}(g)^{3}/\textsc{GetCost}(g)$

7: $Fitness\leftarrow Fitness+GeneFit$

8:return $Fitness$

The fitness of each chromosome is defined as the sum of fitnesses of each gene
composing it. Each gene goes through a _2-Opt_ operation to minimise its cost.
This operation is a simple local search method that is based on swapping
segments of a path in hope of reducing the path’s cost. It is iteratively
applied to all segment combinations and stops if no improvement can be
achieved. It was first presented in [38] as a method to solve the travelling
salesman problem. Then any duplicate vertices that have been already visited
by other genes are removed. Its fitness is calculated as the total reward to
the third power divided by the gene cost as suggested in [39] for the single
agent orienteering problem.

After the population generation is completed, a repetitive process is applied.
This process is applied until either the best chromosome is stable for 10
generations or a maximum amount of generations is reached. This process
involves selecting a new population based on the previous one. A percentage of
the old population is passed to the next generation without selection. This
process is described in the literature as _Elitism_ and is known to speed up
the performance and prevent loss of good solutions during the search process
[40]. The rest of the population is selected using _Tournament selection_.
This method starts with an empty population, then repetitively chooses $n$
random individuals from the previous population and adds the fittest to the
new population, until the maximum amount of individuals is chosen.

Following the new population selection a crossover operation is performed to a
percentage of the population. The operation is presented in algorithm 5 and is
a modified version of the one presented in [37]. The crossover operation aims
to diversify the population and explore more of the search space.

Algorithm 5 Chromosome Crossover

1:$Parent_{1},Parent_{2}$

2:$Child_{1},Child_{2}$

3:$\textsc{SortGenes}(Parent_{1}$)

4:$\textsc{SortGenes}(Parent_{2}$)

5:for each $i\in{1,2}$ do

6: $tP_{1}\leftarrow Parent_{1}$

7: $tP_{2}\leftarrow Parent_{2}$

8: while $\textsc{GetNumGenes}(Child_{i})<numRobots\ \textbf{do}$

9: $p$ $\leftarrow$ SelectRandomParent($tP_{1}$, $tP_{2}$)

10: $\textsc{InsertBestGene}(Child_{i},p)$

11: RemoveVertices($tP_{1}$, $tP_{2}$)

12: EvaluateAndSortGenes($tP_{1}$, $tP_{2}$)

13: $\textsc{EvaluateChromosome}(Child_{i})$

14:return $Child_{1},Child_{2}$

For this operation to take place, two chromosomes are randomly chosen as the
parents of the crossover and they are replaced by the two produced offsprings.
Initially the genes of each parent are sorted by fitness. Then, the two
children are constructed iteratively. The first step of the iteration involves
copies of the parents to be created. Next, as long as the number of genes in
each child is less than the number of robots, the following process takes
place. A parent is selected randomly and its best gene is removed and inserted
to the child. Then, the vertices of this gene are removed from any other genes
of the temporary parents. After that the genes of the parents are evaluated
again and sorted based on their fitness. When the maximum amount of genes is
inserted in the child, it is accordingly evaluated. After both offsprings are
generated, they are returned in the place of their parents.

Continuing the population crossover, a percentage of the population is
evolving by mutating. This aims at attempting to improve the fitness of each
gene that composes the mutating chromosome and thus improve the fitness of the
chromosome itself. The mutation process resembles the one presented in [4] and
can be seen in algorithm 6.

Algorithm 6 Chromosome Mutation

1:$Chromosome,MaxCost,NumMut,AddProb$

2:$Chromosome$

3:$Genes\leftarrow\textsc{GetGenes}(Chromosome)$

4:for each $g\in Genes$ do

5: for $i\leftarrow 1,NumMut$ do

6: $p\leftarrow\textsc{GetRandom}(0,1)$

7: if $p\leq AddProb$ then

8: if $\textsc{GetCost}(g)\geq 0.95*MaxCost$ then

9: $v\leftarrow\textsc{SelectRandom}(g,1)$

10: $n\leftarrow\textsc{GetMaxFreeNeighbour}(v)$

11: $gn\leftarrow\textsc{Replace}(g,v,n)$

12: if $\textsc{Fitness}(g)\leq\textsc{Fitness}(gn)$ then

13: $g\leftarrow gn$

14: else

15: if $FV\neq[\varnothing]$ then

16: $v\leftarrow\textsc{SelectRandom}(FV,1)$

17: $idx\leftarrow\textsc{FindBestInsertion}(g,v)$

18: $\textsc{Insert}(g,v,idx)$

19: else

20: $v\leftarrow\textsc{GetMinLossVertex}(g)$

21: $g\leftarrow\textsc{Remove}(g,v)$

22:$\textsc{EvaluateChromosome}(Chromosome)$

23:return $Chromosome$

For each gene in the chromosome a total number of $NumMut$ mutations will take
place. In the case studied in this work $NumMut$ was set to $10$. A mutation
can happen in two ways, by either adding a vertex in the gene, or by removing
one. This is governed by the $AddProb$ probability and was set to be $0.9$. In
the first case there are two behaviours depending on the cost of the gene that
is mutated. If the cost is $95\%$ of the maximum cost, an improvement of the
solution is attempted by swapping a vertex of the solution with one of its
free neighbours. A random vertex in the solution is chosen and is iteratively
swapped with all its free neighbours. The swap with the maximum improvement in
the gene’s fitness is kept as the performed mutation. In the second behaviour,
a free vertex is randomly chosen. Then the best insertion place is found in
the path. It is consequently inserted if it does not violate the maximum cost
constraint. The best insertion is found using the heuristic value described in
the GRASP method in [22]. In the second case, where a removal is performed,
the minimum loss vertex for this gene’s fitness is found and removed. Likewise
with the insertion, the minimum loss vertex is found using the same heuristic
value. After all the $NumMut$ mutations happen for all the genes, the
chromosome is evaluated using the method already described and is returned.

##  IV Experimental Results

This section describes the experimental results obtained through simulations.
The planning is required for a simulated sampling mission that is to be
performed by a team of vehicles. The mission is composed by a set of sampling
points, arranged in a grid. The team is given specific starting and finishing
points. Additionally, each vehicle is allowed to use a limited amount of
budget to complete the mission. The overall goal of the team is to maximise
the utility obtained without exceeding the budget limitations.

The experiments are divided into two sections. In the first section a method
from the literature is used to tune the different parameters that govern the
performance of the heuristic algorithm. In the second section the tuned
parameters are used and proposed heuristic is compared against the optimal
MIQP solution. The two methods are compared based on the utility obtained by
the team and the time it took to generate a plan. The grid was chosen to be 9
by 9, as choosing a larger would make it impractical to collect results with
the exact method. To further limit the planning time of the exact MIQP method,
it was stopped once the solution found was less than $5\%$ different from the
optimal. These results are referred with the term $CTOP5\%$. For the
heuristic, given that it is a random process, 1000 experiments were run for
each case and statistical results are presented.

All the experiments were run on a system having an Intel i5-7600 CPU, 16GB of
RAM and running on Ubuntu 16.04. The MIQP instances were solved using the
Gurobi 7.0.2 commercial solver [41]. The heuristic was coded in C++11 and
compiled using g++ version 5.4.0. Source code for both approaches can be found
online
111https://github.com/lounick/lwga,222https://github.com/lounick/task_scheduling.

###  IV-A Algorithm Tuning

As described in section III, a GA uses a set of parameters to manage its
behaviour. A tuning process is essential to find the best parameters for the
algorithm. In the literature one can find methods for parameter tuning. In
[42] the general framework for parameter tuning is presented and various
methods are discussed. In the scope of this work a modified version of the
work presented in [43] was used.

The parameter tuning method that [43] presents is based on an evolutionary
approach, using chess rating for evaluation. Each chromosome in this method is
represented by a vector of integer and real values representing the set of
parameters. To ensure that the tuned parameters are generalising well, each
parameter set is tested against a set of problems. For tuning the presented
GA, twelve problem instances were created. They were created by altering the
problem size, the vertex distribution, and the number of vehicles. Problem
instances of 5x5, 7x7 and 9x9 were used. The vertices were distributed in a
grid and a noisy grid form. The noisy grid was created based on the canonical
grid and adding some uniform random noise to the position of each vertex. The
number of vehicles were 3 and 5. The budget was set to $75\%$ of the maximum
budget for all the problems. The rest of the approach is presented in
algorithm 7.

Algorithm 7 CRS-Tuning algorithm

1:$Problems$

2:$BestConfiguration$

3:$C\leftarrow\textsc{InitialisePopulation}$

4:for $trial\leftarrow 1,MaxTrials$ do

5: for each $c\in C$ do

6: for each $pinProblems$ do

7: for ( do$g\in 1,NumGames$)

8: $s$ $\leftarrow$ $GA(c,p)$

9: $S_{c,p,g}$ $\leftarrow$ $s$

10: EvaluateConfigurations($C,S$)

11: $parents$ $\leftarrow$ SelectParents($C$)

12: $newC$ $\leftarrow$ GenerateNewPop(parents)

13: $C$ $\leftarrow$ $newC$

14:select $BestConfiguration$

15:return $BestConfiguration$

The first step of the tuning algorithm is to generate a population of
chromosomes. Each chromosome is constructed by randomly choosing values for
the parameters from a specified range. For the tuning procedure performed for
this work, the chromosome population was composed by 100 individuals. Then the
algorithm iteratively runs for $MaxTrials$. For this work $MaxTrials$ were set
to ten. Initially, each configuration is used to run each problem $NumGames$
times and their scores are recorded. Here $NumGames$ was set to ten as well.

Following, the configurations are evaluated using the _Glicko2_ chess rating
system [44]. The evaluation is performed by each configuration comparing its
results against all the other configurations results and collecting $0,0.5$
and $1$ points for loss, draw and win respectively. These points, as well as
the opponents rating, are used to update the fitness of each configuration. In
[43] the evaluation is performed based only on the scores of the games. For
the tuning performed in this paper a bonus of $0.1$ was given to the faster
solution, while the slower solution got a penalty of $-0.1$. This was done to
differentiate the cases where two solutions were equally good utility-wise but
one was slower to compute than the other.

The next step involved selecting the parents that would generate the next
generation of configurations. The ten best configurations were selected as
parents as well as any other parents that performed close to the first ten.
The choice was performed based on their chess rating intervals. The parents
had a maximum size of half the total population. These new parents were a
portion of the new population. The rest of the population was generated by
applying uniform crossover and mutation operations with $50\%$ and $80\%$
percent respectively, as suggested by [43]. Finally the best configuration was
selected and returned.

The tuning method was run for both population generation methods presented in
section III. The resulted configurations as well as the parameter ranges can
be seen in Table I.

TABLE I: Parameter ranges and tuned values for the two different methods Parameter | Value Range | Random | NNRASP
---|---|---|---
Population Size | $[25,50,\dots,500]$ | $300$ |  $250$
Generations Number | $[5,10,\dots,50]$ | $40$ | $50$
Tournament Size | $[3,4,\dots,10]$ | $6$ | $5$
CX Probability | $[0.0,0.1,\dots,0.9]$ | $0.7$ | $0.9$
Mutation Probability | $[0.0,0.1,\dots,0.9]$ | $0.6$ | $0.7$
Elitist Percentage | $[0.01,0.02,\dots,0.20]$ | $0.19$ |  $0.03$

###  IV-B Results

The performance of the heuristic method is compared against the optimal
solution in two different ways. The first way examines how well the heuristic
method performs for different numbers of robots. For that the maximum budget
was allocated to each robot. The maximum budget was calculated by taking the
maximum budget required for a single robot to complete the mission and
dividing it by the number of robots. The results can be seen in figure 1 and
Table II.

Figure 1: Utility box plots for different numbers of robots. For lower number
of robots a sophisticated population initialisation method tends to give
better results. As the number goes higher the results come closer with a
slight superiority of the NNRASP method. The blue line represents the $5\%$
gap optimal solution obtained by solving the MIQP problem.

Figure 1 presents statistical results regarding the utility for each size of
the robotic team using the GA based planning. It is compared against a
solution produced by solving the MIQP problem using the Gurobi solver. As it
was already mentioned, the solver was stopped once the found solution was
guaranteed to be at most $5\%$ worse than optimal. This solution is
represented by the horizontal blue line. It can be seen that the GA heuristic
with the NNRASP generation method outperforms the MIQP solution in all the
cases. The random generation method falls a bit behind, with part of the
results produced being worse than the $5\%$ optimal.

TABLE II: Average time(s) for different algorithms and vehicle numbers Algorithm | Number of Robots
---|---
2 | 3 | 4 | 5
CTOP 5% | 396.36 | 695.60 | 2830.41 | 1474.18
GA-Random | 0.787703 | 0.729135 | 0.666236 | 0.623074
St.Dev. | 0.087695 | 0.0738167 | 0.0602966 | 0.0579891
GA-NNRASP | 0.792475 | 0.844305 | 0.817855 | 0.826688
St.Dev. | 0.106818 | 0.120596 | 0.112062 | 0.106136

Comparing the average planning time the GA heuristic heavily outperforms the
MIQP exact method, as it can be seen in Table II. It can be seen that the time
required for the optimal solution is at the best case 500 times higher. The
GA-Random method performed slightly better than the GA-NNRASP. Example paths
for a three robot team can be seen in figures 3(a), 3(b) and 3(c). The first
figure shows the $1\%$ optimal path calculated by the MIQP solver, while the
second is a $5\%$ optimal solution used for comparison with the GA. The
heuristic solution can be seen in the third image. It is obvious that the
solution quality of the GA is far better than the $5\%$ optimal.

Figure 2: Utility box plots for three robots having different budget. For
lower budgets a sophisticated population initialisation performs equally good
with the random one. As the budget goes higher the results show the benefits
of the NNRASP method. The blue line represents the $5\%$ gap optimal solution
obtained by solving the MIQP problem.

The second performance metric examines how the heuristic method performs in
cases where the budget is limited. For this the number of vehicles is chosen
to be three and is kept static. Then the budget of each vehicle is reduced by
a specific percentage. Figure 2 and table III present the results for this
metric.

(a) $1\%$ optimal MIQP solution for three vehicles. The gathered utility is
79.662 and the time to calculate is 5847.81 seconds.

(b) $5\%$ optimal MIQP solution for three vehicles. The gathered utility is
76.433 and the time to calculate is 690.74 seconds.

(c) GA heuristic solution for three vehicles. The gathered utility is 78.770
and the time to calculate is 0.8597 seconds.

Figure 3: Comparison of the paths generated for three vehicles having full
budget.

In figure 2 the utility gained for each budget is shown. It can be seen that
the GA method performs at least equally good as the $5\%$ exact method. That
happened only in the case where the robots had $25\%$ of the maximum energy.
In all the other cases the GA performed better. Comparing the two different
population generation methods, one can see that the more sophisticated
approach gives statistically better results. The random method performs close
to the $5\%$ optimal.

TABLE III: Average time(s) for different algorithms and different budgets Algorithm | Budget
---|---
100% | 75% | 50% | 25%
CTOP 5% | 695.60 | 267.94 | 3611.76 | 3764.12
GA-RANDOM | 0.721606 | 0.625391 | 0.362355 | 0.141678
St.Dev. | 0.0742465 | 0.0607859 | 0.0346406 | 0.0147053
GA-NNRASP | 0.837791 | 0.763285 | 0.389283 | 0.165087
St.Dev. | 0.120398 | 0.102063 | 0.0483728 | 0.0244343

Table III the timing results of the limited budget experiments are presented.
As with all the previous cases the GA-Random method performed better. In the
worst case it was more than 300 times faster than the $5\%$ MIQP solution. The
GA with NNRASP population generation method performed constantly a bit slower
than the Random generation method.

##  V Conclusions

This paper studies the online task scheduling problem for sensing missions
performed by multiple vehicles. Such problems can be solved using the
correlated team orienteering problem. The CTOP is designed for such problems
as it takes into account the correlation of sensed information between various
sensing points, while keeping the resource usage constrained. This work
presents a heuristic for solving the CTOP that is based on a genetic
algorithm. Initially, it uses a method from the literature to tune the GA
parameters. The heuristic is then compared against a mixed integer quadratic
programming solution. It is shown that in terms of utility the heuristic
performs usually better than a near optimal MIQP solution. In terms of
solution time complexity the heuristic performs at least 300 times better than
the MIQP. The average execution time for the studied instances was in the
worst case less than a second. These results make the proposed GA heuristic
suitable for online usage. Future work would involve testing the heuristic in
larger, more complex instances where the exact solution is inapplicable.
Finally, it would be worth exploring the coordination capabilities of a
heterogeneous team with different capabilities and objectives using a solution
produced by a single instance of the proposed heuristic.

## References

  * [1] J. Yu, M. Schwager, and D. Rus, “Correlated orienteering problem and its application to informative path planning for persistent monitoring tasks,” in _Intelligent Robots and Systems (IROS 2014), 2014 IEEE/RSJ International Conference on_ , Sept 2014, pp. 342–349.
  * [2] J. Yu, M. Schwager, and D. Rus, “Correlated orienteering problem and its application to persistent monitoring tasks,” _IEEE Transactions on Robotics_ , vol. 32, no. 5, pp. 1106–1118, 2016.
  * [3] N. Tsiogkas, V. De Carolis, and D. M. Lane, “Energy-constrained informative routing for auvs,” in _OCEANS 2016-Shanghai_. IEEE, 2016, pp. 1–5.
  * [4] N. Tsiogkas, V. De Carolis, and D. Lane, “Towards an online heuristic method for energy-constrained underwater sensing mission planning,” in _Intelligent Robots and Systems (IROS), 2017 IEEE/RSJ International Conference on_. IEEE, 2017.
  * [5] W. F. Truszkowski, M. G. Hinchey, J. L. Rash, and C. A. Rouff, “Autonomous and autonomic systems: A paradigm for future space exploration missions,” _IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)_ , vol. 36, no. 3, pp. 279–291, 2006.
  * [6] B. Allotta, R. Costanzi, A. Ridolfi, C. Colombo, F. Bellavia, M. Fanfani, F. Pazzaglia, O. Salvetti, D. Moroni, M. A. Pascali, _et al._ , “The arrows project: adapting and developing robotics technologies for underwater archaeology,” _IFAC-PapersOnLine_ , vol. 48, no. 2, pp. 194–199, 2015.
  * [7] T. Tomic, K. Schmid, P. Lutz, A. Domel, M. Kassecker, E. Mair, I. L. Grixa, F. Ruess, M. Suppa, and D. Burschka, “Toward a fully autonomous uav: Research platform for indoor and outdoor urban search and rescue,” _IEEE robotics & automation magazine_, vol. 19, no. 3, pp. 46–56, 2012\.
  * [8] J. Bellingham and J. Willcox, “Optimizing auv oceanographic surveys,” in _Autonomous Underwater Vehicle Technology, 1996. AUV ’96., Proceedings of the 1996 Symposium on_ , Jun 1996, pp. 391–398.
  * [9] H. M. La, W. Sheng, and J. Chen, “Cooperative and active sensing in mobile sensor networks for scalar field mapping,” _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , vol. 45, no. 1, pp. 1–12, 2015.
  * [10] J. C. Kinsey, D. R. Yoerger, M. V. Jakuba, R. Camilli, C. R. Fisher, and C. R. German, “Assessing the deepwater horizon oil spill with the sentry autonomous underwater vehicle,” in _Intelligent Robots and Systems (IROS), 2011 IEEE/RSJ International Conference on_. IEEE, 2011, pp. 261–267.
  * [11] G. Laporte and S. Martello, “The selective travelling salesman problem,” _Discrete applied mathematics_ , vol. 26, no. 2-3, pp. 193–207, 1990.
  * [12] I.-M. Chao, B. L. Golden, and E. A. Wasil, “A fast and effective heuristic for the orienteering problem,” _European Journal of Operational Research_ , vol. 88, no. 3, pp. 475–489, 1996.
  * [13] P. Vansteenwegen, W. Souffriau, and D. Van Oudheusden, “The orienteering problem: A survey,” _European Journal of Operational Research_ , vol. 209, no. 1, pp. 1–10, 2011.
  * [14] A. Gunawan, H. C. Lau, and P. Vansteenwegen, “Orienteering problem: A survey of recent variants, solution approaches and applications,” _European Journal of Operational Research_ , vol. 255, no. 2, pp. 315–332, 2016.
  * [15] H. Tang and E. Miller-Hooks, “A tabu search heuristic for the team orienteering problem,” _Computers & Operations Research_, vol. 32, no. 6, pp. 1379–1407, 2005.
  * [16] P. Vansteenwegen, W. Souffriau, G. V. Berghe, and D. Van Oudheusden, “Iterated local search for the team orienteering problem with time windows,” _Computers & Operations Research_, vol. 36, no. 12, pp. 3281–3290, 2009\.
  * [17] P. Vansteenwegen, W. Souffriau, G. V. Berghe, and D. Van Oudheusden, “A guided local search metaheuristic for the team orienteering problem,” _European journal of operational research_ , vol. 196, no. 1, pp. 118–127, 2009.
  * [18] L. Ke, C. Archetti, and Z. Feng, “Ants can solve the team orienteering problem,” _Computers & Industrial Engineering_, vol. 54, no. 3, pp. 648–665, 2008.
  * [19] S.-W. Lin and F. Y. Vincent, “A simulated annealing heuristic for the team orienteering problem with time windows,” _European Journal of Operational Research_ , vol. 217, no. 1, pp. 94–107, 2012.
  * [20] N. Labadie, R. Mansini, J. Melechovskỳ, and R. W. Calvo, “The team orienteering problem with time windows: An lp-based granular variable neighborhood search,” _European Journal of Operational Research_ , vol. 220, no. 1, pp. 15–27, 2012.
  * [21] B. Aghezzaf and H. E. Fahim, “The multi-constraint team orienteering problem with time windows in the context of distribution problems: A variable neighborhood search algorithm,” in _Logistics and Operations Management (GOL), 2014 International Conference on_. IEEE, 2014, pp. 155–160.
  * [22] W. Souffriau, P. Vansteenwegen, G. V. Berghe, and D. Van Oudheusden, “A path relinking approach for the team orienteering problem,” _Computers & operations research_, vol. 37, no. 11, pp. 1853–1859, 2010.
  * [23] N. Labadie, J. Melechovskỳ, and R. W. Calvo, “Hybridized evolutionary local search algorithm for the team orienteering problem with time windows,” _Journal of Heuristics_ , vol. 17, no. 6, pp. 729–753, 2011.
  * [24] A. Singh, A. Krause, C. Guestrin, and W. J. Kaiser, “Efficient informative sensing using multiple robots,” _Journal of Artificial Intelligence Research_ , pp. 707–755, 2009.
  * [25] J. Binney, A. Krause, and G. S. Sukhatme, “Informative path planning for an autonomous underwater vehicle,” in _Robotics and automation (icra), 2010 IEEE international conference on_. IEEE, 2010, pp. 4791–4796.
  * [26] L. Heng, A. Gotovos, A. Krause, and M. Pollefeys, “Efficient visual exploration and coverage with a micro aerial vehicle in unknown environments,” in _Robotics and Automation (ICRA), 2015 IEEE International Conference on_. IEEE, 2015, pp. 1071–1078.
  * [27] J. Binney and G. S. Sukhatme, “Branch and bound for informative path planning,” in _Robotics and Automation (ICRA), 2012 IEEE International Conference on_. IEEE, 2012, pp. 2147–2154.
  * [28] S. Frolov, B. Garau, and J. Bellingham, “Can we do better than the grid survey: Optimal synoptic surveys in presence of variable uncertainty and decorrelation scales,” _Journal of Geophysical Research: Oceans_ , vol. 119, no. 8, pp. 5071–5090, 2014.
  * [29] P. Tokekar, J. Vander Hook, D. Mulla, and V. Isler, “Sensor planning for a symbiotic uav and ugv system for precision agriculture,” in _Intelligent Robots and Systems (IROS), 2013 IEEE/RSJ International Conference on_. IEEE, 2013, pp. 5321–5326.
  * [30] N. Azi, M. Gendreau, and J.-Y. Potvin, “An exact algorithm for a vehicle routing problem with time windows and multiple use of vehicles,” _European Journal of Operational Research_ , vol. 202, no. 3, pp. 756–763, 2010.
  * [31] M. Popović, G. Hitz, J. Nieto, I. Sa, R. Siegwart, and E. Galceran, “Online informative path planning for active classification using uavs,” in _Robotics and Automation (ICRA), 2017 IEEE International Conference on_. IEEE, 2017, pp. 5753–5758.
  * [32] K.-C. Ma, L. Liu, and G. S. Sukhatme, “An information-driven and disturbance-aware planning method for long-term ocean monitoring,” in _Intelligent Robots and Systems (IROS), 2016 IEEE/RSJ International Conference on_. IEEE, 2016, pp. 2102–2108.
  * [33] R. Pěnička, J. Faigl, P. Váňa, and M. Saska, “Dubins orienteering problem,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 2, pp. 1210–1217, 2017.
  * [34] R. Pěnička, J. Faigl, P. Váňa, and M. Saska, “Dubins orienteering problem with neighborhoods,” in _Unmanned Aircraft Systems (ICUAS), 2017 International Conference on_. IEEE, 2017, pp. 1555–1562.
  * [35] S. Arora and S. Scherer, “Randomized algorithm for informative path planning with budget constraints,” in _Robotics and Automation (ICRA), 2017 IEEE International Conference on_. IEEE, 2017, pp. 4997–5004.
  * [36] M. Mitchell, _An introduction to genetic algorithms_. MIT press, 1998.
  * [37] A. Singh and A. S. Baghel, “A new grouping genetic algorithm approach to the multiple traveling salesperson problem,” _Soft Computing-A Fusion of Foundations, Methodologies and Applications_ , vol. 13, no. 1, pp. 95–101, 2009\.
  * [38] G. A. Croes, “A method for solving traveling-salesman problems,” _Operations research_ , vol. 6, no. 6, pp. 791–812, 1958.
  * [39] J. Karbowska-Chilinska, J. Koszelew, K. Ostrowski, and P. Zabielski, “Genetic algorithm solving orienteering problem in large networks.” in _KES_ , 2012, pp. 28–38.
  * [40] H. Li and D. Landa-Silva, “An elitist grasp metaheuristic for the multi-objective quadratic assignment problem,” in _Evolutionary multi-criterion optimization_. Springer, 2009, pp. 481–494.
  * [41] I. Gurobi Optimization, “Gurobi optimizer reference manual,” 2015. [Online]. Available: http://www.gurobi.com
  * [42] A. E. Eiben and S. K. Smit, “Parameter tuning for configuring and analyzing evolutionary algorithms,” _Swarm and Evolutionary Computation_ , vol. 1, no. 1, pp. 19–31, 2011.
  * [43] N. Veček, M. Mernik, B. Filipič, and M. Črepinšek, “Parameter tuning with chess rating system (crs-tuning) for meta-heuristic algorithms,” _Information Sciences_ , vol. 372, pp. 446–469, 2016.
  * [44] M. E. Glickman, “Example of the glicko-2 system,” _Boston University_ , 2012\.