<figure><img src="image_1.png" /></figure>

HAL
open science

# Algorithm Selection for Multi-Objective Optimization

Alexandre D. B. Jesus

► To cite this version:

Alexandre D. B. Jesus. Algorithm Selection for Multi-Objective Optimization. Computer Science [cs]. University of Coimbra; University of Lille, 2022. English. NNT: tel-03936989

HAL Id: tel-03936989

https://hal.science/tel-03936989

Submitted on 13 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
UNIVERSITÉ DE LILLE

École doctorale MADIS

Unité de recherche CRISTAL

# Algorithm Selection for Multi-Objective Optimization

Thèse présentée par

ALEXANDRE DANIEL BORGES DE JESUS

En vue de l'obtention du grade de docteur de l'Université de Lille

Soutenue le 9 décembre 2022 devant le jury composé de:

<table><tr><td>Président du jury</td><td>Henrique Madeira</td><td>Professeur à Universidade de Coimbra, Portugal</td></tr><tr><td>Rapporteurs</td><td>Inês Lynce<br/>Michael Emmerich</td><td>Professeure à Instituto Superior Técnico, Portugal<br/>Professeur à Universiteit Leiden, The Netherlands</td></tr><tr><td>Examineurs</td><td>Frédéric Semet<br/>Jin-Kao Hao<br/>Nuno Lourenço</td><td>Professeur à l'Université de Lille, France<br/>Professeur à l'Université de Angers, France<br/>Professeur à Universidade de Coimbra, Portugal</td></tr><tr><td>Directeurs de thèse</td><td>Bilel Derbel<br/>Luís Paquete</td><td>Professeur à l'Université de Lille, France<br/>Professeur à Universidade de Coimbra, Portugal</td></tr></table>
<figure><img src="image_3.png" /></figure>
Algorithm Selection for
Multi-Objective Optimization

by

ALEXANDRE DANIEL BORGES DE JESUS

ajesus@dei.uc.pt

A thesis submitted to the University of Coimbra in partial fulfillment of the
requirements for the doctoral degree in Information Science and Technology,
and to the University of Lille in partial fulfillment of the requirements for the
doctoral degree in Computer Science, under a co-tutelle agreement between
both institutions.

Prepared under the supervision of

LUÍS PAQUETE

Associate Professor

University of Coimbra

BILEL DERBEL

Associate Professor

University of Lille

ARNAUD LIEFOOGHE

Associate Professor

University of Lille
This thesis was financially supported by Portuguese Foundation for Science and Technology under scholarship SFRH/BD/132275/2017 co-funded by the European Social Fund and by the State Budget of the Portuguese Ministry of Education and Science. It was also partially supported by PICS project MOCO-SEARCH co-funded by the French National Center for Scientific Research and the Portuguese Foundation for Science and Technology.

This document was typeset with LATEX and BIBLATEX, in conjunction with various other TEX packages that were graciously made freely available by their authors in the Comprehensive TEX Archive Network (CTAN) and that are part of the TEX Live distribution.

© 2022 Alexandre Jesus
Abstract

Multi-objective optimization problems, which consider multiple objective func-
tions to be optimized, can arise in many real-life scenarios, e.g., when trying
to minimize both the cost and the time needed for traveling between two lo-
cations. In the last few decades, several algorithms have been proposed to
solve multi-objective optimization problems. These algorithms can have very
distinct behaviors, and their performance is often significantly affected by the
problem instance to be solved, the time budget available, or the desirable so-
lution quality. As such, which algorithm performs best often depends on the
scenario that is being considered.

This gives rise to the algorithm selection problem, which is concerned
with the automatic selection of the best algorithm for a given scenario. In
this thesis, we investigate the case of automatically selecting the best multi-
objective optimization algorithm to solve a previously unseen problem instance,
taking into account that the available time budget and desirable solution quality
may be uncertain, and are only known when selecting the algorithm. We make
several contributions in this line.

First, we propose theoretical and empirical models to characterize the any-
time performance of an algorithm, i.e., how solution quality improves over
time, for previously unseen problem instances. Then, considering these mod-
els, we develop an offline selection methodology to select the best algorithm
for a previously unseen problem instance given a utility function that describes
the desirable time budget and solution quality. We also propose an online se-
lection methodology that can swap between multi-objective branch and bound
strategies to improve anytime performance. Lastly, we propose a scalariza-
tion technique and a branch and bound search strategy for multi-objective
optimization problems to achieve a better anytime performance than previ-
ous approaches. Each contribution is backed by an experimental study on a
multi-objective knapsack problem, and the results highlight the quality of the
proposed models, selection methodologies, and algorithms.
<figure><img src="image_4.png" /></figure>
Resumo

Problemas de otimização multi-objetivo, que consideram múltiplas funções ob-
jetivo a otimizar, podem surgir em diversos cenários reais, por exemplo, quando
se quer minimizar tanto o custo como o tempo de uma viagem entre dois locais.
Nas últimas décadas, vários algoritmos foram propostos para resolver proble-
mas multi-objetivo. Estes algoritmos podem ter comportamentos distintos, e o
seu desempenho é tipicamente afetado pela instância do problema a resolver, o
tempo disponível para resolver o problema, e a qualidade da solução desejada.
Como tal, qual o algoritmo que tem melhor desempenho depende do cenário
em consideração.

Isto dá origem ao problema de seleção de algoritmos, que considera a es-
colha automática do algoritmo com melhor desempenho para um dado cenário.
Nesta tese, investigamos a seleção automática do melhor algoritmo para re-
solver instâncias do problema nunca antes vistas, tendo em conta que o tempo
disponível e a qualidade da solução desejada podem ser incertos, e apenas
conhecidos aquando da seleção. Fazemos várias contribuições nesta direção.

Em primeiro lugar, propomos modelos teóricos e empíricos para caracteri-
zar o desempenho de algoritmos *anytime*, ou seja, modelos que caracterizem a
evolução da qualidade da solução devolvida pelo algoritmo ao longo do tempo,
para instâncias do problema nunca antes vistas. Em segundo lugar, tendo em
conta os modelos propostos, desenvolvemos uma metodologia de seleção *of-
line* para selecionar o melhor algoritmo para uma instância do problema nunca
antes vista, dada uma função de utilidade que descreve o tempo disponível e
a qualidade da solução desejada. Também propomos uma metodologia de se-
leção *online* capaz de mudar de estratégias *branch and bound* multi-objetivo
de forma a melhorar o desempenho do algoritmo ao longo do tempo. Por fim,
propomos uma técnica de escalarização e uma estratégia de *branch and bound*
para problemas de otimização multi-objetivo para obter um melhor desem-
penho ao longo do tempo comparativamente a abordagens já existentes. Cada
contribuição é acompanhada por um estudo experimental de um problema
*knapsack* multi-objetivo, sendo que os resultados destacam a qualidade dos
modelos, metodologias de seleção, e algoritmos propostos.
<figure><img src="image_5.png" /></figure>
# Résumé

Les problèmes d'optimisation multi-objectifs, pour lesquels plusieurs objectifs doivent être optimisés, peuvent survenir dans de nombreux scénarios réels, par exemple lorsque l'on essaie de minimiser à la fois le coût et le temps nécessaires pour se déplacer entre deux emplacements. Au cours des dernières décennies, de nombreux algorithmes ont été proposés afin de résoudre des problèmes d'optimisation multi-objectifs. Ces algorithmes peuvent avoir des comportements très distincts et leurs performances sont souvent affectées de manière significative par l'instance du problème à résoudre, le budget temps disponible pour la résolution, ou encore la qualité de solution souhaitée. Ainsi, l'algorithme qui fonctionne le mieux dépend souvent du scénario envisagé.

Cela donne lieu au problème de sélection d'algorithme, qui concerne la sélection automatique du meilleur algorithme pour un scénario donné. Dans cette thèse, nous étudions le cas de la sélection automatique du meilleur algorithme d'optimisation multi-objectifs pour résoudre une instance de problème non-rencontrée auparavant, en tenant compte du fait que le budget temps disponible et la qualité de solution souhaitée peuvent être incertains, et ne sont connus qu'à l'étape de la sélection de l'algorithme. Nous apportons plusieurs contributions dans cette voie.

Dans un premier temps, nous proposons des modèles théoriques et empiriques pour caractériser la performance "anytime" d'un algorithme, c'est-à-dire comment la qualité de solution s'améliore au fil du temps, pour des instances de problème non-rencontrées auparavant. Ensuite, en considérant ces modèles, nous développons une méthodologie de sélection hors ligne afin de sélectionner le meilleur algorithme étant donné une fonction d'utilité qui décrit le budget temps et la qualité de solution souhaités. Nous proposons également une méthodologie de sélection en ligne qui peut basculer entre des stratégies de branch and bound multi-objectifs pour améliorer les performances "anytime". Enfin, nous proposons une technique de scalarisation et une stratégie de branch and bound pour l'optimisation multi-objectifs afin d'obtenir une meilleure performance "anytime" que les approches précédentes. Chaque contribution est étayée par une étude expérimentale sur un problème de sac à dos
multi-objectifs, et les résultats mettent en évidence la qualité des modèles, des méthodologies de sélection et des algorithmes proposés.
**Acknowledgements**

As I am working on the final details of this thesis, I cannot help but reflect on the countless hours, nights, discussions, gatherings, that have led me to this moment, and realize that I was fortunate to have so many people around to support me, both technically and emotionally, without whom this work would not have been possible.

I owe my first thanks to my supervisors, Luís Paquete, Arnaud Liefooghe, and Bilel Derbel, who were always there to guide me with their expertise and support me in this journey, who always treated me very kindly, and who helped me grow a great deal. To Bilel, thank you for all the help and advice, for making me feel a part of the team during my year in Lille, and for all the opportunities while I was there. To Arnaud, thank you for suggesting this interesting topic, and for making me feel welcome in Lille. Thank you also for always being there to answer any questions despite the distance most of the time, and thank you for always being available to discuss anything.

To Luís, I owe a great deal to you. Thank you for instilling in me the curiosity for optimization from the first day I step foot in one of your classrooms, nearly 12 years ago, and for continuously challenging me during my bachelor and masters at various stages, all of which eventually led me here. Thank you for your patience, for the countless discussions we had, for always being quick to give great feedback, for all the opportunities, and for all the pieces of advice, all of which were tremendously impactful and allowed me to grow as a researcher, as a teacher, and as a person.

I also owe a special thanks to all the people with whom I was fortunate to share an office in Coimbra and in Lille. In particular, thank you to Andreia Guerreiro, Geoffrey Pruvost, Gonçalo Lopes, João Macedo, Maryam Abbasi, Nicolas Berveglieri, Pedro Correia, Pedro Rodrigues, Stephan Helfrich, Tiago Carneiro. I am indebted to all of you for your friendship, for the great discussions and work environment, and for making this such a pleasant journey. A special thanks also to Carlos Fonseca, whose curiosity, rigor, and passion for research are contagious and have set a high standard for me.

Thank you to Luís Rodrigues and Luís Coimbra for inviting me to be a
part of NEDUC, and thank you to all NEDUC members who work to make the
University of Coimbra a better place for all its PhD students.

Thank you also to all my friends who have been there throughout this
journey and throughout my life. Thank you for always being there to help
build such great memories, and for always being there to offer a kind word and
support when needed.

However, none of this would have been possible without the support of my
loving family, my parents, my brother, and my aunt, who have always been
there to support me, to provide for me, to advise me, to teach me, and who
have always encouraged me to follow my passions.

And Fátima, my partner in everything, thank you for being my home, for
always being by my side, and for always believing in me and encouraging me
without hesitation. I am forever in your debt.

Lastly, to Coco, thank you for being a constant source of joy from the
moment you entered my life.
# Table of Contents

<table><tr><td>Abstract</td><td>iii</td></tr><tr><td>Resumo</td><td>v</td></tr><tr><td>Résumé</td><td>vii</td></tr><tr><td>Acknowledgements</td><td>ix</td></tr><tr><td>Table of Contents</td><td>xi</td></tr><tr><td>1 Introduction</td><td>1</td></tr><tr><td>1.1 Motivations</td><td>1</td></tr><tr><td>1.2 Goals & Scope</td><td>2</td></tr><tr><td>1.3 Contributions</td><td>4</td></tr><tr><td>1.4 Publications & Software</td><td>6</td></tr><tr><td>1.5 Outline</td><td>8</td></tr><tr><td>2 Background</td><td>9</td></tr><tr><td>2.1 Multi-Objective Optimization</td><td>9</td></tr><tr><td>2.2 Quality Indicators</td><td>11</td></tr><tr><td>2.2.1 Hypervolume Indicator</td><td>12</td></tr><tr><td>2.2.2 Binary Hypervolume Indicator</td><td>13</td></tr><tr><td>2.2.3 The ε-indicator</td><td>14</td></tr><tr><td>2.3 Solving the Multi-Objective Binary Knapsack Problem</td><td>15</td></tr><tr><td>2.3.1 Scalarization Techniques</td><td>16</td></tr><tr><td>2.3.2 Branch-and-Bound</td><td>18</td></tr><tr><td>2.3.3 Dynamic Programming</td><td>22</td></tr><tr><td>2.3.4 Pareto Local Search</td><td>24</td></tr><tr><td>2.4 Anytime Algorithms</td><td>25</td></tr><tr><td>2.5 Algorithm Selection</td><td>28</td></tr><tr><td>2.6 Conclusion</td><td>29</td></tr><tr><td>3 Theoretical Model of Anytime Performance</td><td>31</td></tr><tr><td>3.1 Theoretical Model of Anytime Performance</td><td>31</td></tr></table>
<table><tr><td>3.1.1</td><td>Estimating the non-dominated set</td><td>33</td></tr><tr><td>3.1.2</td><td>Analytical methodology</td><td>35</td></tr><tr><td>3.1.3</td><td>Algorithmic methodology</td><td>39</td></tr><tr><td>3.1.4</td><td>Empirical study</td><td>42</td></tr><tr><td>3.2</td><td>Guided ε-constraint</td><td>49</td></tr><tr><td>3.3</td><td>Discussion</td><td>55</td></tr><tr><td>4</td><td>Empirical Models of Anytime Performance</td><td>57</td></tr><tr><td>4.1</td><td>Frameworks for Empirical Models</td><td>58</td></tr><tr><td>4.1.1</td><td>Framework I</td><td>58</td></tr><tr><td>4.1.2</td><td>Framework II</td><td>59</td></tr><tr><td>4.1.3</td><td>Framework III</td><td>61</td></tr><tr><td>4.2</td><td>Experimental Study</td><td>62</td></tr><tr><td>4.2.1</td><td>Instances</td><td>63</td></tr><tr><td>4.2.2</td><td>Anytime Performance Measure</td><td>64</td></tr><tr><td>4.2.3</td><td>Model Implementation Details</td><td>65</td></tr><tr><td>4.2.4</td><td>Experimental Setup</td><td>68</td></tr><tr><td>4.2.5</td><td>Results</td><td>69</td></tr><tr><td>4.3</td><td>Discussion</td><td>104</td></tr><tr><td>5</td><td>Offline Algorithm Selection</td><td>107</td></tr><tr><td>5.1</td><td>Anytime Performance Measure</td><td>108</td></tr><tr><td>5.2</td><td>Anytime Algorithm Selection</td><td>110</td></tr><tr><td>5.3</td><td>Experimental Study</td><td>111</td></tr><tr><td>5.3.1</td><td>Experimental Setup</td><td>112</td></tr><tr><td>5.3.2</td><td>Results</td><td>113</td></tr><tr><td>5.4</td><td>Discussion</td><td>124</td></tr><tr><td>6</td><td>Online Algorithm Selection</td><td>127</td></tr><tr><td>6.1</td><td>Indicator-based Branch-and-Bound</td><td>128</td></tr><tr><td>6.1.1</td><td>Framework</td><td>128</td></tr><tr><td>6.1.2</td><td>Implementation Details</td><td>130</td></tr><tr><td>6.1.3</td><td>Experimental Study</td><td>131</td></tr><tr><td>6.1.4</td><td>Discussion</td><td>145</td></tr><tr><td>6.2</td><td>Branch-and-Bound Online Selection</td><td>146</td></tr><tr><td>6.2.1</td><td>Methodology</td><td>146</td></tr><tr><td>6.2.2</td><td>Experimental Study</td><td>147</td></tr><tr><td>6.2.3</td><td>Discussion</td><td>160</td></tr><tr><td>6.3</td><td>Conclusion</td><td>160</td></tr><tr><td>7</td><td>Conclusion</td><td>163</td></tr></table>
<table><tbody><tr><td>7.1</td><td>Concluding Remarks</td><td>163</td></tr><tr><td>7.1.1</td><td>Theoretical Models</td><td>163</td></tr><tr><td>7.1.2</td><td>Empirical Models</td><td>164</td></tr><tr><td>7.1.3</td><td>Offline Algorithm Selection</td><td>165</td></tr><tr><td>7.1.4</td><td>Online Algorithm Selection</td><td>166</td></tr><tr><td>7.2</td><td>Future Work</td><td>167</td></tr><tr><td>7.2.1</td><td>Theoretical Models</td><td>167</td></tr><tr><td>7.2.2</td><td>Empirical Models</td><td>168</td></tr><tr><td>7.2.3</td><td>Offline Algorithm Selection</td><td>169</td></tr><tr><td>7.2.4</td><td>Online Algorithm Selection</td><td>170</td></tr><tr><td colspan="2"><strong>Acronyms</strong></td><td><strong>171</strong></td></tr><tr><td colspan="2"><strong>Notation</strong></td><td><strong>173</strong></td></tr><tr><td colspan="2"><strong>List of References</strong></td><td><strong>175</strong></td></tr><tr><td colspan="2"><strong>A Empirical Model Results</strong></td><td><strong>185</strong></td></tr><tr><td>A.1</td><td>PLS – 2 objectives</td><td>186</td></tr><tr><td>A.2</td><td>PLS – 3 objectives</td><td>189</td></tr><tr><td>A.3</td><td>PLS – 5 objectives</td><td>192</td></tr><tr><td>A.4</td><td>BHV-DP – 2 objectives</td><td>195</td></tr><tr><td>A.5</td><td>BHV-DP – 3 objectives</td><td>198</td></tr><tr><td>A.6</td><td>BHV-DP – 5 objectives</td><td>201</td></tr><tr><td>A.7</td><td>GEPS – 2 objectives</td><td>204</td></tr></tbody></table>
<figure><img src="image_6.png" /></figure>
# Chapter 1

## Introduction

### 1.1 Motivations

Multi-Objective Optimization (MOO) deals with the study of problems that have more than one objective to optimize. Such problems may arise in many real-world applications. For example, when considering the problem of finding the optimal path between two cities, several objectives may be considered simultaneously, such as minimizing the traveling time, minimizing the monetary cost of the trip, and maximizing the number of attractions in the path.

One difficulty that commonly arises when tackling MOO problems is the conflicting nature of the objectives. This typically implies that there exists no single solution that optimizes all objectives simultaneously. Instead, there is a set of *efficient solutions*, namely the *efficient set*, that represent the optimal trade-offs between the objectives. A solution is said to be *efficient* if there is no other solution that is better or equal for all objectives, and strictly better for at least one objective. Going back to the previous example, the use of highways instead of national roads typically minimizes the time needed to travel between two cities. However, toll costs increase the overall cost of the trip. On the contrary, opting for national roads can minimize the total cost of the trip, but is likely to increase the travel time. As such, for a non-trivial trip, there may be a large number of efficient solutions that consider different stretches of highway and national roads.

Under the assumption that no preference regarding the objectives is known a priori, every efficient solution is potentially relevant to a Decision Maker (DM). Therefore, a common goal in MOO is to find the entire efficient set. Unfortunately, finding the entire set in a reasonable amount of time can be infeasible for many problems [15]. As such, it is often more useful to find a good approximation to the efficient set, namely an *approximation set*, within a
reasonable amount of time. However, it is not always clear what is a reasonable amount of time or a good approximation set when designing an algorithm. For example, consider an algorithm for a real-time system that can take more or less time to respond to an event depending on the urgency of the situation.

*Anytime algorithms* [11, 78] return an approximation for any time budget, thus offering a trade-off between execution time and approximation quality. As a result, anytime algorithms are particularly suited for situations like the one described above, where the *anytime preferences* of the *DM* with respect to the desired approximation set quality or available time budget are uncertain. Interestingly, many commonly used *MOO* algorithms are inherently anytime since they keep and extend an *archive* of feasible solutions, which can be returned at any time.

Although many anytime *MOO* algorithms have been proposed over the years, it is often not clear which algorithm should be selected to solve a previously unseen problem instance. This selection should take into account the anytime preferences of the *DM* since they may impact the choice. For example, consider a real-time system that receives problem instances to solve. For some instances, the system needs to return a solution within a short amount of time, e.g., 1 second. For other instances, the system has more time, e.g., 60 seconds. For yet other instances, the exact time depends on external factors that are not yet fully established, but it is known that the algorithm will be interrupted between two time points, e.g., between 1 and 60 seconds. Depending on the time available to solve the problem instance, a different algorithm may be chosen.

The problem of automatically selecting an algorithm is commonly known as the *Algorithm Selection Problem (ASP)*, which was introduced by Rice [67] and has gathered the attention of many researchers, on different areas, in recent years [41, 44, 69]. However, research for the *ASP* has mostly focused on selecting between single-objective optimization algorithms, and existing methodologies are not easily extendable to the multi-objective case [41]. Moreover, existing algorithm selection methodologies often consider that the anytime preferences of the *DM* are known in advance, i.e., when designing or training the selection methodology. However, in scenarios like the one described in the previous paragraph, anytime preferences are only known when performing the selection, or during the execution of the selected algorithm itself.

## 1.2 Goals & Scope

In this thesis, we are interested in studying the *ASP* for selecting between *MOO* algorithms to solve a previously unseen problem instance. We consider,
two main perspectives for algorithm selection, offline and online. In the offline case, our goal is to select the best anytime algorithm to solve a previously unseen problem instance, before attempting to solve it. Moreover, we consider that the anytime preferences of the DM are only known when calling the algorithm selection methodology. In the online case, we want to select and change algorithms while solving a previously unseen problem instance. In this case, we consider that the anytime preferences of the DM are unknown, thus the goal is to optimize the quality of the returned approximation set over time, namely to improve anytime performance.

Algorithm selection methodologies often make use of *theoretical* or *empirical models* to predict the performance of algorithms on previously unseen instances. In this thesis, we consider the design of such models for predicting the anytime performance of *MOO* algorithms. We distinguish between theoretical and empirical models in the sense that latter are built using empirically collected data from previous runs of an algorithm on known problem instances, whereas the former rely only on the theoretical knowledge of how an algorithm operates.

We consider that the anytime performance of a MOO algorithm is characterized by a function that maps the time taken by an algorithm, to the quality of the approximation set that would be returned if the algorithm was interrupted at that time. For measuring time, we consider the CPU-time taken by the algorithm. However, since CPU-time depends on the empirical environment, we consider the number of iterations performed by the algorithm when designing theoretical models. For measuring the quality of an approximation set as a scalar value, we consider the *hypervolume indicator* [80], a commonly used quality indicator in MOO that measures the multi-dimensional area dominated by the image of an approximation set in the objective space, with respect to a given reference point.

As such, the main research questions considered for this thesis are the following:

**R1.** Can we design theoretical models to predict the anytime performance of particular MOO algorithms for previously unseen instances with respect to the number of iterations performed by the algorithm, taking into account the theoretical behavior of the algorithm?

**R2.** Can we build empirical models to predict the anytime performance of MOO algorithms for previously unseen instances from a set of training instances?
**R3.** Can we design an offline algorithm selection methodology to select between MOO algorithms for solving a previously unseen instance using the previously designed empirical models of anytime performance?

**R4.** Can we design an online algorithm selection methodology to select and swap between MOO algorithms such that anytime performance is optimized?

**1.3 Contributions**

The main contributions of this thesis, which directly connect to the research questions presented above, are as follows:

**C1.** We propose, and empirically analyze, two variants of theoretical model of anytime performance for bi-objective optimization algorithms that find, at each iteration, a solution maximizing the accumulated hypervolume. Anytime performance is considered in terms of the number of iterations performed by the algorithm, and the hypervolume of the archive of collected solutions. The first variant is given by an analytical methodology that while simple, requires several assumptions. In particular, it assumes that the non-dominated set can be approximated by a piecewise linear function with two symmetric segments and no discontinuities. Furthermore, it requires that the objective space is scaled down such that the non-dominated set is contained in the unit square [0, 1]², and that the reference point for the hypervolume is set to the origin. The second variant is given by an algorithm that generalizes the ideas of the first variant to require less assumptions. In particular, this variant only requires that the non-dominated set can be approximated by a piecewise linear function.

**C2.** We propose three empirical models to predict the anytime performance of a MOO algorithm for a previously unseen problem instance. For building these models we assume that there is a set of training problem instances that can be solved by the algorithm. These models differ in their assumptions, and in their predictive power. In particular, the first model has no other assumptions besides having a set of training instances that can be characterized by a set of instance features, and having collected anytime performance data for the algorithms on those instances. However, it can only predict anytime performance for which anytime performance data has been collected, e.g., if anytime performance data was only collected up to a CPU-time of 10 seconds, it cannot predict anytime
performance for a CPU-time of 20 seconds. On the other hand, the sec-
ond model can predict beyond the collected anytime performance data
used for training, but assumes that the anytime performance of an algo-
rithm can be characterized by a (non-linear) mixed effects model. The
third model, assumes that a theoretical model of anytime performance
that maps the number of iterations to the archive quality exists for the
current algorithm. Its predictive accuracy depends on the quality of this
theoretical model, and on the accuracy of predicting how much compu-
tational time each iteration takes. Lastly, we perform an empirical study
on the first model.

**C3.** We propose, and empirically analyze, an offline algorithm selection framework to select the best anytime MOO algorithm to solve a previously unseen problem instance, under the assumption that the anytime preferences of the decision maker are only (partially) known when performing the selection and can be characterized by a utility function with respect to time and quality. This framework considers empirical models of anytime performance, such as the ones described in contribution C2, to predict the anytime performance of the algorithms on the previously unseen problem instance. Then, a scalar measure of anytime performance is computed for each algorithm, which takes into account the predicted anytime performance and the utility function describing the anytime preferences of the DM. Finally, the algorithm that optimizes this measure of anytime performance is selected.

**C4.** We propose, and empirically analyze, an online algorithm selection methodology to automatically select and swap between MOO algorithms in order to optimize anytime performance. In particular, this methodology selects between different branch-and-bound approaches depending on the active nodes and archive of collected solutions.

The contributions above also led to new anytime MOO approaches:

**C5.** We propose an $\varepsilon$-constraint scalarization approach for bi-objective optimization problems guided by the theoretical model of anytime performance from contribution C1.

**C6.** We propose an indicator-based branch-and-bound framework for MOO that follows a best-first search strategy. This framework considers the use of a binary quality indicator to select the best node to be explored at each iteration. This approach is considered in the context of the online selection methodology of contribution C4.
1.4 Publications & Software

The following publications, in chronological order, are a direct result from the
work carried out for this thesis:

**P1.** A. D. Jesus, L. Paquete, and A. Liefooghe. "A Model of Anytime Algorithm Performance for Biobjective Optimization Problems". In: *Proceedings LeGO - 14th International Global Optimization Workshop*. LeGO 2018. Vol. 2070. AIP Conference Proceedings. AIP, 2019, p. 020049. doi: 10.1063/1.5090016 [39]

In this work, presented at the 14th International Workshop on Global Optimization (LeGO) organized at the Leiden University, we propose and empirically analyze the analytical theoretical model of contribution C1.

**P2.** A. D. Jesus, L. Paquete, and A. Liefooghe. "A Model of Anytime Algorithm Performance for Bi-Objective Optimization". In: *Journal of Global Optimization* 79 (2020), pp. 329-350. doi: 10.1007/s10898-020-00909-9 [38]

In this work, we propose and empirically analyze the algorithmic theo-
retical model of contribution C1, as well as, the ε-constraint approach of
contribution C5.

**P3.** A. D. Jesus, A. Liefooghe, B. Derbel, and L. Paquete. "Algorithm Selection of Anytime Algorithms". In: *Proceedings of the 2020 Genetic and Evolutionary Computation Conference*. GECCO 2020. Association for Computing Machinery, 2020, pp. 850–858. doi: 10.1145/3377930.3390185 [36]

In this work, we propose and empirically analyze the first empirical
model of contribution C2, as well as, the offline algorithm selection frame-
work of contribution C3.

**P4.** A. D. Jesus, L. Paquete, B. Derbel, and A. Liefooghe. "On the Design and Anytime Performance of Indicator-based Branch and Bound for Multi-objective Combinatorial Optimization". In: *Proceedings of the 2021 Genetic and Evolutionary Computation Conference*. GECCO 2021. Association for Computing Machinery, 2021, pp. 234–242. doi: 10.1145/3449639.3459360 [37]

In this work, we propose and empirically analyze the indicator-based
branch-and-bound of contribution C6.
**P5.** A. D. Jesus, L. Paquete, A. Liefooghe, and B. Derbel. "Techniques to Analyze the Anytime Behavior of Algorithms for Multi-Objective Optimization". 31st European Conference on Operational Research (EURO 2021). 2021 [40]

In this talk, we present an R package to measure, and visually analyze, the anytime performance of anytime MOO algorithms.

The following work, is not directly related to this thesis, but analyzes techniques that were used in the implementation of the algorithms for this thesis:

**P6.** D. M. Dias, A. D. Jesus, and L. Paquete. "A Software Library for Archiving Nondominated Points". In: *Proceedings of the 2021 Genetic and Evolutionary Computation Conference Companion*. GECCO 2021. Association for Computing Machinery, 2021, pp. 53–54. DOI: 10.1145/3449726.3462737 [13]

In this work we present, and empirically analyze, a C++ library, available at [14], to perform operations on sets of points in the objective space. Currently, it includes functions to filter the non-dominated points of a set, and to update sets of non-dominated points.

The following software resulted from the development of this thesis and has been made publicly available:

**S1.** A. D. Jesus. *mooutils*. Version v0.1.0. Zenodo, 2022. DOI: 10.5281/zenodo.6855879 [35]

A C++ library of multi-objective optimization utilities, including: functions to compare solutions in terms of dominance and lexicographic orders, data structures to keep mutually non-dominated solutions, data structures and functions for quality indicators such as the hypervolume, different solution queues for search algorithms that keep and process (partial) solutions in specific orders.

**S2.** A. D. Jesus. *mobkp*. Version v0.1.1. Zenodo, 2022. DOI: 10.5281/zenodo.6857821 [33]

A C++ library for the *Multi-Objective Binary Knapsack Problem (MOBKP)* containing several state-of-the-art algorithms, including the algorithms of contributions C5 and C6, as well as, a novel data generation procedure described in contribution C2.

**S3.** A. D. Jesus. *moco_abm*. Version v0.2.0. Zenodo, 2019. DOI: 10.5281/zenodo.3548869 [34]
A Rust implementation of the algorithmic anytime performance model of contribution C1.

S4. A. D. Jesus. *apm*. Version v0.1.1. Zenodo, 2022. doi: 10.5281/zenodo.6857541 [31]

A C++ library containing the theoretical models of anytime performance of contribution C1.

S5. A. D. Jesus. *anytime*. Version v0.0.2. Zenodo, 2022. doi: 10.5281/zenodo.6856120 [30]

An R library for the analysis of anytime performance presented in publication P5.

Lastly, the code for the experiments carried out in this thesis is publicly avail-
able at [32].

**1.5 Outline**

The remainder of this thesis is organized as follows. In Chapter 2, we provide the general concepts and definitions related to MOO, anytime performance and algorithm selection. In Chapter 3, we present the theoretical models and experimental analysis of contribution C1, as well as, the ε-constraint algorithm of contribution C5. In Chapter 4, we present the empirical models and experimental analysis of contribution C2. In Chapter 5, we present the offline algorithm selection framework and experimental analysis of contribution C3. In Chapter 6, we give the online selection methodology and experimental analysis of contribution C4, as well as the branch-and-bound algorithm of contribution C6. In Chapter 7, we summarize and discuss the main findings of this work, and possible directions for future work.
# Chapter 2

## Background

In this chapter, we provide the technical background for this thesis. It is organized as follows. In Section 2.1, we give the concepts and notations related to MOO. In Section 2.2, we present the quality indicators that are used to measure the quality of an approximation set, and discuss their properties. In Section 2.3, we introduce a benchmark MOO problem and present several approaches to solve it. In Section 2.4, we provide the concepts related to anytime algorithms and anytime performance. Lastly, in Section 2.5, we present the ASP and discuss some techniques that have been employed to solve this problem.

### 2.1 Multi-Objective Optimization

A MOO problem is characterized by having several, typically conflicting, objectives to be optimized. Assuming, without loss of generality, that all objectives are to be maximized, a MOO problem with $m$ objectives is defined as:

$$ \underset{x \in \mathcal{X}}{\operatorname{argmax}} f(x) = (f_1(x), f_2(x), \dots, f_m(x)) \qquad (2.1) $$

where $\mathcal{X}$ denotes the set of feasible solutions, and $f_i(x) \to \mathbb{R}$ denotes the $i$-th objective function. The space containing $\mathcal{X}$ is called the *decision space*. If the solutions follow a combinatorial structure we say that the problem is a *Multi-Objective Combinatorial Optimization (MOCO) problem*. The image $f(x)$ of a solution $x \in \mathcal{X}$ is called an *objective vector*. The set of all objective vectors is denoted by $\mathcal{Y} = \{f(x) : x \in \mathcal{X}\}$. The space containing $\mathcal{Y}$ is called the *objective space*.

To compare any two points in the objective space we consider the following *binary dominance relations*.

**Definition 2.1** (Binary dominance relations). Let $y^1, y^2 \in \mathbb{R}^m$. We say that:
* $y^1 \ge y^2$ ($y^1$ weakly dominates $y^2$) iff $y_i^1 \ge y_i^2$ for all $i \in \{1, \dots, m\}$.

* $y^1 > y^2$ ($y^1$ dominates $y^2$) iff $y^1 \ge y^2$ and $y^1 \ne y^2$.

* $y^1 > y^2$ ($y^1$ strictly dominates $y^2$) iff $y_i^1 > y_i^2$ for all $i \in \{1, \dots, m\}$.

* $y^1$ and $y^2$ are mutually non-dominated iff $y^1 \nleq y^2$ and $y^2 \nleq y^1$.

There is often no unique solution that simultaneously optimizes all the objectives in a MOO problem. Instead, there are several efficient solutions that represent the trade-offs between the objectives.

**Definition 2.2 (Efficiency and non-dominance).** A solution $x \in \mathcal{X}$ is efficient, and its objective vector $f(x)$ is non-dominated, iff there is no other solution $x' \in \mathcal{X}$ such that $f(x') > f(x)$ holds. The set of all efficient solutions is called the efficient set and denoted by $\mathcal{X}_E$. The image of the efficient set in the objective space, $\mathcal{Y}_N = \{f(x) : x \in \mathcal{X}_E\}$, is called the non-dominated set.

In MOCO it is often useful to distinguish between supported and non-supported solutions [16].

**Definition 2.3 (Supported solutions).** An efficient solution $x \in \mathcal{X}_E$ is said to be supported if it is also a solution to the weighted sum problem

$$ \underset{x \in \mathcal{X}}{\operatorname{argmax}} \quad \lambda f(x) \qquad (2.2) $$

where $\lambda \in \mathbb{R}_{>0}^m$ denotes a vector of positive weights. Otherwise, $x$ is said to be non-supported. The image in the objective space of a (non-)supported solution is called a (non-)supported objective vector. The sets of all supported solutions and objective vectors are denoted by $\mathcal{X}_S$ and $\mathcal{Y}_S$ respectively.

**Definition 2.4 (Extreme supported solutions).** A solution $x \in \mathcal{X}_E$ is said to be an extreme supported solution if it is a supported solution and its objective vector $f(x)$ is an extreme point of the convex hull of $\mathcal{Y}_N$. The image of an extreme supported solution is denoted an extreme supported objective vector. The sets of all extreme supported solutions and objective vectors are denoted by $\mathcal{X}_{ES}$ and $\mathcal{Y}_{ES}$, respectively.

We also consider the following lexicographical relations.

**Definition 2.5 (Lexicographical relations).** Let $y^1, y^2 \in \mathbb{R}^m$. We say that:

* $y^1 >_{\text{lex}} y^2$ ($y^1$ is lexicographically greater than $y^2$) iff $y_j^1 > y_j^2$ subject to $j = \min\{i : y_i^1 \ne y_i^2\}$.
* $y^1 \ge_{\text{lex}} y^2$ ($y^1$ is lexicographically greater than or equal to $y^2$) iff $y^1 = y^2$ or $y^1 >_{\text{lex}} y^2$.

The problem of finding the lexicographically optimal solution, assuming, without loss of generality, maximizing objectives, is denoted as:

$$ \underset{x \in \mathcal{X}}{\operatorname{arglexmax}} (f_1(x), \dots, f_m(x)) \qquad (2.3) $$

Note that every lexicographically optimal solution is also efficient [16].

In the following, we consider several definitions related to sets of solutions or points in the objective space. First, we consider the definition of a minimal set.

**Definition 2.6 (Minimal set).** Let $X^1, X^2 \subseteq \mathcal{X}$ denote two sets of feasible solutions. We say that $X^1$ is a minimal set of $X^2$ iff $X^1 \subseteq X^2$ and for all $x^2 \in X^2$ there exists $x^1 \in X^1$ such that $f(x^1) = f(x^2)$.

Lastly, to compare between two sets of points in the objective space we define the following binary set dominance relations.

**Definition 2.7 (Binary set dominance relations).** Let $Y^1, Y^2 \subseteq \mathbb{R}^m$ denote two sets of points in the objective space. We say that:

* $Y^1 \ge Y^2$ ($Y^1$ weakly dominates $Y^2$) iff for all $y^2 \in Y^2$ there exists $y^1 \in Y^1$ such that $y^1 \ge y^2$.

* $Y^1 > Y^2$ ($Y^1$ dominates $Y^2$) iff $Y^1 \ge Y^2$ and $Y^2 \not\ge Y^1$.

* $Y^1 > Y^2$ ($Y^1$ strictly dominates $Y^2$) iff for all $y^2 \in Y^2$ there exists $y^1 \in Y^1$ such that $y^1 > y^2$.

* $Y^1$ and $Y^2$ are mutually non-dominated iff $Y^1 \not\ge Y^2$ and $Y^2 \not\ge Y^1$.

## 2.2 Quality Indicators

The need to quantify the quality of a set as a scalar value often arises when comparing sets, e.g., to determine which of two mutually non-dominated sets is better. To this end, quality indicators have been proposed in the literature [43, 81]. We consider the following two classes of quality indicators.

**Definition 2.8 (Unary quality indicator).** Let $A \subset \mathbb{R}^m$ be a set of points in the objective space. A unary quality indicator is a function $I(A) \to \mathbb{R}$ that maps the set to a real value.
**Definition 2.9 (Binary quality indicator).** Let $A, B \subset \mathbb{R}^m$ be two sets of points in the objective space. A binary quality indicator is a function $I(A, B) \to \mathbb{R}$ that assigns a real value to set A with respect to set B.

In the following, we describe one unary and two binary quality indicators and their properties. These indicators were chosen for this thesis since they are commonly used in MOO, and because their order-preserving properties are of particular interest for the analysis of performance of anytime algorithms.

### 2.2.1 Hypervolume Indicator

A frequently used unary quality indicator is the hypervolume indicator [80], or simply hypervolume, which corresponds to the multi-dimensional area dominated by a set of points in the objective space, with respect to a reference point. Formally, for a set $A \subset \mathbb{R}^m$ and a reference point $r \in \mathbb{R}^m$, the hypervolume is defined by:

$$I_H(A) = \Lambda(\{q \in \mathbb{R}^m | \exists a \in A: r \ge q \ge a\}) \quad (2.4)$$

where $\Lambda$ denotes the Lebesgue measure.

One important property of the hypervolume indicator is that it is strictly order-preserving with respect to the binary dominance relations introduced in Section 2.1 under a mild assumption [81]. In particular, let $A, B \subset \mathbb{R}^m$ denote two sets of points in the objective space such that every point in A and B strictly dominates the reference point $r \in \mathbb{R}^m$. Then, $A \ge B$ implies that $I_H(A) \ge I_H(B)$, and $A > B$ implies that $I_H(A) > I_H(B)$.

This property gives a relevant implication for our work. In particular, if an anytime algorithm iteratively adds solutions to an archive set, then the hypervolume for that archive will not decrease. If the objective vector of the solution added to the archive is not weakly dominated by any other point in the archive, and if it strictly dominates the reference point, then the hypervolume for the archive will increase. Moreover, among all subsets of feasible objective vectors, the hypervolume is maximal for any set that contains the non-dominated set.

Another relevant property is that the hypervolume indicator is (strictly) order-preserving with respect to a (strictly) order-preserving scaling of the objectives [43]. This means that we can, for example, normalize the objectives without affecting the order of the approximation sets given by the hypervolume indicator. However, note that the order of the sets given by the hypervolume is affected by the setting of the reference point.
### 2.2.2 Binary Hypervolume Indicator

A binary quality indicator related to the hypervolume is the binary hypervolume indicator [79], or simply binary hypervolume, which corresponds to the multi-dimensional area dominated by one set but not by another, that is, given two sets $A, B \subset \mathbb{R}^m$, the binary hypervolume is given by:

$$I_H(A, B) = I_H(A \cup B) - I_H(B) \qquad (2.5)$$

A relevant property of this indicator for this work, in particular for the indicator-based branch-and-bound approach described in Section 6.1, is that, under some mild assumptions, the binary hypervolume is (strictly) order-preserving when fixing the second parameter.

**Proposition 2.1.** Let $A, B, C \subset \mathbb{R}^m$ denote three sets of points in the objective space such that every point in $A, B$, and $C$, strictly dominates the reference point $r \in \mathbb{R}^m$, and $A$ dominates $C$. Under these assumptions, the binary hypervolume is order-preserving for a fixed reference set $R$, that is:

$$A \geq B \implies I_H(A, C) \geq I_H(B, C) \qquad (2.6)$$

*Proof.* Since $A \geq B$ and $A > C$ then it holds that:

$$A \cup C \geq B \cup C \qquad (2.7)$$

Then, since all points in $A, B$, and $C$, strictly dominate the reference point and because the hypervolume is order-preserving, it holds that:

$$A \cup C \geq B \cup C \implies I_H(A \cup C) \geq I_H(B \cup C) \qquad (2.8)$$

Lastly:

$$I_H(A \cup C) \geq I_H(B \cup C) \implies \qquad (2.9)$$

$$I_H(A \cup C) - I_H(C) \geq I_H(B \cup C) - I_H(C) \implies \qquad (2.10)$$

$$I_H(A, C) \geq I_H(B, C) \qquad (2.11)$$

□

Under the same assumptions, it similarly holds that the relation $A > B$ implies $I_H(A, C) > I_H(B, C)$.

It can similarly be shown that when fixing the first parameter, the binary hypervolume indicator is order-reversing with respect to the dominance of the second parameter, which is relevant to improve the performance of the branch-and-bound approach proposed in Section 6.1.
**Proposition 2.2.** Let A, B, C ⊂ ℝ^m denote three sets of points in the objective space. Then binary hypervolume is order-reversing for a fixed first parameter, that is:

$$B \geq C \implies I_H(A, B) \leq I_H(A, C) \qquad (2.12)$$

*Proof.* By definition:

$$I_H(A, B) \leq I_H(A, C) \implies \qquad (2.13)$$

$$I_H(A \cup B) - I_H(B) \leq I_H(A \cup C) - I_H(C) \implies \qquad (2.14)$$

$$I_H(A \cup B) - I_H(A \cup C) \leq I_H(B) - I_H(C) \qquad (2.15)$$

Since $B \geq C$, this is equal to:

$$I_H(A \cup B \cup C) - I_H(A \cup C) \leq I_H(B \cup C) - I_H(C) \qquad (2.16)$$

which, given the order-preserving property of the unary hypervolume, is equivalent to:

$$I_H(B \cup C) - I_H(C) - \epsilon \leq I_H(B \cup C) - I_H(C) \qquad (2.17)$$

where $\epsilon \geq 0$ corresponds to the intersection of the multi-dimensional area that is covered by both A and B, but not by C. □

Note that by fixing one of the parameters, we could define the indicator as unary. However, we will not be setting a fixed parameter for the whole duration of an algorithm. Instead, we will perform some comparisons with one of the two parameters fixed at specific stages of the algorithm which is why these properties are relevant, but otherwise we keep changing both parameters as the algorithm progresses.

### 2.2.3 The $\epsilon$-indicator

Another commonly used binary quality indicator is the $\epsilon$-indicator [81]. For two sets A, B ⊂ ℝ^m>0 the (multiplicative) $\epsilon$-indicator corresponds to the smallest factor that can be applied to all points in A such that it weakly dominates B, formally:

$$I_{\epsilon}(A, B) = \max_{r \in B} \min_{a \in A} \max_{i \in \{1, \dots, m\}} r_i / a_i \qquad (2.18)$$

where m is the number of objectives.

Like in the previous section, a relevant property for this indicator is that it is order-reversing when fixing the second parameter.

**Proposition 2.3.** Let A, B, C ⊂ ℝ^m>0 denote three sets of points in the objective space. Then, the $\epsilon$-indicator is order-reversing when fixing the second parameter, that is:

$$A \geq B \implies I_{\epsilon}(A, C) \leq I_{\epsilon}(B, C) \qquad (2.19)$$
*Proof.* If $A \ge B$ then it holds that:

$$ \forall b \in B \exists a \in A \forall i \in \{1, \dots, m\} a_i \ge b_i \qquad (2.20) $$

Then, it holds that:

$$ \forall c \in C \forall b \in B \exists a \in A \forall i \in \{1, \dots, m\} c_i/a_i \le c_i/b_i \qquad (2.21) $$

which implies that $I_{\varepsilon}(A, C) \le I_{\varepsilon}(B, C)$. $\square$

<figure><img src="image_7.png" /></figure>

It also holds that the $\varepsilon$-indicator is order-preserving when fixing the first parameter.

**Proposition 2.4.** Let $A, B, C \subset \mathbb{R}_{>0}^m$ denote three sets of points in the objective space. Then, the $\varepsilon$-indicator is order-preserving when fixing the first parameter, that is:

$$ B \ge C \implies I_{\varepsilon}(A, B) \ge I_{\varepsilon}(A, C) \qquad (2.22) $$

*Proof.* If $B \ge C$ then it holds that:

$$ \forall c \in C \exists b \in B \forall i \in \{1, \dots, m\} b_i \ge c_i \qquad (2.23) $$

Then, it holds that:

$$ \forall a \in A \forall c \in C \exists b \in B \forall i \in \{1, \dots, m\} b_i/a_i \ge c_i/a_i \qquad (2.24) $$

which, from the definition of the $\varepsilon$-indicator, implies that $I_{\varepsilon}(A, B) \ge I_{\varepsilon}(A, C)$. $\square$

<figure><img src="image_8.png" /></figure>

## 2.3 Solving the Multi-Objective Binary Knapsack Problem

Throughout this thesis, we will use the MOBKP as a benchmark problem in our experimental studies. This problem has been commonly used as a benchmark problem in many MOCO studies and many algorithms have been proposed to solve it, which makes it an interesting problem for algorithm selection and anytime performance modeling. In the following, we introduce the MOBKP, and discuss several algorithmic approaches that have been proposed to solve it and which we will consider in this thesis.

Given a set of $n$ items, each of which with an associated weight and $m > 1$ values, the MOBKP is the problem of finding the subset of items that maximizes the sum for each of the $m$ values, such that the sum of the weights of the chosen
items does not exceed a given capacity $W \in \mathbb{R}_{>0}$. More formally, it can be
defined as follows:

$$
\begin{equation}
\begin{array}{ll}
\underset{x \in \{0,1\}^n}{\operatorname{argmax}} & f(x) = \left( f_1(x) = \sum_{i=1}^{n} v_i^1 x_i, \dots, f_m(x) = \sum_{i=1}^{n} v_i^m x_i \right) \\
\text{s.t.} & \sum_{i=1}^{n} w_i x_i \le W
\end{array}
\tag{MOBKP}
\end{equation}
$$

where $v_i^j \in \mathbb{R}_{>0}$, $i \in \{1, \dots, n\}$, $j \in \{1, \dots, m\}$, denotes the value for item $i$
and objective $j$, and $w_i \in \mathbb{R}_{>0}^n$ denotes the weight for item $i \in \{1, \dots, n\}$. The
binary variable $x_i$ denotes whether or not an item is chosen for the knapsack,
$x_i = 1$ or $x_i = 0$ respectively.

When solving a MOCO problem the goal is typically to find the efficient
set using an *exact approach*. However, for large instances, finding the efficient
set in a reasonable amount of time may not be feasible. In such cases, the goal
is to find a set of feasible solutions that gives a good approximation in a rea-
sonable amount of time using an *heuristic approach*. In the following sections
we describe several approaches for the MOBKP that are used throughout this
thesis. A library containing implementations for these approaches is available
at [33]. Note that many of these approaches can be adapted to other problems
as well.

**2.3.1 Scalarization Techniques**

Scalarization techniques transform a MOO problem into a series of scalarized
single-objective problems [16].

Dichotomic Weighted Sum

One commonly used scalarization technique for the bi-objective case ($m = 2$)
is the Dichotomic Weighted Sum (DWS) method by Aneja and Nair [1]. This
approach works by performing a dichotomic search on the weight values for a
weighted-sum scalarization given by:

$$
\begin{equation}
\begin{array}{ll}
\underset{x \in \{0,1\}^n}{\operatorname{argmax}} & \lambda_1 f_1(x) + \lambda_2 f_2(x) \\
\text{s.t.} & \displaystyle\sum_{i=1}^{n} w_i x_i \le W
\end{array}
\tag{2.25}
\end{equation}
$$

where $\lambda_1, \lambda_2 \in \mathbb{R}_{>0}$ represent the weight given to each objective. Note that, this
scalarized problem corresponds to a single-objective binary knapsack problem,
which can be solved by using an appropriate algorithm. The DWS method can
be described as follows [1, 19]:
* Step 1: Compute lexicographical optimal solutions $x^1$ and $x^2$ for problems $\text{arglexmax}_{x \in \mathcal{X}} (f_1(x), f_2(x))$ and $\text{arglexmax}_{x \in \mathcal{X}} (f_2(x), f_1(x))$, respectively. Let $Q = \{(1, 2)\}$, $R = \{x^1, x^2\}$, and $k = 3$.

* Step 2: If $Q$ is empty, stop the algorithm and return $R$. Otherwise, choose an element $(r, s) \in Q$ and let $Q = Q \setminus \{(r, s)\}$.

* Step 3: Let $\lambda_1 = f_2(x^s) - f_2(x^r)$, $\lambda_2 = f_1(x^r) - f_1(x^s)$. Let $x^k$ be the optimal solution for the scalarized problem of Equation (2.25) with $\lambda_1$ and $\lambda_2$.

* Step 4: If $f(x^k) \neq f(x^r)$ and $f(x^k) \neq f(x^s)$, then let $R = R \cup \{x^k\}$, $S = S \cup \{(r, k), (k, s)\}$, $k = k + 1$. Go to Step 2.

At the end of the algorithm above the output set $R$ is a subset of $X_S$. If there are multiple optima at step 3 and the algorithm chooses one that maximizes $f_1(x)$ [5], then the output set $R$ is also a minimal set of $X_{ES}$ [1]. Since $R$ cannot be guaranteed to match the efficient set, we consider this to be a heuristic approach. A generalization of this dichotomic scheme for more than two objectives is given by Przybylski et al. [66].

## $\epsilon$-constraint approaches

Another class of commonly used scalarization approaches are $\epsilon$-constraint approaches. These approaches solve a sequence of scalarized problems where all objectives but one are turned into constraints. As an example, a scalarized $\epsilon$-constraint problem for the MOBKP, where all objectives except the first one are turned into constraints, is defined as:

$$ 
\begin{array}{l}
\underset{x \in \{0,1\}^n}{\operatorname{argmax}} f_1(x) \\
\text{s.t.} \quad \sum_{i=1}^{n} w_i x_i \leq W \\
\qquad f_2(x) > \epsilon_2 \\
\qquad \dots \\
\qquad f_m(x) > \epsilon_m
\end{array}
\qquad (2.26) $$

where $\epsilon_j$ denotes the constraint value for objective $j \in \{2, \dots, m\}$. We remark that the image of a solution to this problem can be dominated, but not strictly dominated, by an objective vector in the non-dominated set $Y_N$ [16].

A commonly used $\epsilon$-constraint method for two objectives is given by Srinivasan and Thompson [70]. This method starts by solving a scalarized problem with an arbitrarily small constraint $\epsilon_2$ on the second objective. Then, at each
iteration, given an optimal solution $x^*$ to the scalarized problem in the previous iteration, it sets $\epsilon_2 = f_2(x^*)$ as the constraint value for the new scalarized problem. A generalization of this idea for two or more objectives is the *Adaptive $\epsilon$-constraint (AEPS)* method given by Laumanns et al. [45]. This approach works by partitioning the whole objective space into an hypergrid with respect to the objective functions $f_2$ to $f_m$ based on the efficient solutions already identified. For each grid cell, up to $m$ scalarized problems are solved to find an efficient solution. If an efficient solution is found, then the hypergrid is updated based on the objective values of that solution, otherwise the region for that grid cell is added to a set of already searched regions. The algorithm continues until all grid cells have been explored. We remark that, if the solver for the scalarized $\epsilon$-constraint problems returns all optimal solutions for each problem, then this algorithm can find the efficient set. However, if the algorithm can only return one optimal solution for each problem, as is often the case, then the algorithm finds a minimal set of the efficient set.

### 2.3.2 Branch-and-Bound

Branch-and-Bound (BB) approaches implicitly visit all feasible solutions, which are added to an archive as the algorithm progresses, by recursively dividing the decision space such that each subdivision involves a subproblem. This recursive division is commonly represented by a search tree over the decision variables such that each *node* corresponds to a *partial solution* with some decision variables being fixed, and a *leaf* corresponds to a complete solution. The subproblem of each *node* gives rise to a lower and upper bound set. A *lower bound set* is a non-empty set of feasible solutions to the subproblem. An *upper bound set* is a non-empty set of pairwise mutually non-dominated points in the objective space that weakly dominates the non-dominated set of the subproblem [17, 65].

In the following we present two standard BB frameworks. The first, given in Algorithm 1, corresponds to an *eager* BB framework where the lower and upper bound sets of the nodes are computed as soon as that node is found. The second, given in Algorithm 2, corresponds to a *lazy* BB framework where the lower and upper bound sets are computed when the node is selected for exploration. For both approaches we start by defining a set of active nodes $Q$ initially containing the *root node* $r$ of the search tree. For the MOBKP, this node is the partial solution with no fixed variables. Then, archive $S$, which is used to keep an up to date set of pairwise mutually non-dominated solutions, is initialized. In the case of the *eager* framework the archive is initialized with the lower bound set of the *root node*. For the case of the *lazy* framework the
archive is initially empty.

Then, at each iteration, a node is selected from the set of active nodes $Q$ to be explored. In the eager framework, this exploration consists of computing the branching nodes, and adding them to the set of active nodes $Q$ if the upper bound set of the branching node is not strictly dominated by the archive $S$. Otherwise, the branching node can be discarded. If the branching node is not discarded, then the lower bound set of the branching node is also used to update the archive set. In the lazy framework, we start by evaluating whether or not the upper bound set of the current node is strictly dominated by the archive $S$. If this is true, then that node can be discarded. Otherwise, the lower bound set of the node is used to update the archive, and the branching nodes are added to the set of active nodes $Q$.

<caption>**Algorithm 1:** Eager Branch-and-Bound Framework</caption>

Input : Root node $r$
Output : Solution set $s$
1 $S \leftarrow GetLowerBound(r)$
2 $Q \leftarrow \{r\}$
3 **while** $Q \neq \emptyset$ **do**
4 &nbsp;&nbsp;&nbsp;&nbsp;node $\leftarrow SelectNode(Q)$
5 &nbsp;&nbsp;&nbsp;&nbsp;$Q \leftarrow Q \setminus \{node\}$
6 &nbsp;&nbsp;&nbsp;&nbsp;**for** branch $\in$ GetBranches(node) **do**
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $S \neq$ GetUpperBound branc h **then**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$S \leftarrow S \cup GetLowerBound branc h$)
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$Q \leftarrow Q \cup \{branch\}$)
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
11 &nbsp;&nbsp;&nbsp;&nbsp;**end**
12 **end**
13 **return** $S$

We remark that most approaches in the literature, e.g., [52, 65, 75], follow an eager approach. Moreover, we note that the computation of a lower bound set is not strictly required since the leaf nodes correspond to all complete feasible solutions. However, by using a lower bound set, the algorithm can add more complete solutions to the archive earlier, which can lead to more branches being discarded earlier. This can often outweigh the computational cost of finding the lower bound set and leads to better anytime performance, and execution time.

In the following sections we discuss which strategies can be employed by BB approaches to: i) select the next node to be explored; ii) compute the
<caption>**Algorithm 2:** Lazy Branch-and-Bound Framework</caption>

Input : Root node r
Output : Solution set s
1 S ← ∅
2 Q ← {r}
3 while Q ≠ ∅ do
4   node ← SelectNode(Q)
5   Q ← Q \ {node}
6   if S ≠ GetUpperBound(node) then
7     S ← S ∪ GetLowerBound(node)
8     for branch ∈ GetBranches(node) do
9       Q ← Q ∪ {branch}
10    end
11 end
12 end
13 return S

branching nodes; and iii) compute the lower and upper bound sets.

### Node Selection

Two commonly used search strategies for selecting the next node to be explored are Depth First Search (DFS) [52, 75] and Breadth First Search (BFS) [75], which can be implemented in Algorithms 1 and 2 by using a last-in first-out queue and a first-in first-out queue for the set of active nodes Q, respectively. We remark, that the BFS strategy often suffers from memory issues since it may need to keep a large amount of nodes in memory. On the other hand, the DFS strategy can easily get stuck deep in the tree by attempting to explore nodes that are already dominated by those in the archive.

Therefore, a third strategy, namely Best First Search (BeFS), is to select the most promising node in the set of active nodes Q according to some heuristic. In single-objective optimization this often corresponds to finding the node with the best, e.g., largest in the case of a maximizing objective, upper bound, since it suggests that searching in that direction will lead to a solution with a better objective value. For MOO this is not as straightforward since we are often dealing with incomparable upper bound sets. This issue is further discussed in Chapter 6. A fourth and final strategy, namely Best Depth First Search (BeDFS), is to select the most promising node among those at the deepest level of the search tree.
## Branching Strategies

We consider the *dichotomic branching* strategy, which consists of selecting the next variable that has not yet been fixed, and fixing it to either 0 or 1.

One important aspect is how to define what is the next variable to be fixed. One option is to consider a fixed order for the variables. In particular, we considered an arbitrary order, as well as the $O^{sum}$, $O^{max}$, and $O^{min}$ orders proposed by Bazgan et al. [3]. The latter orders are based on the orders $O^j$ that are induced by the increasing ratios $v_i^j/w_i$ of the items $1 \le i \le n$ for each objective $1 \le j \le m$. In particular, let $r_i^j$ denote the rank of item $i$ for order $O^j$. Then, $O^{sum}$ denotes the increasing order induced by the sums of the ranks of each item for orders $O^j$, that is the increasing order induced by $r_i^1 + \dots + r_i^m$ for each item $i$. The order $O^{max}$ gives the increasing order induced by the maximum rank of each item for orders $O^j$. To break ties, the maximum rank for this order is computed as:

$$r_{\max}(i) = \max_{1 \le j \le m} r_i^j + \frac{1}{mn} \sum_{j=1}^{m} r_i^j \quad (2.27)$$

Lastly, $O^{min}$ denotes the increasing order of the minimum rank of each item for orders $O^j$. To break ties, the minimum rank is computed as:

$$r_{\min}(i) = \min_{1 \le j \le m} r_i^j + \frac{1}{mn} \sum_{j=1}^{m} r_i^j \quad (2.28)$$

Our preliminary experiments revealed that the order $O^{sum}$ very often gave better results than the $O^{min}$ and $O^{max}$ orders. As such, we will not show results for these two orders throughout this thesis. Nonetheless, they are implemented along with the algorithms in [33].

Alternatively, it is possible to consider a dynamic order that varies between nodes. However, this is something we did not explore in this thesis. Cerqueus et al. [6] provides a recent overview of other branching orders and strategies for the MOBKP when $m = 2$, some of which can be generalized for more dimensions.

## Bound Sets Computation

To compute the lower bound set of a node, i.e., a partial solution to the MOBKP, we first consider the extensions given by Dantzig's greedy algorithm [9] for the single-objective problems given by each objective $j$, $1 \le j \le m$. In particular, assume that the items with index 1 to $k$ have already been fixed for the partial solution of the current node and that the items with index $i$, $k < i \le n$, are sorted in non-decreasing order of the ratio $v_i^j/w_i$ for objective $j$. Then, an
<figure><img src="image_9.png" /><figcaption>Figure 2.1: Illustration for two objectives of the image of the lower bound set in the objective space {$ℓ^1, ℓ^2, ℓ^w$} and upper bound set {$u$} for a partial solution with objective vector $y$.</figcaption></figure>

extension for the single-objective problem given by objective $j$ consists of adding all items from $k + 1$ up to, but not including, $b$ such that adding item $b$ would break the weight constraint. The objective vector for such an extension is given by:

$$ℓ^j = \left( \sum_{i=1}^{k} x_i v_i^1 + \sum_{i=k+1}^{b} v_i^1, \dots, \sum_{i=1}^{k} x_i v_i^m + \sum_{i=k+1}^{b} v_i^m \right) \qquad (2.29)$$

In addition, to refine the lower bound set with a solution that offers a more balanced trade-off among the objectives, we consider the extension that is found with Dantzig's greedy algorithm for the items with index $i$, $k < i \le n$, sorted by the non-decreasing ratio $(v_i^1 + \dots + v_i^m) / w_i$. The objective vector of this extension is denoted by $ℓ^{sum}$. We give an illustration of $ℓ^1, ℓ^2$, and $ℓ^{sum}$ for two objectives on the left-hand side of Figure 2.1.

As for the upper bound set, we consider the same set used by Lukata and Teghem [52]. In particular, for the single-objective problem given by each objective $j$, $1 \le j \le m$, value $u_j$ denotes the upper bound by Martello and Toth [54]. Then, a point $u$ is given by:

$$u = (u_1, \dots, u_j) \qquad (2.30)$$

which characterizes an upper bound set {$u$}. This upper bound set is illustrated on the right-hand side of Figure 2.1.

### 2.3.3 Dynamic Programming

Another class of exact approaches is Dynamic Programming (DP). Several DP approaches have been proposed for the MOBKP, many of which extend the
algorithm by Nemhauser and Ullmann [58] for more objectives. Klamroth and Wiecek [42] extended the recurrence formulations of previous DP approaches into a more general framework for the MOBKP. More recently, Bazgan et al. [3] presented a DP approach based on several complementary dominance relations. These dominance relations allow to discard partial solutions that cannot lead to efficient solutions. Figueira et al. [19] proposed some extensions to this algorithm for the bi-objective ($m = 2$) case, by considering a global lower bound set and different techniques to discard partial solutions based on their upper bound. Delort and Spanjaard [12] also considered the use of bound sets to discard partial solutions that cannot lead to efficient solutions for the bi-objective case.

In this work, we consider the Bazgan-Hugot-Vanderpooten Dynamic Programming (BHV-DP) algorithm by Bazgan et al. [3]. This algorithm consists of $n$ stages, such that at any stage $k$ the algorithm generates a set of feasible solutions $X_k$ where the first $k$ items have been set to either 0 or 1, and the remaining items are all set to 0. At each stage, the states are generated according to the following recursion:

$$X_k = \text{Dom}(X'_k) \qquad (2.31)$$

$$X'_k = X_{k-1} \cup \left\{ (\dots, x_{k-1}, 1, x_{k+1}, \dots) : \sum_{i=1}^{n} w_i x_i \le W, x \in X_{k-1} \right\} \qquad (2.32)$$

where the initial state $X_0 = \{(0, \dots, 0)\}$ contains the solution with no chosen items for the knapsack and $\text{Dom}(\cdot)$ denotes one or more dominance relations that can be used to discard solutions that cannot lead to efficient solutions at later stages. The following dominance relations are proposed in [3] for discarding solutions:

* **Dom1:** Let $x, x' \in X'_k$. Then, $x$ can be discarded if $x'$ derives from $x$ at stage $k$, i.e., if $x' = (\dots, x_{k-1}, 1, x_{k+1}, \dots)$, and if $\sum_{i=1}^{n} w_i x_i \le W - \sum_{i=k}^{n} w_i$.

* **Dom2:** Let $x, x' \in X'_k$. Then, $x$ can be discarded if $f(x') > f(x)$ and $k = n$, or if $f(x') > f(x)$, $k < n$ and $\sum_{i=1}^{n} w_i x'_i \le \sum_{i=1}^{n} w_i x_i$.

* **Dom3:** Let $x, x' \in X'_k$. Then, $x$ can be discarded if its upper bound set, $\text{UB}(x)$, is strictly dominated by the image of the lower bound set of $x'$, i.e., $\{f(x'') : x'' \in \text{LB}(x')\} > \text{UB}(x)$.

Note that, dominance relation **Dom2** is an extension for the multi-objective case of the dominance relation by Nemhauser and Ullmann [58]. For dominance relation **Dom3**, the authors consider a lower bound set with two solutions that are found by Dantzig's greedy algorithm [9] with respect to the orders $O^{\text{max}}$
and $O^{sum}$. The upper bound set is computed using the upper bound extensions for each objective function by Martello and Toth [54], which has been discussed in the previous section.

For the bi-objective case, we also consider the *B-DP1* extension by Figueira et al. [19], denoted by Figueira-Paquete-Simões-Vanderpooten Dynamic Programming (FPSV-DP). This approach further considers the following dominance relation:

* **Dom4:** Let $x \in X'_{k}$. Then, x can be discarded if $Y_{ES} > UB(x)$, where $Y_{ES}$ denotes the set of extreme supported solutions that can be computed using the DWS approach.

A possible extension for more dimensions would be to consider the algorithm by Przybylski et al. [66] to compute the $Y_{ES}$ for more dimensions, or to consider an alternative set of feasible solutions.

### 2.3.4 Pareto Local Search

Local search approaches keep an archive of solutions that are iteratively explored. At each iteration, a solution that has not yet been explored is selected from the archive, and its neighboring solutions are considered for inclusion in the archive. A neighborhood is defined according to one or more neighborhood mappings of the form $N: X \to P(X)$, where $P(X)$ is the power set of X, that map a solution to all its neighbors. A local search approach stops once the neighborhood of every solution in the archive has been explored. A Pareto Local Search (PLS) approach [60] considers that the archive contains only mutually non-dominated solutions. Therefore, a neighboring solution is only added to the archive if it is not dominated by any solution currently in the archive. Moreover, when a solution is added to the archive every solution in the archive that becomes dominated by the newly inserted solution is removed.

Algorithm 3 describes a PLS approach for the MOBKP following the 1-flip-exchange neighborhood described in Liefooghe et al. [48]. This algorithm takes as input a set $X_0$ of mutually non-dominated solutions. This set can, for example, be generated randomly, generated by a heuristic algorithm, or simply contain a single solution corresponding to the empty knapsack. Sets $X_u$ and $X_a$ denote the set of unexplored mutually non-dominated solutions and the archive of mutually non-dominated solutions, respectively. The Select method selects a solution from the set of unexplored solutions $X_u$. Different selection mechanisms may be considered. For this thesis, we consider that one solution is selected randomly. The 1-flip-exchange neighborhood considers two neighborhood mappings. The first mapping, denoted by FlipNeighborhood in
the algorithm, gives all solutions that result from exchanging any item in the current solution. The second mapping, denoted by ExchangeNeighborhood in the algorithm, gives all solutions that result from exchanging two items in the current solution. Following the algorithm description in [48], the algorithm only explores the second mapping if no new non-dominated solution was found for the first. It is worth noting that, for this neighborhood exploration definition, the PLS approach may not be able to find the efficient set [48].

<caption>**Algorithm 3: Pareto Local Search**</caption>

**Input:** $X_0$ (Set of initial mutually non-dominated solutions)
1 &nbsp;&nbsp;&nbsp;&nbsp;$X_u \leftarrow X_0$
2 &nbsp;&nbsp;&nbsp;&nbsp;$X_a \leftarrow X_0$
3 &nbsp;&nbsp;&nbsp;&nbsp;**while** $X_u \neq \emptyset$ **do**
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*x* $\leftarrow$ Select($X_u$)
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$X_u \leftarrow X_u \setminus \{x\}$
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;aux $\leftarrow$ 0
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** *x'* $\in$ FlipNeighborhood(*x*) $\cup$ ExchangeNeighborhood(*x*) **do**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** *x'* $\in$ ExchangeNeighborhood(*x*) $\land$ aux = 1 **then**
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;break
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** IsDominated(*x'*, $X_a$) **then**
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;continue
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*$X_u$* $\leftarrow$ RemoveDominated($X_u$, *x'*)
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*$X_a$* $\leftarrow$ RemoveDominated($X_a$, *x'*)
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*$X_u$* $\leftarrow$ $X_u \cup \{x'\}$
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*$X_a$* $\leftarrow$ $X_a \cup \{x'\}$
18 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** *x'* $\in$ FlipNeighborhood(*x*) **then**
19 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;aux $\leftarrow$ 1
20 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
21 &nbsp;&nbsp;&nbsp;&nbsp;**end**
22 &nbsp;&nbsp;&nbsp;&nbsp;**end**
23 &nbsp;&nbsp;&nbsp;&nbsp;**return** $X_a$

## 2.4 Anytime Algorithms

Anytime algorithms [11, 78], in the context of MOO, can return a set of feasible solutions for any time budget. This means that a DM can trade-off execution
time with solution quality. This is relevant in various scenarios, for example,
in real-time systems where the available time for the algorithm varies between
calls, or in composite systems where the allocated time for each algorithm
varies between instances in order to improve the result of the overall system.

A distinction is often made between two categories of anytime algorithms,
*interruptible* and *contract* [78]. Interruptible algorithms do not require a priori
knowledge of the available time budget and can return a set of feasible solutions
when interrupted during their execution. By contrast, contract algorithms
require that the time budget is known a priori, and may or may not return a
solution if interrupted before that time ends. *MOO* approaches that keep an
archive of feasible solutions, such as the ones described in the previous section,
fall into the interruptible category since the archive can always be returned
when the algorithm is interrupted. As such, in this thesis we will focus on
interruptible algorithms.

To study the performance of a single run of an anytime algorithm we define a *performance trace* that describes the trade-off between solution quality and execution time. In the following, a run refers to a single execution of an algorithm on a particular instance.

**Definition 2.10 (Performance trace).** The *performance trace* of a run *r* of an algorithm *a* is given by a function *Q*<sub>*a*,r</sub> : R<sub>≥0</sub> → R that maps execution time to the quality of the solution set that would be returned if the algorithm was interrupted at that time.

A performance trace is said to be monotonic if function $Q_{a,r}(t)$ is non-
decreasing for increasing values of $t$. If every run of an anytime algorithm gives
a monotonic performance trace then that algorithm is said to have *monotonic
behavior*. Otherwise, the algorithm is said to have *non-monotonic behavior*.
Note that algorithms with non-monotonic behavior can typically be made
to have monotonic behavior by keeping the archive with the best quality in
memory. However, this may not always be possible or desirable. For example,
if the objective functions are computationally expensive, the algorithm may
use approximate objective functions that are not order-preserving with respect
to the true objective functions, which may lead to exchanging better solutions
with worse ones. As such, in the following, we want to consider definitions
for anytime performance that are meaningful for both monotonic and non-
monotonic behaviors.

To characterize the anytime performance of an algorithm for a particular
instance, we consider the definition of a *performance profile*.

**Definition 2.11 (Performance Profile).** The *performance profile* of an algorithm
<figure><img src="image_10.png" /><figcaption>Figure 2.2: Performance traces (left) and profile (right) of three runs of an algorithm.</figcaption></figure>

a on a problem instance $\iota$ is given by a function:

$$P_{a,\iota}(t, q) \rightarrow [0, 1] \qquad (2.33)$$

that denotes the probability of a run of algorithm a on problem instance $\iota$ finding a solution with quality greater than or equal to $q$ at time $t$.

We also consider the definition of an *empirical performance profile* which summarizes the anytime performance of an algorithm over a set of runs.

**Definition 2.12** (Empirical Performance Profile). For an algorithm a and a set of runs $R$, an *empirical performance profile* is given by a function:

$$\hat{P}_{a,R}(t, q) = \frac{1}{|R|} \sum_{r \in R} \begin{cases} 1 & \text{if } Q_{a,r}(t) \ge q \\ 0 & \text{otherwise} \end{cases} \qquad (2.34)$$

that denotes the proportion of runs in $R$ that were able to achieve a solution with quality greater than or equal to $q$ at time $t$.

The definition of similar (empirical) performance profiles as a conditional probability function has been previously considered by Hoos and Stützle [29], denoted by *runtime distribution*, and by Chiarandini [7] through the use of the *empirical attainment function* [25]. Those profiles consider the probability of finding a solution of quality greater than or equal to $q$ at, or before, execution time $t$. This definition is equivalent to ours if the performance traces of an algorithm are all monotonic, but is different when considering algorithms with non-monotonic performance traces.

Figure 2.2 illustrates the performance traces of three runs of an algorithm on the left, and the corresponding empirical performance profile over those
runs on the right. We remark that this gives a similar plot as the empirical
attainment function [50].

Lastly, we remark that aggregating performance traces into a single profile
should be carefully done [29, 78]. For example, if the considered instances have
different quality domains, then their aggregation may not be meaningful. One
possibility to mitigate this issue is, for example, to consider a relative measure
of quality.

**2.5 Algorithm Selection**

The ASP [67] concerns the selection of an algorithm to solve a previously
unseen instance, such that a given performance measure is optimized. As
discussed in Chapter 1, we distinguish between offline algorithm selection, i.e.,
the problem of selecting an algorithm before attempting to solve the given
problem instance, and online algorithm selection, i.e., the problem of selecting
and swapping algorithms while solving the given problem instance.

Offline Algorithm Selection

Most existing algorithm selection methodologies fall into the category of offline
algorithm selection. Moreover, existing methodologies often consider selecting
between single-objective optimization algorithms rather than MOO algorithms,
and do not consider the anytime preferences of the DM at the time selection.
We refer to Kotthoff [44] and Kerschke et al. [41] for two recent reviews of
existing offline algorithm selection methodologies.

These methodologies are generally divided into two categories, *regression* and *classification*. Regression approaches predict the performance measure for each algorithm and select the algorithm with maximal value, see for example, Leyton-Brown et al. [46] and Xu et al. [77]. Classification approaches instead select an algorithm without predicting the performance measure, for example, by learning which algorithm gives the best performance measure for different kinds of instances and utility functions, see for example, Vilas Boas et al. [74] and Polyakovskiy et al. [64].

Note that, most offline algorithm selection methodologies make use of *in-
stance features* to characterize a problem instance as a set of measurable val-
ues. For single-objective optimization problems many, general-purpose and
problem-specific, instances features have been proposed, e.g., Nudelman et al.
[59] and Mersmann et al. [55, 56]. For MOO, there has been less research. How-
ever, there are some recent promising developments with respect to general-
purpose instance features, e.g., Daolio et al. [10] and Liefooghe et al. [47, 49].
## Online Algorithm Selection

Online algorithm selection methodologies are less common than offline, but existing methodologies similarly focus on selecting between single-objective optimization algorithms. Some methodologies monitor the performance of the algorithms to detect which ones are more promising and should be given more time to run, e.g., Gagliolo et al. [23] and Gagliolo and Schmidhuber [22]. Other approaches, consider selecting between algorithms for exploring a search tree, e.g., Arbelaez et al. [2] selects the best search strategy at checkpoints in the search tree. Up to our knowledge, no online algorithm selection methodology has focused on selecting between (anytime) MOO algorithms. We refer to Kotthoff [44] for a review of existing online algorithm selection methodologies.

## Related Topics

Lastly, two related areas of research are meta-level control and composition of anytime algorithms [78]. Meta-level control, in the context of anytime algorithms, relates to the problem of deciding when to interrupt an anytime algorithm and act on the returned solution [28, 71]. Composition of anytime algorithms, denotes the problem of building a complex system comprised of one or more anytime algorithm that maximizes the overall performance of the system [68]. Up to our knowledge, these problems have not been considered in the context of MOO problems.

## 2.6 Conclusion

In this chapter, we presented the theoretical background for multi-objective optimization, and discussed several algorithms for the MOBKP, a benchmark problem that will be used throughout this thesis.

We also discussed the concepts related to anytime algorithms and anytime performance since most MOO algorithms naturally fall into this category of algorithms. Another important aspect to note is that, for non-trivial problem instances, MOO exact approaches cannot commonly find the efficient set in a feasible amount of time, and even heuristic algorithms may take too long to complete in some cases. Moreover, preferences with respect to time budget and target quality may not be known when developing the algorithms or training the model, but rather only when selecting the best algorithm to solve a particular problem instance. As such, it is often more relevant to look at MOO algorithms as anytime algorithms, and consequently to model and study their anytime performance rather than simply looking at their performance at the
end of execution or at a fixed time budget or target quality. This motivates us to look at the anytime performance of the MOO algorithms when tackling the ASP.

Lastly, we provided some context regarding the ASP. In particular, we high-
lighted that the ASP has been studied for the selection of single-objective al-
gorithms, often considering fixed time budgets or target qualities, over several
decades. However, researchers have only recently started looking at the ASP
in the context of MOO algorithms. Moreover, these recent works have only
considered offline selection with time budgets that were set before training the
selection model. Still, these recent developments, and the positive results that
are known for the selection of single-objective algorithms, show that this is a
relevant problem to address.
# Chapter 3

## Theoretical Model of Anytime Performance

Theoretical models can be used to predict the performance of an algorithm on previously unseen instances without the need to empirically evaluate the algorithm on a set of training instances beforehand. In this chapter, we give a theoretical model to approximate the ideal anytime performance of bi-objective optimization algorithms that (attempt to) collect an efficient solution at each iteration [38, 39]. Examples are the scalarization techniques introduced in Section 2.3.1. Anytime performance is considered with respect to the number of iterations and the quality of the archive of collected solutions measured with the hypervolume. We also propose an $\varepsilon$-constraint approach guided by this theoretical model [38] to define the scalarized subproblems that are to be solved at each iteration. The experimental results show that the theoretical model can be used to predict the ideal anytime performance, as well as the anytime performance of the proposed $\varepsilon$-constraint approach.

This chapter is organized as follows. In Section 3.1, we present the theoretical model, and carry out an experimental study to analyze its quality against the ideal anytime performance. In Section 3.2, we give the $\varepsilon$-constraint approach guided by the theoretical model, and carry out an experimental study to analyze its anytime performance. In Section 3.3, we conclude with a summary of the main findings in this chapter and discuss possible directions for future research.

### 3.1 Theoretical Model of Anytime Performance

Consider an iterated algorithm to solve a MOO problem that finds an (efficient) solution to the problem at each iteration, and keeps an archive of all mutually
non-dominated solutions found. We describe this algorithm by a function A that returns the set of solutions found after $i > 0$ iterations:

$$A(i) = A(i - 1) \cup x^i \qquad (3.1)$$

where $A(0) = \emptyset$, and $x^i \in \mathcal{X}$ denotes the efficient solution found at iteration $i$. The performance trace of this algorithm with respect to the number of iterations and the hypervolume of the image of the archive in the objective space, is given by:

$$Q(i) = I_H(\{f(x) : x \in A(i)\}) \qquad (3.2)$$

where $f(x)$ denotes the multi-objective function being optimized.

Under the assumption that the number of iterations of the algorithm is unknown, we define an ideal algorithm $A^*$ with complete knowledge of the efficient set $\mathcal{X}_E$ that finds, at each iteration, the solution that contributes the most to the hypervolume of the archive, that is:

$$A^*(i) = A^*(i - 1) \cup \underset{x \in \mathcal{X}_E}{\operatorname{argmax}} I_H(\{f(x') : x' \in A^*(i - 1) \cup \{x\}\}) \qquad (3.3)$$

We assume, for the sake of simplicity, that there exists a single optimum for the problem in the right-hand side of the equation above. Nonetheless, if there are multiple optima, the formulation can be extended to consider the list of all sets that maximize the hypervolume at each iteration. The performance traces of this ideal algorithm gives the ideal performance trace, which can be described by the following function:

$$Q^*(i) = I_H(\{f(x) : x \in A^*(i)\}) \qquad (3.4)$$

Note that this definition of an ideal anytime performance trace does not necessarily give the maximal hypervolume that can be achieved with $i$ efficient solutions. However, a definition that gives, for each iteration, the maximal hypervolume is often incompatible with an algorithm that collects, at each iteration, a single efficient solution, since the sets of solutions for two consecutive iterations can differ by more than one solution.

Finding the ideal performance trace requires knowing the efficient set, which is unrealistic when solving a previously unseen problem instance. Instead, we define a theoretical model to approximate the ideal performance trace for bi-objective optimization problems that does not require such knowledge. This model works in two phases. In the first phase, we define a piecewise linear approximation to the non-dominated set under some assumptions. In particular, we assume that the objective values of the lexicographic optimal solutions are known and that the non-dominated set can be well approximated
by the positive quadrant of a superellipse. Although it is not expected that the non-dominated set matches the positive quadrant of a superellipse exactly, our findings suggest that this gives a good approximation in practice for many problems with linear sum objective functions. Emmerich and Deutz [18] have also considered the use of superellipses for the generation of MOO test problems.

In the second phase, we give two methodologies to compute the performance trace of an oracle that collects, at each iteration, a point in the piecewise linear approximation that maximizes the hypervolume. This oracle is characterized by the increase in hypervolume at each iteration, denoted by $C(i) \to \mathbb{R}$, which corresponds to the hypervolume contribution of the newly found point to the archive. As a result, the performance trace for this oracle, which approximates the ideal performance trace, is given by:

$$ \tilde{Q}^*(i) = \sum_{i=1}^{n} C(i) \qquad (3.5) $$

The first methodology that we give to compute the values of $C(i)$, is an analytical formulation that requires several assumptions to hold. In particular, we assume that the piecewise linear approximation is convex and that it is given by two symmetric linear segments. Moreover, we assume that the reference point for the hypervolume indicator is given by the *nadir point* of the non-dominated set, i.e., the point given by the minimum coordinate values for points in the non-dominated set. The second methodology is an algorithm to compute $C(i)$ for any piecewise linear approximation, and for any reference point setting. Despite being more general, this algorithmic approach is computationally expensive. However, empirical results show that the algorithmic methodology is still very fast in practice, so this is unlikely to be a problem for most applications.

The remainder of this section is organized as follows. In Section 3.1.1, we define the piecewise linear approximation to the non-dominated set. In Sections 3.1.2 and 3.1.3, we present the analytical and algorithmic methodologies respectively. In Section 3.1.4, we present an empirical study comparing the performance trace given by our theoretical model against the ideal performance trace.

### 3.1.1 Estimating the non-dominated set

Let us assume that the non-dominated set for a bi-objective optimization problem with maximizing objectives, scaled down to the unit square $[0, 1]^2$, can be well approximated by the positive quadrant of a superellipse centered in the origin with both semi-diameters of length one. Formally, this is given by the
<figure><img src="image_11.png" /><figcaption>Figure 3.1: Example of the superellipse curve for two distinct values of the curvature parameter $d \in \{0.5, 2\}$.</figcaption></figure>

following parametric equation:

$$y_1^d + y_2^d = 1 \qquad (3.6)$$

where $d > 0$ is a parameter that controls the curvature of the superellipse and $y_1, y_2 \in [0, 1]$. Alternatively, the values of $y_1$ and $y_2$ can also be defined with respect to a parameter $\theta \in [0, \pi/2]$ as follows

$$y_1 = \cos^{2/d} \theta \qquad (3.7)$$

$$y_2 = \sin^{2/d} \theta \qquad (3.8)$$

Note that, for $d < 1$ the resulting approximation is non-convex, whereas for $d \ge 1$ it is convex. This is illustrated in Figure 3.1 where we show the approximation for two distinct values of $d$.

For our theoretical model, we consider an approximation to the positive quadrant of the superellipse defined by a piecewise linear approximation $G$ with $\ell > 0$ linear segments defined between consecutive points in the curve approximation, such that each point $g_i, i \in \{1, \dots, \ell + 1\}$ can be defined by:

$$g_i = \left( \cos^{2/d} \theta_i, \sin^{2/d} \theta_i \right) \qquad (3.9)$$

where $\theta_i < \theta_{i+1}$ and the values of $\theta_i$ are evenly spaced in the interval $[0, \pi/2]$. We found that this simple setting of $\theta_i$ provides a good approximation for the purposes of our model. However, other techniques are available to sample these points, e.g., see [62].

In Figures 3.2 and 3.3, we plot the superellipse curve given by Equation (3.6) for two distinct curvature parameters $d \in \{0.5, 2\}$, as well as the corresponding piecewise linear approximations $G$ for a varying number of segments $\ell$.
<figure><img src="image_12.png" /><figcaption>Figure 3.2: Example of the superellipse curve (continuous line) for $d = 0.5$ and the corresponding piecewise linear approximation (dashed line) for a varying number of segments $l \in \{1, 2, 3, 4\}$.</figcaption></figure>

As expected, the piecewise linear approximation better approximates the superellipse curve as the number of segments increases. Moreover, for $d < 1$ the piecewise approximation weakly dominates the superellipse curve and will therefore overestimate the area dominated by the curve, whereas for $d > 1$ it will underestimate it.

### 3.1.2 Analytical methodology

Under some assumptions, we can define a closed formula for the performance trace of an oracle that collects, at each iteration, a point in the piecewise linear approximation $G$ that maximizes the hypervolume. In particular, we assume
<figure><img src="image_13.png" /><figcaption>Figure 3.3: Example of the superellipse curve (continuous line) for $d = 2$ and the corresponding piecewise linear approximation (dashed line) for a varying number of segments $\ell \in \{1, 2, 3, 4\}$.</figcaption></figure>
that the piecewise linear approximation is convex, i.e., $d \ge 1$, and that it has two linear segments such that $\theta_1 = 0$, $\theta_2 = \pi/4$, and $\theta_3 = \pi/2$. Moreover, the reference point for the hypervolume indicator is given by the *nadir point* of the piecewise linear approximation, i.e., $r = (0, 0)$. Under these assumptions, the piecewise linear approximation can be defined as:

$$G = \{(y_1, g(y_1)) \mid y_1 \in [0, 1]\} \qquad (3.10)$$

$$g(y_1) = \begin{cases} \frac{p-1}{p}y_1 + 1, & 0 \le y_1 \le p \\ \frac{p}{p-1}y_1 + \frac{p}{1-p}, & p < y_1 \le 1 \end{cases} \qquad (3.11)$$

where

$$p = \cos^{2/d} \frac{\pi}{4} = \sin^{2/d} \frac{\pi}{4} \qquad (3.12)$$

The first point collected by the oracle for such a piecewise linear approximation is found by solving the problem:

$$\underset{y_1 \in [0,1]}{\operatorname{argmax}} y_1 \cdot g(y_1) \qquad (3.13)$$

This problem has a single global optimum at coordinates $(p, p)$. As such, the first point collected by the oracle gives a contribution $C(1) = p^2$.

Subsequent points can be found by considering the regions that are not dominated by any of the points previously collected. These regions are said to be *uncovered*. Regions dominated by at least one of the previously collected points are said to be *dominated*. After the first point $(p, p)$ has been collected, there are two uncovered regions each defined by a right triangle such that the hypotenuse corresponds to the non-dominated set of that uncovered region. In Figure 3.4, we illustrate the uncovered and dominated regions, colored in gray and black respectively, after collecting one point (left-hand side) and three points (right-hand side).

To find the largest hypervolume given by the dominated area of a point in an uncovered region defined by a right triangle such that the set of non-dominated points is given by the hypotenuse and the reference point is on the right angle, we solve the following problem:

$$\begin{array}{ll} \max & y_1 \cdot y_2 \\ \text{s.t.} & y_2 = -\frac{c_2}{c_1} \cdot y_1 + c_2 \\ & y_1 \in [0, c_1] \\ & y_2 \in [0, c_2] \end{array} \qquad (3.14)$$

where $c_1$ and $c_2$ denote the length of the catheti of the right triangle on the $y_1$ and $y_2$ axes, respectively. This problem has a global optimum of $c_1 \cdot c_2/4$
<figure><img src="image_14.png" /><figcaption>Figure 3.4: Example of the dominated (in black) and uncovered (in gray) regions after collecting one (left) and three (right) points with the oracle.</figcaption></figure>

at point ($y_1 = c_1/2, y_2 = c_2/2$). Note that, the hypervolume given by the dominated area of this point corresponds to half the area of the uncovered region. Moreover, after collecting this point and considering its dominated area, there are two smaller uncovered regions that are also right triangles, each of which with catheti that are half the length of those in the original triangle, and consequently covering a quarter of the area.

After collecting the first point, the resulting uncovered regions are two equivalent right triangles with catheti of size $(1-p)$ and $p$, and an area of $(1-p) \cdot p/2$. Then, from Equation (3.14), it follows that the second and third contributions for points collected by our oracle are $C(2) = C(3) = (1-p) \cdot p/4$. After excluding the regions dominated by the two points providing these contributions, one for each uncovered region, there are four equivalent uncovered regions as illustrated on the right-hand side of Figure 3.4. After collecting the points that maximize the hypervolume in all equivalent uncovered regions, the number of uncovered regions doubles. As such, a general equation for $C(i)$ is given by

$$C(i) = \begin{cases} p^2 & i = 1 \\ \frac{(1-p) \cdot p}{4[\log_2 i]} & i > 1 \end{cases} \qquad (3.15)$$

which shows a logarithmic rate of convergence. Computing the performance trace of the oracle, given by $\tilde{Q}^*(j)$ until iteration $j > 0$, which is the result of the sum of $C(i)$, $0 < i \le j$, can be done in $O(j)$ if the value of $[\log_2 i]$ is kept updated in constant time.
<figure><img src="image_15.png" /><figcaption>Figure 3.5: On the left, an uncovered region defined by a simple polygon (gray area). On the middle and on the right, the individual segments (black line) of the polygonal chain, and their respective right triangles (gray area).</figcaption></figure>

### 3.1.3 Algorithmic methodology

For the general case, i.e., for any number of linear segments $l > 0$, any curvature parameter $d > 0$, and any reference point $r \in \mathbb{R}^2$, the previous analytical model does not hold. Instead, we present an algorithm to compute the contributions $C(i)$, for iterations $i > 0$, which works by keeping an updated set of the uncovered regions. Note that, only one point from a single uncovered region is collected at each iteration, and that the uncovered regions share no area. As such, to avoid repeated calculations, the algorithm keeps a cache of the points in each uncovered region that contribute the most to the hypervolume of the archive. Moreover, to efficiently find the largest contribution, the uncovered regions can be kept in a priority queue with respect to the largest contribution given by a point in that region.

For the purposes of our algorithm, an uncovered region that results from the piecewise linear approximation can be described by a simple polygon with a right angle on the point with the smallest coordinates, i.e., the reference point for that uncovered region, such that the non-dominated set is given by the simple polygonal chain opposite to that point. As a result, each segment of the polygonal chain is part of the hypotenuse of a right triangle with its right angle on the reference point. In Figure 3.5 we give an illustration of an uncovered region with a polygonal chain made up of two segments, and the corresponding right triangles for each segment.

We know from Equation (3.14) that the point with the largest hypervolume for a right triangle with the reference point on the right angle is given by $z = (c_1/2, c_2/2)$, where $c_1, c_2$ are the lengths of the catheti on axes $y_1, y_2$, respectively. Moreover, the result for that optimization problem is strictly monotonically increasing for $0 \le y_1 \le c_1/2$, and strictly monotonically decreasing for $c_1/2 \le y_1 \le c_1$. As such, the point in a segment of the polygon chain with
<figure><img src="image_16.png" /><figcaption>Figure 3.6: First three steps of the algorithmic methodology for a curvature parameter $d = 0.5$, a number of linear segments $\ell = 2$, and a reference point $r = (0, 0)$.</figcaption></figure>

the largest hypervolume is the point in that segment closest to $z$. Finally, we can iterate over all segments of the polygonal chain to find the point with the largest hypervolume for the whole uncovered region. In Figure 3.6, we show the first three steps of the algorithm for a non-convex piecewise approximation with $\ell = 2$ linear segments.

Note that there may be multiple points that give the same hypervolume at each iteration. This is the case for the first step illustrated in Figure 3.6 where the points (0.5, 1/6) and (1/6, 0.5) both give the same hypervolume for a reference point (0, 0). If there are only two points, and the new uncovered regions that result from selecting either point are equivalent, then we can select either point since the progression of $C(i)$ will be the same regardless of the choice. We say that two uncovered regions are equivalent, if they have the same, possibly mirrored, shape. This is the case for the first step illustrated in Figure 3.6, where selecting the first point leads to uncovered regions that are mirrored by those that result from selecting the second point. By contrast, if the resulting uncovered regions are not equivalent, then we need to consider both choices until a point is found at a later iteration that results in a greater contribution $C(i)$ for either of those choices. We can do this by keeping a list of priority queues to be processed at each iteration, such that each priority queue denotes an alternative choice. The number of priority queues in this list may grow exponentially since we have at most two choices at each iteration. However, in our experimental analysis we never observed more than one or two priority queues in the list at any given time, which suggests that an exponential scenario is unlikely to occur.

If more than two points give the same contribution, then we show that we only need to consider choosing between the two lexicographically extreme points. Let $Z \subseteq G$ be a list of at least three distinct lexicographically sorted
points that give the same contribution $C(i)$ for an iteration $i > 0$. Let $Z^{\min}$ and $Z^{\max}$ be the lexicographically smallest and largest points respectively. Since the points in $Z$ are mutually non-dominated and lexicographically sorted, then $Z^{\min}$ has the smallest first coordinate, and $Z^{\max}$ the smallest second coordinate. Also, let $r \in \mathbb{R}^2$ be the reference point for the hypervolume such that it is strictly dominated by every point in the polygonal chain. Then, for any point $y$ in the polygonal chain $G$ that can be selected in the succeeding iteration and has a larger first coordinate than the previously selected point, it holds that:

$$ \forall z \in Z \setminus \{Z^{\min}\} \ y_1 > z_1 \implies (y_1 - Z_1^{\min})(y_2 - r_2) > (y_1 - z_1)(y_2 - r_2) \quad (3.16) $$

since $Z_1^{\min} < z_1$. This implies that, if the succeeding point has a larger first coordinate than the previously selected point, then choosing $Z^{\min}$ in the previous iteration would yield a strictly larger contribution for $y$. Likewise, if point $y$ has a smaller first coordinate than the previously selected point, it holds that:

$$ \forall z \in Z \setminus \{Z^{\max}\} \ y_1 < z_1 \implies (y_1 - r_1)(y_2 - Z_2^{\max}) > (y_1 - r_1)(y_2 - z_2) \quad (3.17) $$

since $Z_2^{\max} < z_2$. This implies, that if the succeeding point has a smaller first coordinate than the previously selected point, then choosing $Z^{\max}$ in the previous iteration would yield a strictly larger contribution for $y$. As we have covered every possible succeeding point $y \in G$, and either of the two lexicographically extreme points always give a strictly better hypervolume for $y$, then we can discard the other points.

To illustrate this, consider Figure 3.7 where we show the uncovered regions after selecting one of the three points $Z^1 = (0.08, 0.5)$, $Z^2 = (0.2, 0.2)$, $Z^3 = (0.5, 0.08)$ that give the same hypervolume for a reference point $(0, 0)$. We observe that any subsequent point returned in the uncovered regions A and B resulting from choosing point $Z_2$, will have a strictly larger contribution in the uncovered regions A' and B' that result from choosing the lexicographically largest or smallest points $Z_3$ or $Z_1$ respectively.

In the following we give the complexity of our algorithm for computing $C(i)$. At each iteration, the algorithm needs to find an uncovered region with a point that gives the largest contribution, remove it from the priority queue containing all uncovered regions, split the region into two new uncovered regions and insert those regions into the data structure. As such, for the $i$-th call, there are at most $i$ uncovered regions in the data structure that need to be considered. Since the uncovered regions share no area, we only need to compute the points with the largest contribution in a given region once. Then, we can keep the uncovered regions in a priority queue supported by max-heap data structure. Removing the uncovered region with largest hypervolume
<figure><img src="image_17.png" /><figcaption>Figure 3.7: Uncovered regions that result from choosing one of three points with the same hypervolume.</figcaption></figure>

contribution from such a data structure can be done in $O(\log i)$ time, and insert new uncovered regions in $O(1)$ or $O(\log i)$ time depending on the max-heap data structure chosen.

To find the points in each uncovered region that give the largest hypervolume we need to loop over all segments and find the point for each segment with maximal hypervolume. This can be done in constant time for each segment, as discussed previously, and there are at most $l$ segments in each uncovered region, so this takes $O(l)$ time for each region.

Finally, there may be at most $2^i$ priority queues to be processed at each iteration if multiple points yield the same contribution. As a result our algorithm has a worst-case time complexity given by $O(2^i (i\ell + i \log i))$ for function $C(i)$. To approximate $P^*(j)$ we can compute $C(j)$, and sum all the intermediate values $C(i)$, $1 \le i \le j$, so the worst-case time complexity is also given by $O(2^j (j\ell + j \log j))$. Alternatively, we consider a greedy variant that always selects the lexicographically smallest point when multiple points have the same hypervolume. The worst-case time complexity of this greedy variant is given by $O(j\ell + j \log j)$.

### 3.1.4 Empirical study

A C++ implementation of the exact and greedy algorithmic methodologies is available at [31]. In Figure 3.8, we report the average execution times of this implementation for a varying number of iterations, a varying number of linear segments, and curvature parameter $d \in \{0.5, 2\}$. The reference point is set to $r = (0, 0)$ since it had a negligible impact on the results. The execution times were taken with the benchmark library [4] (v.1.6.0) considering a MinTime of 10 seconds, which broadly means that for each combination of parameters the methodologies were repeatedly executed until CPU time was greater than 10 seconds. The experiments were carried out in a computer with an AMD
<figure><img src="image_18.png" /><figcaption>Figure 3.8: Execution times in seconds of the algorithmic methodology for a varying number of iterations, varying number of linear segments, reference point $r = (0, 0)$, and curvature parameter $d \in \{0.5, 2\}$.</figcaption></figure>

Ryzen CPU 5 1600 3.2Ghz processor. The benchmark code and full results are available at [32].

The results show that both variants are very fast (less than 0.1 seconds on average) for the given parameters. We expect most practical scenarios to have a number of linear segments and iterations within the values reported. There is a small difference between the exact and greedy methodologies since the exact methodology has to take into account the possibility of more than one point giving the same maximal contribution. Moreover, if there are multiple points, which in fact happens for the first iteration of $d = 0.5$, the exact methodology attempts to detect if the uncovered regions left after selecting each of those points are equivalent. In this case, the uncovered regions are equivalent. As such, only a single point needs to be considered and the impact is relatively small. Otherwise, the impact on performance could be greater. However, as previously mentioned, in our experiments we did not observe such cases.

In the following, we analyze the quality of our theoretical model on non-dominated sets that arise from particular MOO problems. First, we consider non-dominated sets arising from the Unconstrained Binary Knapsack Problem (UBKP) scaled down to the unit square $[0, 1]^2$, which is formally given by:

$$ \underset{x \in \{0,1\}^n}{\operatorname{argmax}} f(x) = \left( f_1(x) = \frac{\sum_{i=1}^{n} v_i x_i}{\sum_{i=1}^{n} v_i}, f_2(x) = 1 - \frac{\sum_{i=1}^{n} w_i x_i}{\sum_{i=1}^{n} w_i} \right) \quad (3.18) $$

where $x_i$ is a binary variable denoting whether or not item $i$ has been chosen for the knapsack, $v_i \in [0, 1]$ denotes a value for that item, and $w_i \in [0, 1]$ de-
notes a weight for the item. The goal, is to maximize the sum of the values and minimize the sum of the weights of the chosen items. The values and weights of the items are randomly generated according to a multivariate uniform distribution in the range [0, 1], with correlation $\rho \in [-1, 1]$ following the procedure described in [72]. The non-dominated set can be found by the DP algorithm of Nemhauser and Ullmann [58].

Figures 3.9a to 3.9c shows the non-dominated set for instances with correlation $\rho \in \{-0.8, 0.0, 0.8\}$ and problem size $n = 100$. We see that for greater values of $\rho$ the shape of the non-dominated set almost resembles a line, and that for smaller values of $\rho$ the non-dominated set resembles a more accentuated convex curve. The same figures show that the non-dominated set can be well approximated with the positive quadrant of a superellipse with curvature parameter $d = \log_{\rho} 0.5$, where $p$ is found by considering the solution to the following optimization problem:

$$ \underset{x \in \{0,1\}^n}{\operatorname{argmax}} \min \{f_1(x), f_2(x)\} \qquad (3.19) $$

such that if $x'$ is the optimal solution to this formulation, then the parameter is given by $p = (f_1(x') + f_2(x'))/2$.

We remark that the non-dominated sets generated for the UBKP are always approximated by a convex curve. As such, to validate our model for non-convex curves, we consider a variant of this problem characterized by a tight capacity constraint on the number of items to be included in knapsack. We refer to this variant as the Capacity Constrained UBKP (CCUBKP), which, scaled down to the unit square $[0, 1]^2$, is formally defined as:

$$ \underset{x \in \{0,1\}^n}{\operatorname{argmax}} f(x) = \left( f_1(x) = \frac{\sum_{i=1}^n v_i x_i}{\max_{j \in \{1, \dots, n\}} v_j}, f_2(x) = 1 - \frac{\sum_{i=1}^n w_i x_i}{\max_{j \in \{1, \dots, n\}} w_j} \right) \quad (3.20) $$

$$ \text{s.t.} \quad \sum_{i=1}^{n} x_i \le 1 \qquad (3.21) $$

This problem is trivial to solve since the set of feasible solutions contains the solutions that result from adding a single item to the knapsack, or none. However, it allows us to devise instances such that the non-dominated set is non-convex. In particular, instances for this problem are generated by sampling the value vector $v$ and an auxiliary weight vector $w'$ from equally parameterized Weibull distributions for a fixed scale parameter $\lambda = 1$, and a varying shape parameter $s$. Then, the value vector is sorted in increasing order and the auxiliary weight vector is sorted in decreasing order. The weight vector values are then given by $w_i = (\max_{j \in \{1, \dots, n\}} w'_j) + (\min_{j \in \{1, \dots, n\}} w'_j) - w'_i$, $i \in \{1, \dots, n\}$. The non-dominated set for the CCUBKP is trivially found by enumerating all objective vectors.
Note that, by varying the shape parameter s, we vary the curvature of the non-dominated set. In particular, for smaller values of s the curvature is more accentuated, whereas for larger values of s the curvature is less accentuated. In Figures 3.9d to 3.9f we plot the non-dominated sets that arise in problem instances with shape parameter s ∈ {1.6, 2.0, 3.0} and problem size n = 2000. A larger problem size for the CCUBKP was chosen in order to have a similar number of non-dominated points between the two problems. The plots show that when s increases the non-dominated sets generated for the CCUBKP can be well approximated with the positive quadrant of a superellipse. The curvature parameter for the superellipse approximation is given by $d = \log_p 0.5$, where p is found by considering the solution to the following optimization problem:

$$ \underset{x \in \{0,1\}^n}{\operatorname{argmax}} \min \{f_1(x), f_2(x)\} \qquad (3.22) $$

$$ \text{s.t.} \quad \sum_{i=1}^{n} x_i \le 1 \qquad (3.23) $$

such that if $x'$ is the optimal solution to this formulation, then parameter p is given by $p = (f_1(x') + f_2(x'))/2$.

In Figures 3.10 and 3.11 we compare the performance trace given by our theoretical model against the ideal performance trace for the instances described above. We consider a varying number of linear segments for our theoretical model, and two distinct settings for the reference point. The hypervolume of the collected points is given relative to the hypervolume of the non-dominated set. We observe that the theoretical model gives a very good approximation of anytime performance for $l = 10$ and $l = 100$ linear segments, with the latter showing a slightly better prediction for non-dominated sets that can be approximated by a non-convex curve. We also see that for $l = 2$ linear segments, the theoretical model can only adequately estimate the ideal performance trace when the non-dominated more closely resembles a line, and the reference point is further away from the non-dominated set. Lastly, note that the model generally underestimates the hypervolume for the convex case, and overestimates it for the non-convex case.
<figure><img src="image_19.png" /><figcaption>(a) $\rho = -0.8$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(d) $s = 1.6$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) $\rho = 0.0$</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(e) $s = 2.0$</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(c) $\rho = 0.8$</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(f) $s = 3.0$</figcaption></figure>

Figure 3.9: Non-dominated set (points in gray) and respective superellipse approximation (straight black line) for instances of the UBKP with varying correlation $\rho$ and problem size $n = 100$ (left), and for instances of the CCUBKP with varying shape parameter $s$ and problem size $n = 2000$ (right).
<figure><img src="image_25.png" /><figcaption>Figure 3.10: Results for the theoretical model on non-dominated sets of the UBKP with $n = 100$, reference point $r \in \{(0,0), (-1,-1)\}$, and correlation $\rho \in \{-0.8, 0.0, 0.8\}$.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>(a) $s = 1.6, r = (0, 0)$</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(d) $s = 1.6, r = (-1, -1)$</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(b) $s = 2.0, r = (0, 0)$</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(e) $s = 2.0, r = (-1, -1)$</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(c) $s = 3.0, r = (0, 0)$</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(f) $s = 3.0, r = (-1, -1)$</figcaption></figure>

<caption>Figure 3.11: Results for the theoretical model on non-dominated sets of the **CCUBKP** with $n = 2000$, reference point $r \in \{(0,0), (-1,-1)\}$, and shape parameter $s \in \{1.6, 2.0, 3.0\}$.</caption>
## 3.2 Guided ε-constraint

In this section, we present an ε-constraint approach guided by our theoretical model, namely the Guided ε-constraint (GEPS) method. As discussed in Section 2.3.1 an ε-constraint approach solves a sequence of constrained single-objective problems to find efficient solutions.

To find an efficient solution, at iteration $i > 0$, to a bi-objective optimization problem with maximizing objective functions, our approach first solves the following single-objective problem:

$$ \begin{array}{ll} \underset{x \in \mathcal{X}}{\operatorname{argmax}} & f_1(x) \\ \text{s.t.} & f_2(x) > \varepsilon^i \end{array} \qquad (3.24) $$

where $\varepsilon^i$ denotes a constraint on the second objective for iteration $i$. Note that, the problem in Equation (3.24) may yield a weakly efficient solution [16], in particular there may be a solution that gives an equal value for first objective but a better value the second. As such, to guarantee that we find an efficient solution, we consider a second problem at each iteration, which is given by:

$$ \begin{array}{ll} \underset{x \in \mathcal{X}}{\operatorname{argmax}} & f_2(x) \\ \text{s.t.} & f_1(x) \ge f_1(x') \end{array} \qquad (3.25) $$

where $x' \in \mathcal{X}$ denotes the optimal solution found for the problem of Equation (3.24).

To set the $\varepsilon^i$ constraint at iteration $i > 0$, we start by considering the second coordinate of the $i$-th point returned by our theoretical model. However, if it is known that no new non-dominated point can be found for that constraint, then, we consider the first subsequent point collected by our theoretical model for which a new point may be found. In particular, let $y_2^i$ denote the second coordinate of the non-dominated point used at iteration $i$ to set a constraint $\varepsilon^i$. Then, we know that no new non-dominated point can be found for a constraint between $y_2^i$ and $\varepsilon^i$.

In the following, we report the performance trace for this algorithm against the ideal anytime performance. We use the GNU Linear Programming Kit MILP Solver [53] to solve the single-objective subproblems given by Equations (3.24) and (3.25) for both the UBKP and the CCUBKP. For the theoretical model used to guide the algorithm, we consider the greedy algorithmic methodology with a varying number of linear segments $\ell \in \{2, 10, 100\}$ since it returns at most one point at each iteration. Note that, if the assumptions for the analytical methodology are met, the greedy algorithmic methodology gives the same
result as the analytical methodology. The performance trace is given with
respect to the number of iterations and archive quality expressed in terms of
relative hypervolume, i.e., the ratio between the hypervolume of the archive
and the hypervolume of the non-dominated set. We also report archive quality
in terms of relative hypervolume deviation, that is, the different between the
maximal relative hypervolume, i.e., 1, and the actual relative hypervolume.

In Figures 3.12 and 3.13 we show the results for the first 128 iterations
on instances for the **UBKP** with problem size *n* = 100, correlation values
*ρ* ∈ {-0.8, 0.0, 0.8}, and reference point *r* ∈ {(0, 0), (-1, -1)}. The correspond-
ing non-dominated sets are shown in Figures 3.9a and 3.9c. In Figures 3.14
and 3.15 we show the results on instances for the **CCUBKP** for a problem size
*n* = 2000, varying shape parameter *s* ∈ {1.6, 2.0, 3.0}, and reference point
*r* ∈ {(0, 0), (-1, -1)}. The corresponding non-dominated sets are shown in
Figures 3.9d and 3.9f. Looking at the results, we can see that all three settings
for the number of linear segments are quite good, and the **GEPS** achieves an
anytime performance trace very close to the ideal one. Nonetheless, using 10
or 100 linear segments seems to provide better results, which is to be expected
since the non-dominated set can be better approximated by the piecewise ap-
proximation in those cases. We observed that these results were coherent for
instances with different problem sizes, different parameters, and different ref-
erence points. Moreover, we expect the results to generalize for any problem
where the scaled non-dominated set can be approximated by the quadrant of
a superellipse, such as the **MOBKP**.
<figure><img src="image_32.png" /><figcaption>(a) $\rho = -0.8$, $r = (0, 0)$</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(b) $\rho = 0.0$, $r = (0, 0)$</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(c) $\rho = 0.8$, $r = (0, 0)$</figcaption></figure>

<caption>Figure 3.12: Results for the UBKP with $n = 100$, reference point $r = (0, 0)$, and varying correlation $\rho \in \{-0.8, 0.0, 0.8\}$.</caption>
<figure><img src="image_35.png" /><figcaption>(a) $\rho = -0.8, r = (-1, -1)$</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(b) $\rho = 0.0, r = (-1, -1)$</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(c) $\rho = 0.8, r = (-1, -1)$</figcaption></figure>

<caption>Figure 3.13: Results for the UBKP with $n = 100$, reference point $r = (-1, -1)$, and varying correlation $\rho \in \{-0.8, 0.0, 0.8\}$.</caption>
<figure><img src="image_38.png" /><figcaption>Figure 3.14: Results for the CCUBKP with $n = 2000$, capacity constraint $c = 1$, reference point $r = (0, 0)$, and varying shape parameter $s \in \{1.6, 2.0, 3.0\}$.</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>(a) $s = 1.6, r = (-1, -1)$</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>(b) $s = 2.0, r = (-1, -1)$</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>(c) $s = 3.0, r = (-1, -1)$</figcaption></figure>

<caption>Figure 3.15: Results for the CCUBKP with $n = 2000$, capacity constraint $c = 1$, reference point $r = (-1, -1)$, and varying shape parameter $s \in \{1.6, 2.0, 3.0\}$.</caption>
## 3.3 Discussion

In this chapter, we presented a theoretical model to characterize the performance trace given by the trade-off between the number of iterations and the hypervolume for bi-objective optimization algorithms that find, at each iteration, a solution to the problem, such as scalarization techniques. In the experimental studies carried out, we observed that our theoretical model approximates quite well the ideal anytime performance. We also presented an $\varepsilon$-constraint approach guided by our theoretical model, namely the GEPS, which showed a performance trace very close to the ideal one.

As evidenced by this algorithm, a possible use of our model is to design scalarization techniques with good anytime performance. Moreover, our model can be used to monitor and predict the anytime performance of our algorithm or techniques that similarly attempt to maximize the hypervolume at each iteration [61]. For example, it should be possible to detect if an algorithm that has an anytime performance similar to our model is taking too long to find each solution, and whether or not it will achieve some desirable quality within a certain time budget, which may justify a restart or a switch to a different strategy.

One possible direction for future research, is to extend the model for approximations of the non-dominated set other than a quadrant of a superellipse. In practice, this can already be accomplished by our algorithmic methodology by considering any set of mutually non-dominated, possibly disconnected, linear segments. However, the main challenge is to define a good piecewise linear approximation to the non-dominated set. A related idea is to consider a feedback loop between our model and a scalarization technique that sequentially finds efficient solutions. In particular, if the set of linear segments corresponds to an approximation of the non-dominated set, and throughout the execution of the scalarization technique we find the actual non-dominated points, then we can insert these points back into the model to increase the quality of the approximation in an online fashion, which could be used to better predict where the next points in the non-dominated set will be. Lastly, a clear direction for future research is to extend our model for more than 2 objectives. This seems more challenging in part due to the splitting and tracking of the uncovered regions after a point with maximal hypervolume contribution is found, since the regions may no longer be independent, i.e., finding a point in one region may affect the other regions.

Lastly, for the purposes of algorithm selection in this thesis, it would be more useful to have a model that considers anytime performance in terms of
CPU-time rather than the number of iterations. In the next chapter, we discuss an empirical model to convert the number of iterations given by this theoretical model to CPU-time by considering empirically collected data.
# Chapter 4

## Empirical Models of Anytime Performance

In this chapter, we consider the development of empirical models for predicting the anytime performance of an algorithm, which take into account empirically collected anytime performance data from previous runs of the algorithm. Unlike the theoretical model discussed in the previous chapter, these empirical models can be used to predict anytime performance regardless of time or quality unit, as long as the anytime performance data to train the models is collected in terms of the desired units of time and quality. In particular, this is desirable for the algorithm selection methodologies that we will be considering in the following chapters, as we want to consider time in terms of CPU-time, which was not possible with the theoretical model of the previous chapter.

This chapter is organized as follows. In Section 4.1 we propose three abstract frameworks for the development of empirical models that take into account anytime performance traces collected for an algorithm on a set of training instances. Then, in Section 4.2 we propose an empirical model for predicting the anytime performance of algorithms to the MOBKP that follows the first of these proposed frameworks. We also carry out an experimental study considering three algorithms and instances with 2, 3, and 5 objectives, which shows that the model can often predict the anytime performance correctly. However, we also highlight some cases for which the model is not as accurate as we would hope, and we discuss possible reasons for this. Lastly, in Section 4.3 we summarize the main results of this chapter and discuss possible directions for future work.
# 4.1 Frameworks for Empirical Models

In this section we propose three abstract frameworks to guide the development of empirical models of anytime performance. We recall that we want to develop empirical models with the goal of predicting the anytime performance of an algorithm on a previously unseen problem instance. In particular, let $\iota$ denote a previously unseen problem instance, and let $a$ denote an algorithm to solve $\iota$. We consider that an empirical model to predict the performance of algorithm $a$ on instance $\iota$ can be characterized in terms of a performance profile:

$$ \tilde{P}_{a,\iota}(t, q) \rightarrow [0, 1] \qquad (4.1) $$

that denotes the (approximate) likelihood of algorithm $a$ achieving an approximation with quality greater than or equal to $q$ at time $t$, when solving problem instance $\iota$.

In all the presented frameworks, we assume that there is a set of training instances $\mathcal{I}$ that is known, and that there is a set of runs for algorithm $a$ associated to this set of instances, denoted by $R_{a,\mathcal{I}}$. We generally assume that the instances on the training data are not too different from the instances for which we will be predicting anytime performance. In particular, we assume that when predicting the anytime performance for an instance $\iota$, that there is a set of instances in the training data set that give an anytime performance similar to that of instance $\iota$. This is not to say that we cannot predict anytime performance if there are no such instances in the training data set. However, we expect the quality of the prediction to be better if there are.

Another relevant aspect regarding the training data set is that in the presented frameworks the training set may be updated at any time. This may be particularly relevant in real-world scenarios where one possibility to build the training data set is to continuously update the training data set with anytime traces from previously encountered instances.

Lastly, these frameworks are discussed in general terms, since the implementation details for the model can depend on the prediction scenario, e.g., how much human-interaction is possible and what information can be gathered from the problem instances. That being said, an empirical model to predict the anytime performance of algorithms to the MOBKP problem, following the first framework, is presented in Section 4.2.3.

## 4.1.1 Framework I

Given the set of training instances $\mathcal{I}$ and a previously unseen instance $\iota \notin \mathcal{I}$, a model following this framework starts by selecting a subset of training in-
stances, denoted by $\mathcal{I}_S \subseteq \mathcal{I}$ that are expected to have similar anytime performance to instance $\iota$. Then, given a set of runs $R_{a,\mathcal{I}_S}$ denoting the runs of algorithm $a$ on the selected training instances $\mathcal{I}_S$, the model builds and returns an approximated performance profile given by the empirical performance profile for those runs, that is:

$$ \tilde{P}_{a,\iota}^{\mathcal{I}}(t, q) = \hat{P}_{R_{a,\mathcal{I}_S}}(t, q) \qquad (4.2) $$

The key aspect for this model is how to select a subset of training instances $\mathcal{I}_S$ that are expected to have similar anytime performance to the previously unseen instance $\iota$. One possibility is to manually design a set of rules to perform this selection, based on a set of instance features that characterize a problem instance. This requires that the behavior of algorithm $a$ is well known, and that it is obvious how the different instance features affect anytime performance, which, arguably, is not often the case.

Another possibility, that requires less human interaction, is to manually analyze the anytime performance of an algorithm on the training instances, and to group instances that give similar anytime performance. Then, given a set of instance features that characterize an instance, we can use classification techniques to automatically learn to classify the group of instances to which the previously unseen instance $\iota$ belongs. Alternatively, given a set of measurable features that characterize the anytime performance of an algorithm on an instance, we could also build the groups without any human interaction by using an unsupervised partitioning technique that takes into account those features, such as $k$-means clustering.

Lastly, another possibility is to consider automatic learning techniques that can learn to select similar instances based on a set of instance features without the need for human interaction. For example, given a set of instance features that characterize a problem instance, we can use clustering techniques, such as a $k$-nearest neighbor algorithm, that can select instances that have similar instance features, and which we expect to have similar anytime performance. This is the possibility that we consider later in Section 4.2.3 when implementing our model.

## 4.1.2 Framework II

To predict the approximation quality at time $t$ with the previous framework, we assume that the runs of an algorithm on the training instance set were interrupted at a time greater than or equal to $t$, or that the algorithm finished before time $t$, since otherwise there is no empirical data that can be used to build the predicted performance profile. This can lead to infeasible training
times if we want to predict the approximation quality for large values of t. To address this, our second framework considers the use of parametric models that characterize the anytime performance of an algorithm on the training instances I, to predict the anytime performance on a previously unseen instance $i \notin I$.

In the field of statistics, data collected over time for the same subject, in this case the performance trace of a run, is commonly known as longitudinal data [20]. An issue posed by longitudinal data is that measurements taken from the same subject are not independent. To address this, mixed-effects models where each subject is modeled by a separate curve whose parameters are a random perturbation of a baseline curve can be considered [20]. Note that, since the performance traces of an anytime algorithm are not typically expected to be linear, a non-linear mixed-effects model is likely to be more appropriate. Assuming that the mean response of a subject i at time $t_{ij}$ can be modeled in terms of a non-linear regression function and random error, a non-linear mixed effects model is given by:

$$y_{ij} = f(t_{ij}, \beta_i) + \epsilon_{ij} \quad (4.3)$$

where the evolution of the mean response $y_{ij}$ is modeled by a non-linear function $f(t_{ij}, \beta_i)$ and a random error term $\epsilon_{ij}$. The non-linear function $f(\beta_i, t_{ij})$ takes the time parameter $t_{ij}$ and a set of subject-specific parameters $\beta_i$. These parameters $\beta_i$ are commonly defined to linearly depend on a set of covariates $A_i$, and random effects $b_i$, that is:

$$\beta_i = A_i\beta + b_i \quad (4.4)$$

Alternatively, parameters $\beta_i$ can also be modeled by a non-linear function. We refer to [20] and [63] for a more complete description and discussion of non-linear mixed effects models.

In the context of anytime algorithms, Gagliolo et al. [21] considered non-linear mixed-effects models following an exponential decay growth curve to model the anytime performance of a stochastic single-objective optimization algorithm and predict the optimum objective value, given a set of 400 instances and 25 runs of the algorithm for each instance. There were two experiments that defined subjects for the non-linear mixed effects models differently. In the first, the authors fit a separate model for each instance, taking into account the 25 runs gathered for each instance as the subjects in the model. In the second experiment, the authors fit 25 models, such that for each model a single run for each instance is chosen, and the model takes into account the instances as the subjects.

For the development of an empirical model under this second framework, we consider two possibilities to develop the parametric model. For the first
possibility, we consider the idea of the first experiment by Gagliolo et al. [21], and fit a separate model for each training instance based on the available runs. Since each separate model is for a single instance, instance specific features are not mapped to the response by the model itself. As such, we need to learn to predict the parameters $\beta_i$ of the non-linear function that were fit by the model from a set of instance features, e.g., by using supervised learning techniques from machine learning. Then, given a previously unseen instance $\iota \notin \mathcal{I}$, we map its instance features to a set of predicted parameters $\tilde{\beta}_i$ to be considered in the non-linear function. The estimated performance trace for an algorithm a on a previously unseen instance $\iota$ is given by:

$$ \tilde{P}_{a,\iota}^{\Pi}(t, q) = \begin{cases} 1 & \text{if } f_a(t, \tilde{\beta}_i) \ge q \\ 0 & \text{otherwise} \end{cases} \qquad (4.5) $$

where $f_a(t, \tilde{\beta}_i)$ is the non-linear function considered for algorithm a that takes as parameters the time $t$ and the predicted set of parameters $\tilde{\beta}_i$. If, instead, multiple non-linear mixed-effects models were fit for each instance, e.g., a set of non-linear mixed-effects models fit using quantile regression methods [24], then we could learn to predict the parameters from each model and build the estimated performance profile from multiple curves.

The second possibility we consider is to fit a non-linear mixed-effects model that takes into account the instances as subjects, e.g., by considering the median performance trace for each instance. In this case, instance features can be given as parameters to the non-linear function. Then, when given a previously unseen instance $\iota \notin \mathcal{I}$ we can apply the covariates $A_i$ that were fit by the model to the instance features of $\iota$ to get the predicted median anytime performance trace for $\iota$. As before, we can consider different models for different quantiles to have a more complete anytime performance profile. We expect that it is harder to build a model following this second possibility, since a function needs to be defined to map the instance features to the non-linear function parameters, which may not be trivial in many cases. However, it may allow to better predict instances that are quite different from the training instances.

### 4.1.3 Framework III

This framework assumes a theoretical model that characterizes the performance trace of an algorithm on a previously unseen instance with respect to the number of iterations and archive quality, such as the one presented in Section 3.1. In particular, let this theoretical model be denoted by a function $\tilde{Q}_{a,\iota}(i) \to \mathbb{R}$ that returns the predicted quality of the archive at iteration $i \in \mathbb{Z}_{>0}$
for an algorithm a on instance $\iota$. Moreover, let $\mathcal{I}_S \subseteq \mathcal{I}$ denote a set of selected training instances that, compared to instance $\iota$, take a similar amount of CPU-time going from iteration $i$ to $i + 1$ for every $i \in \mathbb{Z}_{\ge 0}$. Lastly, let $R_{a,\mathcal{I}_S}$ denote the runs of algorithm a on the training instances $\mathcal{I}_S$, and $T_{a,r}(i) \to \mathbb{R}_{\ge 0}$ denote a function that returns the time algorithm a took to reach iteration $i$ on run $r$.

Then, we define a function:

$$ \tilde{Q}_{a,\iota,r}(t) = \tilde{Q}_{a,\iota} \left( \min\{i : T_{a,r}(i) \ge t\} \right) \quad (4.6) $$

that denotes an approximate performance trace for instance $\iota$ assuming the CPU-times obtained for run $r$. Finally, the approximate performance profile obtained by this framework is given by:

$$ \tilde{P}_{a,\iota}^{\text{III}}(t, q) = \frac{1}{|R_{a,\mathcal{I}_S}|} \sum_{r \in R_{a,\mathcal{I}_S}} \begin{cases} 1 & \text{if } \tilde{Q}_{a,\iota,r}(t) \ge q \\ 0 & \text{otherwise} \end{cases} \quad (4.7) $$

The key aspect for this framework is how to select training instances that take a similar amount of time on each iteration, compared to a previously unseen instance $\iota$. As in Framework I, we can either consider a set of rules designed by hand, or use (semi)-automated learning techniques to identify similar instances. However, one important difference is that for this model we need to take into account what happens at each iteration. For example, consider the GEPS approach presented in Section 3.2, which can be modeled by the theoretical model presented in Chapter 3. For an iteration to take a similar amount of time when solving two distinct instances, the scalarized $\epsilon$-constraint problems to be solved should have similar difficulty so that they take a similar amount of time to solve. As such, an approach following this empirical model should look at instance features that affect the performance of the solver to the scalarized problem, e.g., how tight the $\epsilon^i$ constraint is at each iteration.

## 4.2 Experimental Study

In this section, we carry out an experimental study for an empirical model implemented according to Framework I. In particular, we consider the use of this model in predicting the anytime performance of PLS, BHV-DP, and GEPS algorithms on the MOBKP considering different numbers of objectives. This section is organized as follows. In Section 4.2.1, we describe the procedure used to generate the problem instances used in this study. In Section 4.2.3, we give the implementation details of the model. In Section 4.2.4, we describe the setup used for the experimental study. And lastly, in Section 4.2.5, we present and discuss the results of this study.
## 4.2.1 Instances

Benchmark instances for the MOBKP are generated according to the following procedure. For a given number of items $n \in \mathbb{Z}_{>0}$ and objectives $m \in \mathbb{Z}_{\ge 2}$, the item values $v_i^j \in [0, 1], i \in \{1, \dots, n\}, j \in \{1, \dots, m\}$, are generated according to a multi-variate uniform distribution of dimension $m$, defined by a positive-definite symmetric correlation matrix $C$ of size $m \times m$, such that $C_{aa} = 1$ and $C_{ab} = \rho_v$ for all $a, b \in \{1, \dots, m\}$ with $a \neq b$, that is

$$ \begin{pmatrix} 1 & \rho_v & \cdots & \rho_v \\ \rho_v & 1 & \cdots & \rho_v \\ \vdots & \vdots & \ddots & \vdots \\ \rho_v & \rho_v & \cdots & 1 \end{pmatrix} $$

where $\frac{6}{\pi} \sin^{-1} \frac{1}{2(1-m)} < \rho_v < 1$ denotes the correlation between any two distinct value vectors $v^j$ and $v^k$. We follow the procedure described by [73] to generate the value vectors under these conditions. In particular, we start by generating vectors $(V^1, \dots, V^m)$ following a multi-normal distribution of means 0 and correlation matrix $R = 2 \sin(\frac{\pi}{6}C)$. Then, the vectors $v^1 = \phi(V^1), \dots, v^m = \phi(V^m)$, such that $\phi(\cdot)$ is the univariate normal cumulative density function, are uniformly distributed with a correlation matrix $C$.

The weight vector $w$ is generated following a uniform distribution in the range $[0, 1]$ with approximate correlation $\rho_w$ to the auxiliary vector

$$ v' = \left( \sum_{j=1}^{m} v_{1}^{j}, \dots, \sum_{j=1}^{m} v_{n}^{j} \right) \qquad (4.8) $$

that gives the sum of the values for each item. Let $z(v')$ denote a function that centers and scales $v'$ to have mean 0 and standard deviation 1, and let $E \sim N(0, 1)$ be a vector of size $n$ drawn from a normal distribution with mean 0 and standard deviation 1. Then, we consider the following formulation to get an auxiliary vector $w'$ with correlation approximately $\rho_w$:

$$ w' = \rho_w z(v') + E \sqrt{1 - \rho_w^2} \qquad (4.9) $$

Note that, $w'$ is not drawn from a population with uniform distribution. However, we can run a Monte Carlo simulation to repeatedly resample from the same population as $w'$ and get an empirical cumulative density function. Then, we can get a vector $w$ drawn from a approximate uniform distribution with:

$$ w = \tilde{\phi}(w') \qquad (4.10) $$

where $\tilde{\phi}(w')$ denotes the empirical cumulative density function of the population of $w'$. The correlation between $w$ and $v'$ is approximately $\rho_w$.
Lastly, the weight constraint W is defined as a ratio of the sum of the items weights, with respect to a variable $\omega \in [0, 1]$ as follows:

$$W = \omega \sum_{i=1}^{n} w_i \qquad (4.11)$$

We remark that the commonly used benchmark instance classes described by Bazgan et al. [3] for two objectives, can be generalized for two and more objective using the procedure described above, in particular:

A. Uncorrelated instances: $\rho_v = 0$, $\rho_w = 0$, and $\omega = 0.5$.

B. Unconflicting instances: $\rho_v$ is set to a positive value, $\rho_w = 0$, and $\omega = 0.5$.

C. Conflicting instances: $\rho_v$ is set to a negative value, $\rho_w = 0$, and $\omega = 0.5$.

D. Conflicting instances with correlated weight: $\rho_v$ is set to a negative value, $\rho_w$ is set to a positive value, and $\omega = 0.5$.

For this experimental study, we will be considering instances generated with the following parameters:

* $m \in \{2, 3, 5\}$

* $n \in \mathcal{U}(50, 150)$

* $p_v \in \mathcal{U}(\frac{6}{\pi} \sin^{-1} \frac{1}{2(1-m)}, 1)$

* $p_w \in \mathcal{U}(-1, 1)$

* $\omega \in \mathcal{U}(0.3, 0.7)$

where $\mathcal{U}(\cdot)$ denotes the uniform distribution.

### 4.2.2 Anytime Performance Measure

To measure the performance of an empirical anytime performance profile P as a scalar value, we consider an anytime performance measure given by:

$$M(P) = \int_{q_l}^{q_u} \int_{0}^{t_u} P(t, q) \, dt \, dq \qquad (4.12)$$

where $q_l$ and $q_u$ denote the lower and upper bound for solution quality, and $t_u$ denotes the upper bound on time. This is a particular case of the anytime performance measure given by Jesus et al. [36], and related to the anytime performance measure given by López-Ibáñez and Stützle [51]. The former will be presented in more detail in Section 5.1.
### 4.2.3 Model Implementation Details

In this section we describe the implementation details of our model, which follows Framework I. We recall that according to this framework, the goal is to select a subset of problem instances from the training data set $\mathcal{I}$ that are expected to have similar anytime performance to the previously unseen problem instance $\iota \notin \mathcal{I}$ for which we want to predict anytime performance. Then, we take the performance traces collected a priori for the selected instances, and use them to build a performance profile that predicts the anytime performance of the previously unseen instance $\iota$.

To define the implementation of our model, we consider a training data set of problem instances to the MOBKP for which an anytime performance trace resulting from a single run has been collected for each algorithm. We also consider the following set of instance features that characterize a problem instance for the MOBKP as a set of measurable values:

* $\mathcal{F}_1$ — the number of decision variables $n$;

* $\mathcal{F}_2$ — the ratio between the capacity of the knapsack $W$, and the sum of the items' weights, i.e.:

$$ \mathcal{F}_2 = \frac{W}{\sum_{i=1}^{n} w_i} \qquad (4.13) $$

where $w_i$ denotes the weight of item $i \le i \le n$;

* $\mathcal{F}_3$ — the mean Spearman's correlation coefficient between all pairs of value vectors $v^i$, $1 \le i \le m$, i.e.:

$$ \mathcal{F}_3 = \frac{1}{m(m-1)/2} \sum_{i=1}^{m} \sum_{j=i+1}^{m} r_{v^i v^j} \qquad (4.14) $$

where $r_{v^i v^j}$ denotes the Spearman's correlation coefficient between value vectors $v^i$ and $v^j$. This feature relates to parameter $\rho_v$ that was used to generate the problem instance;

* $\mathcal{F}_4$ — the Spearman's correlation coefficient between the weights of the items $w$, and the sum of their values, i.e., $\sum_{i=1}^{m} v^i$. This feature relates to parameter $\rho_w$ that was used to generate the problem instance.

The empirical model then follows a k-nearest neighbor algorithm. In particular, we compute the distance between the features of the problem instances in the training data set and the features of the previously unseen instance $\iota$. The distance between features is considered with respect to the weighted Euclidean norm, i.e.:

$$ \| \mathcal{F}^l - \mathcal{F}^j \|_{\xi} = \sqrt{\sum_{\ell=1}^{4} (\mathcal{F}_{k}^{\ell} - \mathcal{F}_{k}^{j})^2 \xi_{k}} \qquad (4.15) $$
where $\xi_l > 0$ denotes the weight for the $l$-th feature, $\mathcal{F}^l$ denotes the feature set of the problem instance $l$ whose anytime performance we want to model, and $\mathcal{F}^j$ denotes the feature set of a problem instance $j$ in the training data set. To make it easier to reason about the weights we scale all features to the range $[0, 1]$ taking into account their range of possible values, i.e.:

* $\mathcal{F}_1 \in [50, 150]$

* $\mathcal{F}_2 \in [0.3, 0.7]$

* $\mathcal{F}_3 \in [-1/(m-1), 1]$

* $\mathcal{F}_4 \in [-1, 1]$

To define the weights $\xi$ for the distance metric, we consider the Spearman's correlation coefficient between the features and the anytime performance measure for the instances. In particular, let the $\rho_{\mathcal{F}_l}$ denote the Spearman's correlation coefficient for the instance feature $\mathcal{F}_l$, then weight $\xi_l$ is defined as:

$$ \xi_l = \exp(4 |\rho_{\mathcal{F}_l}|) \qquad (4.16) $$

We consider this definition because we want features that have a large absolute correlation value to have a significantly larger impact than features that have a small absolute correlation value. However, we do not want this difference to be too large, which is why we only multiply the absolute correlation value by 4 and not by a larger number, since features with a small absolute correlation value might still be relevant to distinguish between otherwise similar instances. Alternatively, we could set the weight values in a more informed manner, e.g., by using a grid search over possible values for the weights and taking into account the error for the model using different weights. However, this would result in much higher training times, so we opted to use the above definition instead since it gave good results in a preliminary analysis.

After computing the distances between each instance from the training data set and the previously unseen instance $l$, we select the $k$ instances from the training data set that minimize the distance. Then, we aggregate the anytime traces of the algorithm on those instances and build an anytime performance profile, which is expected to predict the anytime performance of instance $l$. Note that, the range of hypervolume values for different problem instances can be very different. As such, for the aggregation to be meaningful, we consider the relative hypervolume for measuring solution quality, which means that solution is defined in the range $[0, 1]$ for all anytime traces. Since we cannot always compute the non-dominated set in a feasible amount of time to have the maximal hypervolume, we use the maximum hypervolume found by any
of the algorithms we tested for a particular instance in order to compute the
relative hypervolume metric for each instance.

To validate our model, we split our known problem instances into two
distinct data sets: training and testing. In particular, we consider 80% of the
instances for the training data set, and 20% of the instances for the testing data
set. Then, we consider both numerical and visual analysis. For the former, we
consider the difference between the anytime measure of the predicted anytime
performance profile, and the anytime measure of the anytime trace gathered for
a problem instance *i* taken from the testing data set. To summarize the quality
of our model over all instances in the testing data set, we consider the *Mean
Absolute Error (MAE)*, which is given by the absolute difference between the
predicted and true anytime measures for each testing problem instance divided
by the number of testing problem instances. For visual analysis, we plot the
predicted performance profile and the true anytime trace.

Selecting an appropriate value for *k* in the *k*-nearest neighbor algorithm is
of the utmost importance to the quality of our model. In particular, if the value
for *k* is too small we are likely to select instances that are very similar in terms
of the instances features we considered. However, we might not have enough
anytime traces to build an anytime performance profile that is sufficiently
robust against random variations on anytime performance that may arise from
the stochastic nature of the algorithm, experimental conditions, or from other
instance features that have not been considered. On the other hand, if the
value for *k* is too large we might select instances that are too different from
the instance whose anytime performance we want to predict, which will result
in an anytime performance profile that is not specific enough.

To set parameter *k* for the model, we consider the *MAE* of our model given
by the *Leave-One-Out Cross-Validation (LOOCV)* method for different values
of *k*. In *LOOCV* we take the training data set, leave out one instance, and
use the remaining instances for training our model. Then, we use the trained
model to predict the anytime performance of the left-out instance. We repeat
this process such that every training instance is left out once, and compute the
*MAE* considering the predicted and true values for each instance. We compute
the *MAE* for each value of *k*, considering the same set of training problem
instances, and choose the value of *k* that gave the minimum *MAE* to build our
final model using the full training data set.

We remark that the choice of LOOCV, rather than other resampling meth-
ods, is due to its small bias when considering a small number of samples [57].
Nonetheless, preliminary experiments considering a 10-fold cross-validation
method yielded similar results. The choice of MAE, instead of Root-Mean-
Square Error (RMSE), which is also commonly used, is due to the fact that we expect that there are a small number of problem instances whose prediction will give large errors and cannot be significantly improved, and MAE is less sensitive to large errors since it does not square the errors.

The implementation of this model is available at [32].

### 4.2.4 Experimental Setup

To test our model we generated 500 instances for each number of objectives $m \in \{2, 3, 5\}$. For each instance, the value of $n$ was randomly sampled from the uniform distribution $\mathcal{U}(50, 150)$, the value of $\rho_v$ was randomly sampled from the uniform distribution $\mathcal{U}(\frac{6}{\pi}\sin^{-1}(\frac{1}{2(1-m)}), 1)$, the value of $\rho_w$ was randomly sampled from the uniform distribution $\mathcal{U}(-1, 1)$, and the value of $\omega$ was randomly sampled from the uniform distribution $\mathcal{U}(0.3, 0.7)$.

For each instance, we record the anytime performance trace, in terms of CPU-time and hypervolume, for one run of the BHV-DP and PLS algorithms for $m \in \{2, 3, 5\}$ objectives, and the GEPS algorithm for $m = 2$ objectives. As for the configuration of the algorithms, for the BHV-DP algorithm we consider the ordering of the items given by $O^{\min}$ since it gave the best overall results during preliminary experiments, which is consistent with the results reported by Bazgan et al. [3]. For the PLS algorithm, we consider the default ordering of the items, an initial solution corresponding to the empty knapsack, i.e., $x_i = 0$ for all $i \in \{1, \dots, n\}$, and a random selection of the next item to process. We chose this configuration since it often achieved a solution with good quality quickly. For the GEPS algorithm we considered a parameter $\ell = 100$, since it showed good results in Chapter 3, and in preliminary experiments for the MOBKP. The algorithms were implemented in C++ and are available at [33].

The algorithms were executed on a computer with two Intel(R) Xeon(R) Silver 4210R CPUs with clock frequency 2.40GHz, 10 cores, and 20 threads each. The algorithms were executed with a CPU time limit of 100 seconds and a memory limit of 8 Gb. Moreover, due to the large amount of time required to run all the experiments, the executions were carried out in parallel, 20 at a time (corresponding to one execution per CPU core). Although this is likely to impact the performance of the algorithms, we consider that it does not affect the results for our models, since both the training and testing instances are generated under the same conditions. In a real-world scenario, training instances should be generated under the same conditions that the new instances are expected to be solved on. The code used to generate the problem instances, and to gather the performance traces of the algorithms, is available at [32].
<figure><img src="image_42.png" /><figcaption>Figure 4.2: Relation between anytime measure and feature values for training problem instances. PLS, 2 objectives.</figcaption></figure>

Figure 4.2 shows the relation between feature values and anytime measure for each problem instance in the training data set. The blue line, given by the Locally Estimated Scatterplot Smoothing (LOESS) method [8], highlights the relationship between the feature values and the anytime measure. Table 4.1 gives the Spearman's correlation coefficient between each feature and the anytime measure of the training instances. These results suggest that anytime measure is most impacted by the correlation features, $F_3$ and $F_4$. The capacity feature, $F_2$, also appears to have an impact on anytime measure, but to a lesser degree. The problem size feature, $F_1$, appears to have no discernible impact on anytime measure. Table 4.2 gives the weights that were set according to Equation (4.16). We recall, that a higher value of $\xi_i$ implies that feature $F_i$ has a larger impact on the computed distance between instances.

Table 4.3 gives the MAE value obtained by LOOCV on the training data

<table><caption>Table 4.1: Spearman's correlation coefficient between the features and anytime measures of the training instances. PLS, 2 objectives.</caption><thead><tr><th>F<sub>1</sub></th><th>F<sub>2</sub></th><th>F<sub>3</sub></th><th>F<sub>4</sub></th></tr></thead><tbody><tr><td>-0.087</td><td>0.289</td><td>0.620</td><td>-0.504</td></tr></tbody></table>
<table><caption>Table 4.2: Weights for the distance metric. PLS, 2 objectives.</caption><thead><tr><th>ξ<sub>1</sub></th><th>ξ<sub>2</sub></th><th>ξ<sub>3</sub></th><th>ξ<sub>4</sub></th></tr></thead><tbody><tr><td>1.415</td><td>3.175</td><td>11.919</td><td>7.512</td></tr></tbody></table>

<figure><img src="image_43.png" /><figcaption>Figure 4.3: Histogram of the prediction error in terms of anytime measure for testing problem instances. The red line corresponds to the median. PLS, 2 objectives.</figcaption></figure>

set for each value of the parameter $k \in \{3, 4, \dots, 12\}$ that denotes how many instances the model selects to build the anytime performance profile that will predict anytime performance. We see that the best, i.e., minimal, MAE value is for $k = 11$. As such, we consider $k = 11$ to train the final model using all instances in the training data set.

<table><caption>Table 4.3: MAE given by LOOCV on the training problem instances. PLS, 2 objectives.</caption><thead><tr><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td>1.799</td><td>1.745</td><td>1.709</td><td>1.686</td><td>1.672</td><td>1.660</td><td>1.644</td><td>1.648</td><td><strong>1.643</strong></td><td>1.645</td></tr></tbody></table>

Figure 4.3 gives the distribution of the prediction errors for testing problem instances considering the final model trained with $k = 11$. In particular, the prediction error is considering in terms of the difference between the anytime measure of the predicted anytime performance profile for an instance in the
<figure><img src="image_44.png" /><figcaption>Figure 4.4: True and predicted anytime performance for 20 randomly selected testing problem instances. PLS, 2 objectives.</figcaption></figure>

testing data set, and the anytime measure of the true anytime performance trace for that instance. We see that most error values are close to zero, and the median (red line) is also very close to zero. We do see one particularly large error, but this is not surprising given the fact that the density plots shown earlier in this section revealed that there were outliers in the testing data set. Overall, we consider that these results indicate that our empirical model can accurately predict the anytime measure of the PLS algorithm for previously unseen problem instances with *m* = 2 objectives.

Figure 4.4 shows the true and predicted anytime performance for 20 ran-
domly selected testing problem instances. The black line corresponds to the
true anytime trace gathered for that problem instance, the red region gives
the range of anytime performance traces of the *k* problem instances selected
to build the anytime performance profile used for the prediction, and the red
line gives the median anytime trace of that predicted anytime performance
<figure><img src="image_45.png" /><figcaption>Figure 4.5: True and predicted anytime performance for testing problem instances that resulted in the largest prediction errors. PLS, 2 objectives.</figcaption></figure>

profile. We observe that the true anytime trace is often within the red region and close to the median anytime trace. Moreover, the red region is often quite small. These observations indicate that our model is often selecting problem instances with similar anytime behavior for building the predicted anytime performance profile. Therefore, it appears that our model can accurately predict the anytime performance of the PLS algorithm for previously unseen instances with *m* = 2 objectives quite well. The true anytime trace and predicted anytime performance of all 100 testing problem instances is presented in Appendix A.1.

Figure 4.5 and Table 4.4 give the predicted and true anytime performance profiles and measures for the 10 testing problem instances that had the largest prediction errors. We see that the top 3 largest prediction errors are for instances with a small true anytime measure as we had anticipated in the beginning of this section. Note that, due to the stochastic nature of the algorithm, and because we only gathered one run for each instance, these large errors might be a consequence of lucky, or unlucky, runs of the algorithm on particular instances. As such, one aspect to be studied in the future is the impact that collecting multiple runs per instance has. Finally, we see that for the last of these instances, the predicted performance profile is already quite close to the true anytime behavior. As such, we consider that our model can accurately predict the anytime performance of the PLS algorithm for problem instances with *m* = 2 objectives.
<table><caption>Table 4.4: True and predicted anytime measure values for testing problem instances that resulted in the largest prediction errors. PLS, 2 objectives.</caption><thead><tr><th>Instance</th><th>Predicted measure</th><th>True measure</th><th>Error</th></tr></thead><tbody><tr><td>207</td><td>90.13</td><td>68.53</td><td>21.60</td></tr><tr><td>335</td><td>89.34</td><td>82.97</td><td>6.37</td></tr><tr><td>343</td><td>93.47</td><td>87.21</td><td>6.26</td></tr><tr><td>142</td><td>95.56</td><td>90.57</td><td>4.99</td></tr><tr><td>51</td><td>89.35</td><td>94.29</td><td>-4.94</td></tr><tr><td>474</td><td>97.20</td><td>92.74</td><td>4.47</td></tr><tr><td>243</td><td>97.02</td><td>92.63</td><td>4.38</td></tr><tr><td>395</td><td>92.72</td><td>88.74</td><td>3.99</td></tr><tr><td>494</td><td>92.27</td><td>88.32</td><td>3.96</td></tr><tr><td>420</td><td>93.61</td><td>97.36</td><td>-3.74</td></tr></tbody></table>

## PLS — 3 Objectives

In this section, we analyze the results of predicting the anytime behavior of the PLS algorithm for problem instances with 3 objectives.

Figure 4.6 gives the density of anytime measure values for training and testing problem instances. We see that the distribution is similar for both types of instances, and that most instances have an anytime measure value greater than 90. Due to the small number of instances on the left tail, we expect that these might lead to erroneous predictions. We note that these results are similar to those of the previous section where we considered instances with 2 objectives. However, there are less significant outliers on the left tail.

Figure 4.7 gives the relation between the anytime measure and feature values for instances in the training data set, and Table 4.5 gives the corresponding Spearman's correlation coefficient between the features and anytime measure values. We see that the problem size feature, $F_1$, has the most meaningful impact on anytime measure. We also see that features $F_3$ and $F_4$ appear to have a relevant impact for some instances, but that there are many instances that always have an anytime measure value of 100 regardless of those features. This

<table><caption>Table 4.5: Spearman's correlation coefficient between the features and anytime measures of the training instances. PLS, 3 objectives.</caption><thead><tr><th>F<sub>1</sub></th><th>F<sub>2</sub></th><th>F<sub>3</sub></th><th>F<sub>4</sub></th></tr></thead><tbody><tr><td>0.559</td><td>0.208</td><td>-0.134</td><td>0.095</td></tr></tbody></table>
<figure><img src="image_46.png" /><figcaption>Figure 4.6: Anytime measure density plot for training and testing problem instances. PLS, 3 objectives.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>Figure 4.7: Relation between anytime measure and feature values for training problem instances. PLS, 3 objectives.</figcaption></figure>
<table><caption>Table 4.6: Weights for the distance metric. PLS, 3 objectives.</caption><thead><tr><th>ξ<sub>1</sub></th><th>ξ<sub>2</sub></th><th>ξ<sub>3</sub></th><th>ξ<sub>4</sub></th></tr></thead><tbody><tr><td>9.373</td><td>2.301</td><td>1.707</td><td>1.461</td></tr></tbody></table>

happens because for some instances, especially for small values of $n$, the BHV-DP algorithm can find the efficient set. Therefore, in those cases, the anytime measure value for the PLS algorithm will be relative to the hypervolume of the non-dominated set. However, for other instances, the BHV-DP cannot find the efficient set within the allotted time, and the best known approximation, in terms of hypervolume, might be given by the PLS algorithm. Therefore, in those case, the anytime measure value for the PLS algorithm will be relative to its own approximation, and will have a value of 1. As a result, it seems that $F_1$ is the feature that best indicates whether or not the PLS algorithm found the best known approximation, and the correlation features, $F_3$ and $F_4$, better explain the anytime measure for instances where it did not. Note that, in the previous case with 2 objectives, the BHV-DP algorithm always finds the efficient set within the allotted time, so we did not see this behavior. Table 4.6, gives the weights obtained from the corresponding correlation coefficients according to Equation (4.16).

Table 4.7 gives the MAE values obtained with the LOOCV method for each value of the parameter $k$ considered. Since the minimum MAE value was found for $k = 9$, we will consider that value to train our final model on the complete training data set.

Figure 4.8 shows the distribution of the prediction errors in terms of anytime measure value for testing problem instances using the final model. We see that the error for most instances, as well as the median error, are close to zero. Still, it seems that there are slightly more errors with negative values, which means that the model is under-estimating the true anytime measure. This is in contrast to the expected result that the largest errors would be related to the outliers in the left tail, which is what we saw in the previous section for instances with 2 objectives. One reason for this, might be related to what was discussed previously that for some instances the best known approximation

<table><caption>Table 4.7: MAE given by LOOCV on the training problem instances. PLS, 3 objectives.</caption><thead><tr><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td>1.801</td><td>1.820</td><td>1.830</td><td>1.816</td><td>1.821</td><td>1.792</td><td><strong>1.784</strong></td><td>1.811</td><td>1.840</td><td>1.866</td></tr></tbody></table>
<figure><img src="image_48.png" /><figcaption>Figure 4.8: Histogram of the prediction error in terms of anytime measure for testing problem instances. The red line corresponds to the median. PLS, 3 objectives.</figcaption></figure>

is given by the *BHV-DP* algorithm, whereas for others it is given by the *PLS* algorithm. As a results of this, for similar instances, it might happen that the *BHV-DP* was able to find the efficient set within the allotted time, but for others it did not. As such, we might have similar instances with distinct relative anytime measure values, i.e., close to 100 if *PLS* found the best known solution, and significantly below 100 if *BHV-DP* found the optimal solution. Despite this, the results show that the model can accurately predict the anytime measure obtained by the *PLS* algorithm for instances with 3 objectives in most cases.

Figure 4.9 gives the true and predicted anytime performance for 20 randomly selected testing problem instances. We see that the results are generally quite good, and that the model can accurately predict the anytime performance for previously unseen instances. Nonetheless, there are some notable outliers, e.g., instances 148, 165, 193, 252, and 305. However, Table 4.8 shows that, with the exception of instance 193, these outliers are in the list of the 10 worst predictions. Thus, we were a bit unlucky with the random sample, and looking at the full results in Appendix A.2, it is clear that the model can very often predict the anytime performance of the *PLS* algorithm for problem instances with 3 objectives.

Figure 4.10 and Table 4.8 show the predicted and true anytime performance and anytime measure for the 10 testing problem instances with the worst pre-
<figure><img src="image_49.png" /><figcaption>Figure 4.9: True and predicted anytime performance for 20 randomly selected testing problem instances. PLS, 3 objectives.</figcaption></figure>

diction error. Some of the largest errors are for instances with a true anytime measure very close to 100, which is consistent with the analysis regarding instances whose best known solution may be found by the PLS or BHV-DP algorithm. There are also some errors for instances with a small true anytime measure, which is to be expected due to the small number of instances with those anytime measure values. Nonetheless, we see that the overall anytime behavior for these instances is not very far from the true anytime trace.
<figure><img src="image_50.png" /><figcaption>Figure 4.10: True and predicted anytime performance for testing problem instances that resulted in the largest prediction errors. PLS, 3 objectives.</figcaption></figure>

<table><caption>Table 4.8: True and predicted anytime measure values for testing problem instances that resulted in the largest prediction errors. PLS, 3 objectives.</caption><thead><tr><td>Instance</td><td>Predicted measure</td><td>True measure</td><td>Error</td></tr></thead><tbody><tr><td>308</td><td>92.04</td><td>99.16</td><td>-7.12</td></tr><tr><td>305</td><td>92.23</td><td>85.22</td><td>7.01</td></tr><tr><td>289</td><td>93.02</td><td>99.98</td><td>-6.96</td></tr><tr><td>165</td><td>90.14</td><td>83.65</td><td>6.48</td></tr><tr><td>282</td><td>97.17</td><td>91.75</td><td>5.42</td></tr><tr><td>252</td><td>96.72</td><td>91.59</td><td>5.13</td></tr><tr><td>159</td><td>90.16</td><td>95.22</td><td>-5.06</td></tr><tr><td>461</td><td>94.93</td><td>99.98</td><td>-5.05</td></tr><tr><td>148</td><td>96.26</td><td>91.56</td><td>4.70</td></tr><tr><td>130</td><td>96.86</td><td>92.19</td><td>4.67</td></tr></tbody></table>
PLS — 5 Objectives

In this section, we analyze the results of predicting the anytime behavior of
the PLS algorithm for problem instances with 5 objectives.

Figure 4.11 shows the density of the anytime measure for testing and train-
ing problem instances. We see a similar density for the two data sets. However,
the testing data set seems to have more instances with an anytime measure
around 90, and less instances with an anytime measure close to 100. As such,
we expect to be able to generally predict anytime performance, but that our
largest errors might be for instances with an anytime measure close to 90.
Compared to the two previous cases where we considered instances with 2
and 3 objectives, the anytime measure values are much closer to 100, i.e., the
maximum possible value. This is a consequence of the best known solution
for instances with 5 objectives being very often found by the PLS algorithm,
instead of the BHV-DP algorithm, since the latter is not often able to find the
efficient set in the given amount of time.

Figure 4.12 gives the relation between the features values and anytime mea-
sure for instances in the training data set, and Table 4.9 gives the corresponding
Spearman’s correlation coefficient. We see that, unlike in the previous cases,
all features have a give a small correlation coefficient, and as such the weights
for the distance metric, shown in Table 4.10, are not too different from each
other.

<figure><img src="image_51.png" /><figcaption>Figure 4.11: Anytime measure density plot for training and testing problem instances. PLS, 5 objectives.</figcaption></figure>
<figure><img src="image_52.png" /><figcaption>Figure 4.12: Relation between anytime measure and feature values for training problem instances. PLS, 5 objectives.</figcaption></figure>

<table><caption>Table 4.9: Spearman's correlation coefficient between the features and anytime measures of the training instances. PLS, 5 objectives.</caption><thead><tr><th>F<sub>1</sub></th><th>F<sub>2</sub></th><th>F<sub>3</sub></th><th>F<sub>4</sub></th></tr></thead><tbody><tr><td>0.130</td><td>0.115</td><td>-0.063</td><td>-0.092</td></tr></tbody></table>

<table><caption>Table 4.10: Weights for the distance metric. PLS, 5 objectives.</caption><thead><tr><th>ξ<sub>1</sub></th><th>ξ<sub>2</sub></th><th>ξ<sub>3</sub></th><th>ξ<sub>4</sub></th></tr></thead><tbody><tr><td>1.680</td><td>1.583</td><td>1.286</td><td>1.445</td></tr></tbody></table>
Table 4.11 gives the MAE values obtained by LOOCV on the training problem instances for values of parameter *k* = 3, 4, ..., 12. The best result was found for *k* = 5, which we use to build the final model.

Figure 4.13 gives the histogram of prediction errors for the testing problem instances. We see that most errors, and the median, are close to zero, and that, as expected, the most significant errors are positive, i.e., a consequence of over-estimating the anytime performance. Nonetheless, we see that the model can accurately predict the anytime measure values very often.

Figure 4.14 shows the true and predicted anytime performance of 20 randomly selected testing problem instances. We see that for most instances the prediction is quite good, but there are some outliers. Most notably, instance 412 has a very significant error, but it is worth noting that it is in the fact the testing instance that has the largest error and a clear outlier. Compared to the previous scenarios with 2 and 3 objectives, we seem to have larger red regions,

<table><caption>Table 4.11: MAE given by LOOCV on the training problem instances. PLS, 5 objectives.</caption><thead><tr><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td>1.530</td><td>1.525</td><td><strong>1.475</strong></td><td>1.521</td><td>1.546</td><td>1.555</td><td>1.559</td><td>1.539</td><td>1.534</td><td>1.536</td></tr></tbody></table>

<figure><img src="image_53.png" /><figcaption>Figure 4.13: Histogram of the prediction error in terms of anytime measure for testing problem instances. The red line corresponds to the median. PLS, 5 objectives.</figcaption></figure>
<figure><img src="image_54.png" /><figcaption>Figure 4.14: True and predicted anytime performance for 20 randomly selected testing problem instances. PLS, 5 objectives.</figcaption></figure>

which corresponds to the region of anytime traces of the training instances selected for prediction. This suggests that it is harder to find instances with similar anytime behavior for a larger number of objectives, which is not too surprising. The results for every testing problem instance can be found in Appendix A.3.

Figure 4.15 and Table 4.12 give the true and predicted anytime measure and performance for the testing problem instances with the largest prediction error. As expected, most of the largest errors are for instances with a true measure close to 90. Nonetheless, we see that the model can already predict the anytime performance accurately for the last 5 of these instances, and also taking into account the full results, we see that our model can accurately predict anytime behavior often.
<figure><img src="image_55.png" /><figcaption>Figure 4.15: True and predicted anytime performance for testing problem instances that resulted in the largest prediction errors. PLS, 5 objectives.</figcaption></figure>

<table><caption>Table 4.12: True and predicted anytime measure values for testing problem instances that resulted in the largest prediction errors. PLS, 5 objectives.</caption><thead><tr><td>Instance</td><td>Predicted measure</td><td>True measure</td><td>Error</td></tr></thead><tbody><tr><td>412</td><td>95.64</td><td>69.40</td><td>26.23</td></tr><tr><td>461</td><td>95.42</td><td>88.46</td><td>6.95</td></tr><tr><td>317</td><td>99.75</td><td>93.25</td><td>6.50</td></tr><tr><td>465</td><td>96.96</td><td>90.79</td><td>6.16</td></tr><tr><td>13</td><td>95.29</td><td>89.71</td><td>5.59</td></tr><tr><td>78</td><td>96.28</td><td>90.82</td><td>5.46</td></tr><tr><td>322</td><td>95.25</td><td>99.97</td><td>-4.73</td></tr><tr><td>370</td><td>96.53</td><td>91.98</td><td>4.55</td></tr><tr><td>490</td><td>93.24</td><td>97.43</td><td>-4.18</td></tr><tr><td>7</td><td>96.19</td><td>99.97</td><td>-3.79</td></tr></tbody></table>
## BHV-DP – 2 Objectives

In this section, we analyze the prediction of anytime performance of the BHV-DP algorithm for problem instances with 2 objectives.

Figure 4.16 gives the density of anytime measure for problem instances in the testing training data sets. We can see that the anytime measure is very close to 100 for most instances. This is due to the fact that the BHV-DP algorithm can often find the efficient set quickly for the instances and time budget we considered. However, there are some testing instances with an anytime measure below 97.5 that we expect to lead to large errors since there are very little training instances with similar values.

Figure 4.17 shows the relation between the instance features and the anytime measure values, and Table 4.13 gives the Spearman’s correlation coefficients. Instance features $F_1$ and $F_3$ seem to have the most significant impact on anytime measure values, followed by $F_4$. On the other hand, $F_2$ seems to have little or no impact. Table 4.14 gives the weights for the distance metric computed from the correlation values according to Equation (4.16).

Table 4.15 gives the MAE obtained by LOOCV for values of parameter $k = 3, 4, ..., 12$. The best value was found for $k = 9$, which consequently is the value that we will use to build the final model considering all problem instances in the training data set.

<figure><img src="image_56.png" /><figcaption>Figure 4.16: Anytime measure density plot for training and testing problem instances. BHV-DP, 2 objectives.</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>Figure 4.17: Relation between anytime measure and feature values for training problem instances. BHV-DP, 2 objectives.</figcaption></figure>

<table><caption>Table 4.13: Spearman's correlation coefficient between the features and anytime measures of the training instances. BHV-DP, 2 objectives.</caption><thead><tr><th>F<sub>1</sub></th><th>F<sub>2</sub></th><th>F<sub>3</sub></th><th>F<sub>4</sub></th></tr></thead><tbody><tr><td>-0.576</td><td>-0.037</td><td>0.658</td><td>-0.356</td></tr></tbody></table>

<table><caption>Table 4.14: Weights for the distance metric. BHV-DP, 2 objectives.</caption><thead><tr><th>ξ<sub>1</sub></th><th>ξ<sub>2</sub></th><th>ξ<sub>3</sub></th><th>ξ<sub>4</sub></th></tr></thead><tbody><tr><td>10.016</td><td>1.158</td><td>13.907</td><td>4.149</td></tr></tbody></table>

<table><caption>Table 4.15: MAE given by LOOCV on the training problem instances. BHV-DP, 2 objectives.</caption><thead><tr><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td>0.252</td><td>0.249</td><td>0.242</td><td>0.238</td><td>0.234</td><td>0.233</td><td><strong>0.226</strong></td><td>0.231</td><td>0.226</td><td>0.229</td></tr></tbody></table>
<figure><img src="image_58.png" /><figcaption>Figure 4.18: Histogram of the prediction error in terms of anytime measure for testing problem instances. The red line corresponds to the median. BHV-DP, 2 objectives.</figcaption></figure>

Figure 4.18 shows the prediction errors, in terms of anytime measure, over all 100 testing problem instances. The errors are for the most part quite small, which is expected given that most instances have a similar anytime measure value as discussed earlier. We also see some significant positive errors that, as expected, correspond to instances that have a small anytime measure value.

Figure 4.19 gives the true and predicted anytime performance for 20 randomly selected testing problem instances. The results for all 100 testing problem instances can be found in Appendix A.4. We note that the red region of anytime traces used to predict the anytime performance of each instance is generally quite large compared what was seen before when predicting the anytime behavior of the PLS algorithm. This suggests that accurately predicting the anytime performance of the BHV-DP algorithm is more difficult, and we may need to consider more instances or instance features for the model to be more precise. Still, the tendency of anytime behavior is correctly predicted in most cases. As such, we consider that the model can accurately predict the anytime behavior of the BHV-DP algorithm for instances with 2 objectives.

Figure 4.20 and Table 4.16 show the true and predicted anytime measure and anytime performance for the testing problem instances with the largest prediction errors. We see that, as expected, the worst predictions are for instances that resulted in a small true anytime measure value.
<figure><img src="image_59.png" /><figcaption>Figure 4.19: True and predicted anytime performance for 20 randomly selected testing problem instances. BHV-DP, 2 objectives.</figcaption></figure>

<table><caption>Table 4.16: True and predicted anytime measure values for testing problem instances that resulted in the largest prediction errors. BHV-DP, 2 objectives.</caption><thead><tr><td>Instance</td><td>Predicted measure</td><td>True measure</td><td>Error</td></tr></thead><tbody><tr><td>472</td><td>96.06</td><td>88.61</td><td>7.45</td></tr><tr><td>343</td><td>99.42</td><td>95.04</td><td>4.37</td></tr><tr><td>496</td><td>96.19</td><td>92.16</td><td>4.03</td></tr><tr><td>159</td><td>98.89</td><td>95.93</td><td>2.96</td></tr><tr><td>424</td><td>98.87</td><td>96.90</td><td>1.97</td></tr><tr><td>467</td><td>96.66</td><td>95.36</td><td>1.30</td></tr><tr><td>30</td><td>98.71</td><td>97.98</td><td>0.73</td></tr><tr><td>346</td><td>99.21</td><td>99.68</td><td>-0.47</td></tr><tr><td>294</td><td>97.73</td><td>98.08</td><td>-0.35</td></tr><tr><td>440</td><td>99.40</td><td>99.06</td><td>0.34</td></tr></tbody></table>
<figure><img src="image_60.png" /><figcaption>Figure 4.20: True and predicted anytime performance for testing problem instances that resulted in the largest prediction errors. BHV-DP, 2 objectives.</figcaption></figure>

## BHV-DP – 3 Objectives

In this section we analyze the prediction of anytime performance of the BHV-DP algorithm for problem instances with 3 objectives. Figure 4.21 gives the distribution of anytime measure values for the training and testing problem instances. We see a similar distribution between both data sets, but a large range of anytime measure values. This large range of values is explained by the fact that the BHV-DP can sometimes take a long time to find a good approximation on the considered instances. Given that the number of instances in the training data set is not very extensive when taking into account the combination of parameters, having a large range of anytime measure values may suggest that there are not many instances with similar anytime behavior and instance features, which could result in degraded prediction quality.

Figure 4.22 gives the relation between the feature values and anytime measure values for the training problem instances, and Table 4.17 gives the corresponding Spearman's correlation coefficient. The results indicate that all

<table><caption>Table 4.17: Spearman's correlation coefficient between the features and anytime measures of the training instances. BHV-DP, 3 objectives.</caption><thead><tr><th>F<sub>1</sub></th><th>F<sub>2</sub></th><th>F<sub>3</sub></th><th>F<sub>4</sub></th></tr></thead><tbody><tr><td>-0.479</td><td>-0.247</td><td>0.684</td><td>-0.313</td></tr></tbody></table>
<figure><img src="image_61.png" /><figcaption>Figure 4.21: Anytime measure density plot for training and testing problem instances. BHV-DP, 3 objectives.</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>Figure 4.22: Relation between anytime measure and feature values for training problem instances. BHV-DP, 3 objectives.</figcaption></figure>
<table><caption>Table 4.18: Weights for the distance metric. BHV-DP, 3 objectives.</caption><thead><tr><th>ξ<sub>1</sub></th><th>ξ<sub>2</sub></th><th>ξ<sub>3</sub></th><th>ξ<sub>4</sub></th></tr></thead><tbody><tr><td>6.782</td><td>2.684</td><td>15.411</td><td>3.494</td></tr></tbody></table>

<table><caption>Table 4.19: MAE given by LOOCV on the training problem instances. BHV-DP, 3 objectives.</caption><thead><tr><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td><strong>4.908</strong></td><td>5.113</td><td>5.088</td><td>5.173</td><td>5.193</td><td>5.233</td><td>5.298</td><td>5.333</td><td>5.442</td><td>5.467</td></tr></tbody></table>

instance features have some impact on the anytime measure value, with the most significant being $F_3$ and $F_1$. Table 4.18 gives the weight values for the distance metric obtained with Equation (4.16) from the correlation coefficients.

Table 4.19 gives the MAE values obtained by LOOCV for $k = 3, 4, \dots, 12$. The best value was found for $k = 3$, which we use to build the final model.

Figure 4.23 gives the histogram of prediction errors for the testing problem instances. We can see that most errors are close to zero. However, we see a considerable number of instances with significant errors, which is not unexpected given the large range of anytime measure values discussed before. We expect

<figure><img src="image_63.png" /><figcaption>Figure 4.23: Histogram of the prediction error in terms of anytime measure for testing problem instances. The red line corresponds to the median. BHV-DP, 3 objectives.</figcaption></figure>
<figure><img src="image_64.png" /><figcaption>Figure 4.24: True and predicted anytime performance for 20 randomly selected testing problem instances. BHV-DP, 3 objectives.</figcaption></figure>

that by considering a larger number of training instances the results from our model could be significantly improved. Nonetheless, as it stands, the quality of prediction in terms of anytime measure is not very good for this case.

Figure 4.24 gives the true and predicted anytime performance for 20 randomly selected testing instances. The results for every testing problem instance can be found in Appendix A.5. As in the previous case, we see that the region of anytime traces, i.e., the red region, is often large despite the fact that we are selecting *k* = 3 instances. This indicates that the model fails to find instances that have similar anytime performance, which, as discussed earlier, is expected due to the large range of anytime measure values. Moreover, we see that in quite a few cases the prediction is not very good. However, on a more positive note, we see that the overall tendency of anytime behavior can often be predicted. This suggests that by considering more instances the model could better predict anytime behavior.
<figure><img src="image_65.png" /><figcaption>Figure 4.25: True and predicted anytime performance for testing problem instances that resulted in the largest prediction errors. BHV-DP, 3 objectives.</figcaption></figure>

<table><caption>Table 4.20: True and predicted anytime measure values for testing problem instances that resulted in the largest prediction errors. BHV-DP, 3 objectives.</caption><thead><tr><td>Instance</td><td>Predicted measure</td><td>True measure</td><td>Error</td></tr></thead><tbody><tr><td>402</td><td>78.92</td><td>42.91</td><td>36.01</td></tr><tr><td>349</td><td>54.12</td><td>26.67</td><td>27.45</td></tr><tr><td>204</td><td>68.01</td><td>42.91</td><td>25.11</td></tr><tr><td>346</td><td>85.88</td><td>61.20</td><td>24.68</td></tr><tr><td>17</td><td>56.65</td><td>32.91</td><td>23.74</td></tr><tr><td>470</td><td>79.33</td><td>57.35</td><td>21.99</td></tr><tr><td>289</td><td>83.33</td><td>61.72</td><td>21.61</td></tr><tr><td>342</td><td>74.42</td><td>53.32</td><td>21.09</td></tr><tr><td>298</td><td>75.30</td><td>55.63</td><td>19.67</td></tr><tr><td>452</td><td>53.31</td><td>34.62</td><td>18.69</td></tr></tbody></table>
Figure 4.25 and Table 4.20 give the true and predicted anytime measure and performance for the testing problem instances with the largest prediction errors. We see that often, the true anytime trace lies on the boundary of the red region. Moreover, we see that the largest errors happen for instances whose true anytime measure is not very large, despite there being more instances with large anytime measure values as seen earlier in the density plot. This reinforces the idea that by considering more instances the quality of the prediction could be significantly improved.

## BHV-DP — 5 Objectives

In this section, we analyze the prediction of anytime performance of the BHV-DP algorithm for problem instances with 5 objectives. Figure 4.26 gives the distribution of anytime measure values for the training and testing problem instances. We see a similar distribution between both data sets, and that the generated instances cover the whole range of possible anytime measure values. As in the previous case, given the fact that we do not have a very extensive set of training instances, this is likely to result in large prediction errors.

Figure 4.27 gives the relation between the instance features and anytime measure on training problem instances, and Table 4.21 gives the corresponding Spearman's correlation coefficient. We see that all features have some

<figure><img src="image_66.png" /><figcaption>Figure 4.26: Anytime measure density plot for training and testing problem instances. BHV-DP, 5 objectives.</figcaption></figure>
<figure><img src="image_67.png" /><figcaption>Figure 4.27: Relation between anytime measure and feature values for training problem instances. BHV-DP, 5 objectives.</figcaption></figure>

<table><caption>Table 4.21: Spearman's correlation coefficient between the features and anytime measures of the training instances. BHV-DP, 5 objectives.</caption><thead><tr><th>F<sub>1</sub></th><th>F<sub>2</sub></th><th>F<sub>3</sub></th><th>F<sub>4</sub></th></tr></thead><tbody><tr><td>-0.489</td><td>-0.233</td><td>0.719</td><td>-0.322</td></tr></tbody></table>

<table><caption>Table 4.22: Weights for the distance metric. BHV-DP, 5 objectives.</caption><thead><tr><th>ξ<sub>1</sub></th><th>ξ<sub>2</sub></th><th>ξ<sub>3</sub></th><th>ξ<sub>4</sub></th></tr></thead><tbody><tr><td>7.061</td><td>2.542</td><td>17.735</td><td>3.630</td></tr></tbody></table>

impact on anytime measure, but features $F_3$ and $F_1$ are the ones with the most significant impact. Table 4.22 gives the weights calculated with Equation (4.16).

Table 4.23 gives the MAE obtained by LOOCV for $k = 3, 4, \dots, 12$. The best value was found for $k = 4$, which we use to build our final model.

Figure 4.28 gives the prediction errors for the testing problem instances. We see that the median and largest concentration of errors is close to zero. However, there are many large prediction errors. This is not unexpected, since a larger number of objectives will give rise to more distinct behaviors. There-
<table><caption>Table 4.23: MAE given by LOOCV on the training problem instances. BHV-DP, 5 objectives.</caption><thead><tr><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td>8.686</td><td><strong>8.636</strong></td><td>8.859</td><td>8.794</td><td>8.768</td><td>9.022</td><td>8.963</td><td>8.979</td><td>8.917</td><td>9.105</td></tr></tbody></table>

<figure><img src="image_68.png" /><figcaption>Figure 4.28: Histogram of the prediction error in terms of anytime measure for testing problem instances. The red line corresponds to the median. BHV-DP, 5 objectives.</figcaption></figure>

fore, we consider that more instances or features are needed to precisely predict anytime measure for the BHV-DP algorithm on instances with 5 objectives.

Figure 4.29 gives the true and predicted anytime performance of 20 randomly selected testing problem instances. The results for all 100 testing problem instances are shown in Appendix A.6. We see that the red region of anytime traces that make up the predicted anytime performance profile is often quite large, and in a few cases the true tendency of anytime performance deviates significantly from the predicted anytime performance. As such, we consider that the overall prediction is not very good and could be further improved.

Figure 4.30 and Table 4.24 show the true and predicted anytime measure and performance profile for testing problem instances with the largest prediction errors. It is quite clear that the model is selecting instances that have quite different anytime behavior, so more information is needed to accurately predict anytime performance.
<figure><img src="image_69.png" /><figcaption>Figure 4.29: True and predicted anytime performance for 20 randomly selected testing problem instances. BHV-DP, 5 objectives.</figcaption></figure>

<table><caption>Table 4.24: True and predicted anytime measure values for testing problem instances that resulted in the largest prediction errors. BHV-DP, 5 objectives.</caption><thead><tr><td>Instance</td><td>Predicted measure</td><td>True measure</td><td>Error</td></tr></thead><tbody><tr><td>288</td><td>85.92</td><td>35.85</td><td>50.07</td></tr><tr><td>317</td><td>51.82</td><td>98.50</td><td>-46.67</td></tr><tr><td>102</td><td>48.17</td><td>94.18</td><td>-46.01</td></tr><tr><td>400</td><td>56.81</td><td>19.42</td><td>37.39</td></tr><tr><td>70</td><td>62.19</td><td>24.99</td><td>37.20</td></tr><tr><td>465</td><td>61.92</td><td>98.17</td><td>-36.25</td></tr><tr><td>13</td><td>58.40</td><td>93.52</td><td>-35.12</td></tr><tr><td>217</td><td>52.65</td><td>84.51</td><td>-31.86</td></tr><tr><td>188</td><td>29.85</td><td>3.05</td><td>26.80</td></tr><tr><td>226</td><td>51.14</td><td>74.61</td><td>-23.47</td></tr></tbody></table>
<figure><img src="image_70.png" /><figcaption>Figure 4.30: True and predicted anytime performance for testing problem instances that resulted in the largest prediction errors. BHV-DP, 5 objectives.</figcaption></figure>

GEPS — 2 Objectives

In this section, we conclude our experimental analysis by looking at the prediction of anytime performance for the GEPS algorithm on problem instances with 2 objectives.

Figure 4.31 gives the density plot of anytime measure for training and testing problem instances. We can see that virtually all instances have an anytime measure close to 100. For this reason, we show in Figure 4.32, the density of instances with an anytime measure value greater than 99.9. We can see that the density is similar for both data sets. Still, we see a small number of testing problem instances below 99.9, but no such instances in the training data set. As such, we expect our model to give large prediction errors for these.

Figure 4.33 shows the relation between the instance features and anytime measure for the training problem instances, and Figure 4.34 shows the same plot for anytime measure values greater than 99.9. Moreover, Table 4.25 gives the Spearman’s correlation coefficient, on all training instances. We see that instance feature $F_3$ has the most significant impact on anytime measure. Then, features $F_1$ and $F_4$ also have some impact but to a lesser degree, and feature $F_2$ has little impact. Table 4.26 gives the weights calculated with Equation (4.16) according to the correlation coefficients.

Table 4.27 shows the MAE values obtained by the LOOCV method for $k = 3, 4, ..., 12$. Thus, we choose $k = 7$ to build the final model as it has the best MAE value.
<figure><img src="image_71.png" /><figcaption>Figure 4.31: Anytime measure density plot for training and testing problem instances. GEPS, 2 objectives.</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>Figure 4.32: Anytime measure density plot for training and testing problem instances, for anytime measures greater than 99.9. GEPS, 2 objectives.</figcaption></figure>