Nawal Benabbou
email: nawal.benabbou@lip6.fr

Cassandre Leroy
email: cassandre.leroy@lip6.fr

Thibaut Lust
email: thibaut.lust@lip6.fr

An Interactive Regret-Based Genetic Algorithm for Solving Multi-Objective Combinatorial Optimization Problems

We propose a new approach consisting in combining genetic algorithms and regret-based incremental preference elicitation for solving multi-objective combinatorial optimization problems with unknown preferences. For the purpose of elicitation, we assume that the decision maker's preferences can be represented by a parameterized scalarizing function but the parameters are initially not known. Instead, the parameter imprecision is progressively reduced by asking preference queries to the decision maker during the search to help identify the best solutions within a population. Our algorithm, called RIGA, can be applied to any multi-objective combinatorial optimization problem provided that the scalarizing function is linear in its parameters and that a (near-)optimal solution can be efficiently determined when preferences are known. Moreover, RIGA runs in polynomial time while asking no more than a polynomial number of queries. For the multi-objective traveling salesman problem, we provide numerical results showing its practical efficiency in terms of number of queries, computation time and gap to optimality.

Introduction

Preference-based search is an important line of research in Artificial Intelligence with various applications to investments planning, resource allocation, group configuration, recommendation systems and electronic commerce (see e.g., (Junker 2002;[START_REF] Auger | Articulating user preferences in manyobjective problems by sampling the weighted hypervolume[END_REF][START_REF] Roijers | Multi-objective decision making[END_REF]Marinescu, Razak, and Wilson 2017;[START_REF] Zintgraf | Ordered preference elicitation strategies for supporting multi-objective decision making[END_REF]). The increasing complexity of applications encountered in Artificial Intelligence and other areas of Computer Science significantly complicates the task of decision makers (DM) who need to find the best solution to their problems. In many practical search problems, the evaluation of solutions involves several aspects or points of view corresponding to multiple objectives (e.g., time, cost, distance, risk). Multi-objective combinatorial optimization (MOCO) is concerned with problems involving several (conflicting) objectives/criteria to be optimized simultaneously (e.g., minimizing costs while maximizing profits).

In MOCO problems, without preference information, we only know that the best solution for the decision maker (DM) is among the Pareto-optimal solutions1 . While the first contributions mainly focused on the generation of all Pareto optimal solutions, tackling problems with many objectives has required the development of new methods to avoid computational overhead due to the high number of Pareto optimal solutions. Many of these new contributions concern the integration of preferences. However most of them focus on the elaboration of efficient methods to identify the best solution given a preference model, considering preference elicitation as a preliminary stage. In this paper, we address both issues simultaneously. The idea is to collect preference information during the search to focus on promising solutions while reducing the elicitation burden by only asking queries to discriminate between the solutions found during the search.

Preference elicitation on combinatorial domains is an active topic that has been recently studied in various contexts, e.g. in voting problems [START_REF] Benabbou | Solving multi-agent knapsack problems using incremental approval voting[END_REF][START_REF] Benabbou | Solving multi-agent knapsack problems using incremental approval voting[END_REF], in stable matching problems [START_REF] Drummond | Preference elicitation and interview minimization in stable matchings[END_REF], in constraint satisfaction problems [START_REF] Boutilier | Constraint-based optimization and utility elicitation using the minimax decision criterion[END_REF], in Markov Decision Processes [START_REF] Regan | Eliciting additive reward functions for markov decision processes[END_REF][START_REF] Weng | Interactive value iteration for markov decision processes with unknown rewards[END_REF][START_REF] Gilbert | Reducing the number of queries in interactive value iteration[END_REF] and in multi-objective optimization problems (Korhonen 2005;[START_REF] Kaddani | Weighted sum model with partial preference information: application to multi-objective optimization[END_REF][START_REF] Benabbou | Interactive resolution of multiobjective combinatorial optimization problems by incremental elicitation of criteria weights[END_REF]. In this paper, we assume that the DM's preferences can be represented by a parameterized scalarizing function that is initially not known, and we study the potential of incremental elicitation (White III, Sage, and Dozono 1984) in this context. The idea is to select preference queries one at a time, to be as informative as possible, so as to progressively reduce the set of admissible preference parameters until determining a (near-)optimal solution. We focus here on regret-based incremental elicitation, an approach recently designed within the Artificial Intelligence community, as it was proved to be very effective in practice [START_REF] Wang | Incremental Utility Elicitation with the Minimax Regret Decision Criterion[END_REF][START_REF] Boutilier | Constraint-based optimization and utility elicitation using the minimax decision criterion[END_REF][START_REF] Benabbou | Incremental weight elicitation for multiobjective state space search[END_REF]. More precisely, we propose a new interactive solving method that uses regret-based incremental elicitation techniques to solve hard MOCO problems with unknown preference parameters. For this purpose, we focus on genetic algorithms, which are widely applied to nu-merous hard combinatorial optimization problems.

The paper is organized as follows: Firstly, we conduct a literature overview on interactive genetic algorithms and regret-based solving methods. Secondly, we recall the general principle of incremental elicitation driven by the minimax regret decision criterion. Thirdly, a new regret-based interactive genetic algorithm is proposed and illustrated on a small instance of the Multi-objective Traveling Salesman Problem (MTSP). Finally, we provide numerical tests showing its practical efficiency, comparing its performances with that of several existing methods.

Related Works

In the past decade, integrating preferences into genetic algorithms has become increasingly popular, see e.g., (Branke et al. 2008;[START_REF] Xin | Interactive multiobjective optimization: A review of the state-of-the-art[END_REF]. Due to the high number of different interactive methods that has been developed, [START_REF] Xin | [END_REF] have recently established a taxonomy identifying the important factors to differentiate these methods. Four essential design factors are defined: interaction pattern (how the interaction with the DM is scheduled during the run), preference information (how the preference information is obtained from the DM), preference model (utility function, dominance relation or decision rules), and search engine (how the interesting solutions are produced, e.g. mathematical programming techniques or heuristics).

For the new method developed in this paper, the interaction with the DM is done during the run, the preference information is retrieved from pairwise comparisons, the preference model is a utility function and the search engine can be both mathematical programming techniques or heuristics. To the best of our knowlegde, the existing methods that share the same factors are: the Interactive Evolutionary Metaheuristic (IEM) [START_REF] Phelps | An interactive evolutionary metaheuristic for multiobjective combinatorial optimization[END_REF], the interactive Pareto Memetic Algorithm (IPMA) (Jaszkiewicz 2004), the Progressively Interactive Evolutionary Multi-Objective approach using Value Functions (PI-EMOVF) (Deb et al. 2010), the Necessary-preference-enhanced Evolutionary Multi-objective Optimizer (NEMO) [START_REF] Branke | Learning value functions in interactive evolutionary multiobjective optimization[END_REF], the Brain-Computer Evolutionary Multi-objective Optimization Algorithm (BC-EMOA) [START_REF] Battiti | Brain-computer evolutionary multiobjective optimization: A genetic algorithm adapting to the decision maker[END_REF] and the interactive evolutionary multi-objective algorithm with preference model called INSPM [START_REF] Pedro | INSPM: An interactive evolutionary multiobjective algorithm with preference model[END_REF]. However, all but the first two of these methods are adaptations of NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF], which is the reference algorithm for solving continuous multi-objective problems, and they have not been tested in MOCO problems. The first two methods only consider utility functions based on linear or Chebyshev aggregation. In IEM, linear programming techniques are used to learn the parameters of the utility function, whereas an internal genetic algorithm is employed in IPMA. In (Jaszkiewicz 2004), IEM and IPMA were directly compared and IPMA achieved better results on the considered MSTP instances. IPMA follows the classical steps of genetic algorithms while generating pairwise comparison queries periodically to reduce the set of possible utility functions. The frequency of preference queries is controlled by a comparison probability, which is progressively reduced during the search. Instead, our algorithm uses regret-based incremental elicitation techniques developed in the AI community to select informative preference queries and generate promising solutions during the search. As a result, our method generates at most 26 queries on existing MTSP instances with 300 cities and 7 objectives (see the numerical section) whereas IPMA needs between 30 and 60 queries on smaller instances (150 cities and 6 objectives). Moreover, we propose to apply the crossover and mutation operators on the preference parameters instead of solutions to better cover the space of admissible utility functions.

Recently it has been proposed to combine search and regret-based incremental elicitation by asking preference queries during the construction of the (near-)optimal solution [START_REF] Benabbou | Incremental weight elicitation for multiobjective state space search[END_REF]. The basic principle consists in constructing the optimal solution from optimal subsolutions using the available preference information, asking new preference information only when necessary. In this paper, we explore another way by considering nonconstructive algorithms. More precisely, we introduce an interactive genetic algorithm that uses incremental elicitation techniques to select individuals from populations. The main novelty is to make use of metaheuristics instead of constructive algorithms, which enables to tackle preferencebased combinatorial optimization problems for which no efficient exact algorithm is known. Moreover, the combination of preference elicitation and genetic algorithms has several specific advantages. First, preference queries only involve complete feasible solutions. This provides at least two advantages: 1) solutions are easier to compare, and 2) no independence assumption is required in the definition of preferences. The latter point is of special interest when preferences are represented by non-linear decision models because the cost of partial solutions is often a poor predictor of the actual cost of their extensions. Another interest of combining regret-based incremental elicitation and genetic algorithms is to produce interactive methods with performance guarantees: the proposed method 1) generates no more than a polynomial number of queries when the population size and number of generations are polynomial (as preference information is only used to discriminate between solutions within a population), and 2) runs in polynomial time for MOCO problems that can be solved (exactly or approximately) in polynomial time when preference are known.

More recently, it has been proposed to combine search and regret based incremental elicitation in a different way [START_REF] Benabbou | A general interactive approach for solving multi-objective combinatorial optimization problems with imprecise preferences[END_REF]. The idea is to identify informative queries by exploiting the extreme points of the polyhedron representing the admissible preference parameters; this new method is called IEEP for Incremental Elicitation based one Extreme Points. Our algorithm and IEEP are both general in the sense that they can be applied to any MOCO problem, provided that the scalarizing function is linear in its parameters (e.g., a weighted sum, a Choquet integral) and that a (near-)optimal solution can be efficiently determined when preferences are known. However, these two methods give different performance guarantees: IEEP is an exact exponential time method that may generate an exponential number of queries, whereas our method is a poly-time heuristic that asks no more than a polynomial number of queries. These two procedures will be compared in the numerical section.

Background and Notations

In this paper, we consider a general MOCO problem with n objective functions y i , i ∈ {1, . . . , n}, to be minimized. In this problem, any solution x ∈ X is associated with a performance vector y(x) = (y 1 (x), . . . , y n (x)) ∈ R n , where X is the feasible set in the decision space and y i (x) is the evaluation of x on the i-th objective.

Here we assume that the DM's preferences over solutions can be represented by a scalarizing function f ω that is linear in its parameters ω (e.g., weighted sums, Choquet integrals). Formally, x ∈ X is preferred to x ∈ X iff f ω (y(x)) ≤ f ω (y(x )). We also assume that parameters ω are initially not known. Instead, we are given a (possibly empty) set Θ of preference statements of type (a, b) ∈ R n × R n , meaning that the DM prefers vector a to vector b. Given Θ, we consider the set Ω Θ of all parameters ω that are compatible with Θ the set of available preference information. Formally, we have

Ω Θ = {ω : ∀(a, b) ∈ Θ, f ω (a) ≤ f ω (b)}. Since f ω is
linear in its parameters ω, we can assume that Ω Θ is a convex polyhedron without loss of generality.

To make decisions with partial preference information, we use the minimax regret criterion (MMR) defined using pairwise max regrets (PMR) and max regrets (MR) as follows:

Definition 1 For any two solutions x, x ∈ X :

P M R(x, x , Ω Θ ) = max ω∈ΩΘ {f ω (y(x)) -f ω (y(x ))} M R(x, X , Ω Θ ) = max x ∈X P M R(x, x , Ω Θ ) M M R(X , Ω Θ ) = min x∈X M R(x, X , Ω Θ )
The set arg min x∈X M R(x, X , Ω Θ ) is the set of all optimal solutions according to the minimax regret decision criterion. Recommending any of these solutions allow to minimize the worst-case loss; in particular, these solutions are necessarily optimal according to the DM's preferences when M M R(X , Ω Θ ) = 0 holds.

Note that M M R(X , Ω Θ ) can only decrease when adding new preference information in Θ. This observation has led to the following elicitation approach: reduce the parameter imprecision by iteratively asking queries to the DM until the value M M R(X , Ω Θ ) drops below a given threshold δ ≥ 0 representing the maximum allowable gap to optimality [START_REF] Boutilier | Constraint-based optimization and utility elicitation using the minimax decision criterion[END_REF]); if we set δ = 0, then we obtain the (optimal) preferred solutions at the end of the execution. This approach, sometimes referred to as regretbased incremental elicitation, was efficiently used in various decision contexts, such as multicriteria decision making and voting problems (e.g., [START_REF] Lu | Robust approximation and incremental elicitation in voting protocols[END_REF][START_REF] Benabbou | Incremental elicitation of Choquet capacities for multicriteria choice, ranking and sorting problems[END_REF]).

A Regret-Based Incremental Genetic Algorithm

To implement the regret-based incremental elicitation approach, one need to compute the value M M R(X , Ω Θ ) at every iteration step of the procedure. For MOCO problems, this may induce prohibitive computation times since it may require to compute the pairwise max regrets for all pairs of Algorithm 1 : Regret-Based Incremental Genetic Algorithm IN ↓ δ: threshold; Θ: a set of preference statements; f ω : a scalarizing function with unknown parameters ω.

OUT ↑: a solution in X .

--| Initialization of the admissible parameters:

Ω Θ ← {ω : ∀(a, b) ∈ Θ, f ω (a) ≤ f ω (b)} --| Generation of the initial population: P ← ComputeExtremePoints(Ω Θ ) --| Genetic Algorithm: for 1 to M do P ← P ∪ Crossover&Mutation(P, S, µ) while M M R(X P , Ω Θ ) > δ do --| Ask the DM to compare two solutions in X P : (x, x ) ← Query(X P ) --| Update preference information: Θ ← Θ ∪ {(y(x), y(x ))} Ω Θ ← {ω : ∀(a, b) ∈ Θ, f ω (a) ≤ f ω (b)} end while --| Move to the next population: x * ← arg min x∈X P M R(x, X P , Ω Θ ) P ← Selection(x * , P ) end for return x *
distinct solutions in X (see Definition 1). Therefore, we propose instead to combine search and regret-based incremental elicitation to reduce both computation times and number of preference queries. More precisely, we now introduce an interactive genetic algorithm that uses regret-based incremental elicitation techniques to select individuals from populations. Our algorithm, called RIGA for Regret-based Incremental Genetic Algorithm, follows the traditional scheme of genetic algorithms but differs in the following way:

• Our populations P are composed of pairs (ω, x ω ), where ω is a possible instance of the preference parameters and solution x ω is f ω -optimal (or almost). • The crossover and mutation operators are applied on parameter instances not on solutions. More precisely, RIGA takes as input a MOCO problem, a threshold δ ≥ 0, a scalarizing function f ω with unknown parameters ω and an initial set Θ of preference statements.

It proceeds as follows (see Algorithm 1):

Initial Population. To generate the initial population, we have to generate a set of possible preference parameters ω and then determine a solution x ω that is (near-)optimal for the precise scalarizing function f ω using an existing solving algorithm. These parameters could be generated uniformly at random but instead we propose to use the extreme points of polyhedron Ω Θ to better cover the feasible region.

Crossover and Mutation. As already mentioned, we perform crossovers and mutations on parameter vectors (not on solutions), and for every resulting parameter vector ω, we proceed as follows: we determine a solution x ω that is f ω -optimal (or almost) using an existing efficient solving method, and then we add the pair (ω, x ω ) in population P . To obtain a population of the desired size (say S), we create new parameter vectors by means of convex combinations of parameter vectors in P . This crossover operator is of particular interest in optimization problems with imprecise preference parameters because it only generates new admissible parameters from admissible parameters. In our experiments, we create a new parameter vector ω from two parameter vectors ω , ω in P as follows: ω = λω + (1 -λ)ω where λ ∈ (0, 1) is generated uniformly at random. Then, given a mutation rate µ, we perform Gaussian mutations on single objectives which yield very good results, but more sophisticated mutation operators would be worth investigating.

Selection. To create the new generation, we select K promising pairs from population P as follows:

• First, we determine a (near-)optimal solution in population P by means of regret-based incremental elicitation. More precisely, let X P be the set of all solutions in P . While M M R(X P , Ω Θ ) > δ, the DM is asked to compare two solutions x, x ∈ X P and the set Ω Θ of admissible parameters is updated by inserting the linear constraint

f ω (x) ≤ f ω (x ) (or f ω (x) ≥ f ω (x )
depending on her answer). Once M M R(X P , Ω Θ ) drops below δ, we select a solution that is optimal for the minimax regret criterion, i.e. a solution x * in arg min x∈X P M R(x, X P , Ω Θ ).

• Then, we compute the distance in objective space between x and x * for every pair (ω, x) in P and we select the K pairs that minimize the distance to breed the next generation. In our experiments, we use the Euclidean distance but other distances would be worth investigating.

Termination. RIGA stops after M steps and returns a solution arbitrary chosen in arg min x∈X P M R(x, X P , Ω Θ ), where P is the last generated population.

Computation Times and Number of Preference Queries.

Our algorithm has different tunable parameters: S the size of the population (generated by crossovers and mutations), K < S the number of pairs selected for the next generation and M the number of iterations steps. For any MOCO problem with existing solving methods that run in polynomial time (in the problem size), we can ensure that RIGA also runs in polynomial time and asks no more than a polynomial number of queries by carefully setting these tunable parameters. More precisely, if S the size of the population is also polynomial, then the selection step will generate at most S 2 comparison queries in order to determine x * the preferred solution in the population. Then it is sufficient that M the number of iteration steps be also polynomial to obtain the desired performance guarantees. To give an example, for the MTSP with a weighted sum, we can solve PMRoptimization problems using linear programming and apply the Lin-Kerninghan heuristic [START_REF] Lin | An effective heuristic algorithm for the travelingsalesman problem[END_REF] to efficiently solve the corresponding MOCO problem with precise preferences.

For illustration purposes, we now present an execution of our algorithm on a small MTSP instance: Example 1 Consider the instance of the MTSP depicted in Figure 1, which includes 6 vertices and n = 3 additive cost functions to be minimized. In this optimization problem, the set X of feasible solutions is the set of all Hamiltonian cycles, i.e. cycles that include every node exactly once. We assume here that the DM's preferences over Hamiltonian cycles can be represented by a weighted sum with the hidden weight ω * = (0.2, 0.1, 0.7). We now apply RIGA on this instance with δ = 0, S = 4, K = 2, M = 2. We start the procedure with an empty set of preference statements (i.e. Θ = ∅); hence Ω Θ the set of preference parameters is formally defined by:

Ω Θ = {ω = (ω 1 , ω 2 , ω 3 ) ∈ (0, 1) 3 | ω 1 + ω 2 + ω 3 = 1}.
In Figure 2, polyhedron Ω Θ is represented by the triangle ABC in the space (ω 1 , ω 2 ); ω 3 is implicitly defined by the constraint

ω 3 = 1 -ω 1 -ω 2 . ω1 ω2 • ω * • B 1 • A 0 • C 1 Figure 2: Initial set Ω Θ .
Initial population: We generate the initial population by determining one optimal solution for every extreme point of polyhedron Ω Θ . For this problem, we have P = {(ω A , x A ), (ω B , x B ), (ω C , x C )} where ω A = (0, 0, 1), ω B = (0, 1, 0), ω C = (1, 0, 0), y(x A ) = (23,34,26), y(x B ) = (20, 30, 31), and y(x C ) = (17,34,34).

First iteration step: First, we must perform crossovers and mutations to obtain a population of size S = 4; since |P | = 3, we have to generate one more pair. Assume that vector ω = (0.05, 0.45, 0.5) is generated by RIGA by applying the crossover operator on ω B = (0, 1, 0) and ω A = (0, 0, 1) and then performing a Gaussian mutation on the first objective. Then, weighted sum f ω is optimized by RIGA, resulting in the generation of solution x whose cost vector is (21,32,27) and the insertion of (ω , x ) in P . Now the selection step begins. We must ask preference queries to the DM until M M R(X P , Ω Θ ) ≤ δ = 0, where

X P = (x A , x B , x C , x ). Since M M R(X P , Ω Θ ) = 4 > 0,
we ask the DM to compare two solutions in X P , say x A and x B . Since f ω * (y(x B )) = 28.7 > f ω * (y(x A )) = 26.2, the DM answers: "solution x A is better than solution x B ". Then Θ the set of preference statements is updated by adding the pair (y(x A ), y(x B )); thus we have Θ = { ((23, 34, 26), (20, 30, 31))}. Moreover, Ω Θ the set of admissible parameters is restricted by imposing the linear constraint f ω (y(x A )) ≤ f ω (y(x B )), i.e. ω 2 ≤ 5 9 -8 9 ω 1 . Now Ω Θ is represented by the triangle ADE in Figure 3.

ω1 ω2 • ω * • A 0 1 - D • • E + 1 Figure 3: Ω Θ after 1 query.
Since we have M M R(X P , Ω Θ ) = 0.75 > 0, the DM is asked again to compare two solutions in X P , say x A and x C . Here the DM prefers solution x A to solution x C since we have f ω * (y(x C )) = 30.6 > f ω * (y(x A )) = 26.2. Therefore we must perform the following updates: we add the pair ((23, 34, 26), (20, 30, 31)) to Θ and we restrict Ω Θ by imposing the linear constraint f ω (y(x A )) ≤ f ω (y(x C )), i.e. ω 2 ≤ 1 -7

4 ω 1 (see Figure 4 where polyhedron Ω Θ is represented by ADFG).

ω1 ω2 • ω * • A 0 1 - D • • G • F + 1 Figure 4: Ω Θ after 2 queries. Now we have M M R(X P , Ω Θ ) = M R(x A , X P , Ω Θ ) = 0 ≤ δ.
We move to the next population: x * = x A and as K = 2 we need to select one more pair. The closest solution according to the Euclidean distance is x and thus we have P = {(ω A , x a ), (ω , x )}.

Second iteration step: Since |P | = 2 and S = 4, we need to generate two more pairs. Assume that we obtain (0.08, 0.4, 0.52) and (0.01, 0.1, 0.89) by applying the crossover operator on ω A and ω , and then performing a mutation on the resulting vectors. We then optimize the corresponding weighted sums and we obtain solutions x and x whose cost vectors are y(x ) = (21, 32, 27) and y(x ) = (23,34,26). Hence we have

P = ((ω A , x A ), (ω , x ), (ω , x ), (ω , x )). Since M M R(X P , Ω Θ ) = M R(x A , X P , Ω Θ ) = 0 ≤ δ, no
questions are needed at this step. This step ends by setting x * = x a and P ((ω A , x A ), (ω , x )).

Thus, after two generations, RIGA stops by returning the solution x * = x A (corresponding to cycle 1 -4 -3 -2 -5 -6 -1) which is actually the optimal solution according to the DM's preferences. Note that only two preference queries were needed to discriminate between the 60 feasible solutions (among which 10 are Pareto-optimal).

Experimental Results

In this section, we provide numerical results aiming to estimate the performance of RIGA in terms of computation time (given in seconds), number of queries and gap to optimality (expressed in terms of percentage from the optimal solution). In our experiments, we use existing Euclidean instances of the MTSP with 50 and 300 cities, and n = 3 to 7 objectives2 . The DM's Preferences. We first assume that the DM's preferences over solutions can be represented by a weighted sum f ω with imprecise weights, and we start with an empty set of preference statements; therefore we initially have

Ω Θ = {ω ∈ R n + : n i=1 ω i = 1}.
During the execution of RIGA, the answers to queries are simulated using a weighting vector ω randomly generated before running the algorithm. This weighting vector is generated using the procedure presented in (Rubinstein 1982) so as to guarantee a uniform distribution of the weights.

Query Generation Strategy. Recall that, at each iteration step of RIGA algorithm, we ask preference queries to the DM until the inequality M M R(X P , Ω Θ ) ≤ δ holds, where X P is the set of all solutions in the current population P . Here queries are generated using the Current Solution Strategy [START_REF] Boutilier | Constraint-based optimization and utility elicitation using the minimax decision criterion[END_REF]): at each step, the DM is asked to compare a solution x ∈ X P achieving the minimax regret to one of its adversary's choice (i.e. a solution in arg max y∈X P P M R(x, y, Ω Θ )).

Implementation Details. Numerical tests were performed on a Intel Core i7-8550U CPU with 16 GB of RAM, with a program written in C++. PMR-optimizations are performed by CPLEX Optimizer 3 and MTSP instances with precise weights are solved using the Lin-Kernighan heuristic from the Concorde website 4 . To generate the initial population, the library polymake 5 is used to determine the extreme points of polyhedron Ω Θ .

In this section, we will compare the following regretbased solving procedures:

• RIGA: the proposed algorithm. Due to space constraint, we only give the results obtained with M = 15 generations, S = 30 solutions per iteration, K = 5 solutions selected per iteration, and δ = 0.1, which corresponds to the best results achieved by RIGA on these instances.

• RIGA S : RIGA where genetic operators are directly applied to solutions (instead of parameter vectors). Here we use single point crossovers and Gaussian mutations on single objectives.

• IEEP δ : the exact solving method introduced in (Benabbou and Lust 2019) which returns a solution x ∈ X such that MR(x, X , Ω Θ ) ≤ δ holds at the end of the procedure.

• Two-Phase δ : a two-phase method which consists in reducing the size of the Pareto front by constructing a "wellrepresented" set (Jaszkiewicz 2018) and then applying the CSS strategy on this set until the minimax regret drops below threshold δ ≥ 0; here we generate sets of size 3000.

In the following tables, the results are averaged over 100 runs and "/" means that the timeout (45 mins) is exceeded.

Firstly, we compare the performances of RIGA and RIGA S (see Tables 12). We observe that RIGA generates very few preference queries (at most 5 against 26 for RIGA) and ends after only one minute (against three minutes for RIGA). This shows that RIGA S mainly produces solutions that are quite easy to discriminate without preference information (Pareto-dominated solutions); note that RIGA cannot produce Pareto-dominated solutions when using an exact solving method to generate new solutions during crossovers and mutations. Moreover, RIGA S induces an error that is much larger than that of RIGA (about 20% against at most 2%). Therefore, we can conclude that RIGA generates solutions that are more relevant according to the DM's maker preferences. Thus it seems more appropriate to apply genetic operators on possible parameters rather than applying them on feasible solutions for solving MOCO problems with unknown preference parameters.

Secondly, we compare RIGA to IEEP δ . We use two different tolerance thresholds: δ = 0 to obtain the optimal solution and δ = 0.02 to allow for the same error as RIGA (to be as fair as possible). In Table 3 the optimal solution is very costly both in terms of computation time (the timeout is exceeded for n > 5) and number of preference queries (56 queries are needed for n = 5). Then, when comparing Tables 1 and4, we see that IEEP Thirdly, we compare RIGA to Two-Phase δ with δ = 0.02 (see Tables 1 and5). In all the settings, RIGA outperforms Two-Phase δ in terms of computation time, number of queries and error6 . This illustrates the efficiency of interweaving search and regret-based elicitation for the determination of a near-optimal solution in combinatorial domains. Finally, we give the results obtained by RIGA when assuming that the DM's preferences are represented by an ordered weighted averaging (OWA) operator (Yager 1988):

f ω (x) = n i=1 ω i y (i) (x)
where ω = (ω 1 , . . . , ω n ) ∈ R n + and (.) is a permutation sorting the components of y(x) = (y 1 (x), . . . , y n (x)) in decreasing order. Moreover, we use decreasing weights to model preferences for well-balanced solutions (as the larger weights are associated with the worst values). Since the OWA operator is not linear in y (x) , we cannot reduce the multi-objective TSP to a single-objective TSP. Therefore, to approximately solve the problem with known weights, we use a simple local search based on the 2-opt neighborhood function (Croes 1958). In addition, to solve the problem with known weights exactly (to be able to compute the error), we use a well-known linearization of the OWA operator with decreasing weights (see [START_REF] Ogryczak | On solving linear programs with the ordered weighted averaging objective[END_REF]). Since solving the corresponding linear program is very costly in practice (more than two hours for 50 cities and n = 6), we only provide the results for the MTSP instances with 50 cities and n ranging from 3 to 6 objectives. Results are averaged over 100 runs.

In Table 6, we observe that RIGA is very efficient both in terms of computation times and number of queries. Moreover, we see that the error is below 1% for all the instances. This shows that RIGA can also be used to solve multi-objective optimization problems with complex decision models, even for problems such that there is no efficient algorithm for the determination of the optimal solution with known preference parameters.

Conclusion

In this paper, we have proposed a general method based on genetic algorithms and regret-based preference elicitation (called RIGA) for solving hard multi-objective combinatorial optimization problems with imprecise preference parameters. Our interactive method makes use of regretbased incremental elicitation techniques to generate informative preference queries during the search. Moreover, one can ensure that RIGA asks no more than a polynomial number of queries (in the problem size) by simply limiting the size of the populations and the number of generations. Similarly, one can ensure that RIGA runs in polynomial time for MOCO problems that can be (approximately) solved in polynomial time when preferences are precisely known.

We have applied the proposed method on existing instances of the MTSP, considering two different aggregation functions: the weighted sum and the ordered weighted averaging operator. We have compared the performances achieved by RIGA on these instances with that of different regret-based interactive solving methods. In particular, the provided numerical results show that our method achieves better results than IEEP (the state-of-the-art regret-based incremental solving method for the MTSP) on the bigger instances (300 cities and more than 4 criteria). Overall, RIGA was shown to be very efficient in practice: it returned very high quality solutions (the error was always less than 2%) in a few minutes with a relatively small number of queries.

Our algorithm can be applied to any MOCO problem for which preferences are represented by an aggregation function that is linear in its parameters. The next step is to extend our approach to more complex decision models (e.g., Chebyshev functions) so as to enhance descriptive and perspective possibilities.

Figure 1 :

 1 Figure 1: A MTSP instance with 6 vertices and 3 objectives.

Table 1 :

 1 RIGA for the MTSP instances with 300 cities.

	, we observe that computing
	3 https://www.ibm.com/analytics/
	cplex-optimizer
	4 http://www.math.uwaterloo.ca/tsp/concorde
	5 http://polymake.org

Table 2 :

 2 RIGA S for MTSP instances with 300 cities.

Table 3 :

 3 0.02 is better than RIGA on the smaller instances (n = 3 criteria) but RIGA outperforms IEEP 0.02 on the bigger instances; to give an example, for n = 7 criteria, RIGA is significantly faster than IEEP 0.02 (12 times faster) and both procedures need no more than 26 queries to solve the problem. IEEP 0 for MTSP instances with 300 cities.

	n time(s) queries error
	3 156.76	19.77	0
	4 445.69	34.00	0
	5 1817.68	56.30	0
	6	/	/	/
	7	/	/	/
	n time(s) queries error
	3	50.46	6.47	1.07
	4 132.02	10.47	1.16
	5 307.87	15.53	0.91
	6 802.49	19.80	1.06
	7 2336.73	25.67	0.94

Table 4 :

 4 IEEP 0.02 for MTSP instances with 300 cities.

Table 5 :

 5 Two-Phase 0.02 for MTSP instances with 300 cities.

	n time(s) queries error
	3 182.17	11.55	0.67
	4 346.44	18.38	1.16
	5 422.79	24.16	1.48
	6 542.82	31.86	2.78
	7 546.18	36.02	3.78

Table 6 :

 6 RIGA for MTSP instances with 50 cities and an OWA aggregation function.

	n time(s) queries error
	3	42.92	15.50	0.89
	4	39.26	15.80	0.65
	5	40.91	15.47	0.83
	6	44.13	15.63	0.88

A solution is called Pareto-optimal if no other solution is better on all objectives while being strictly better on at least one objective.

https://eden.dei.uc.pt/ ˜paquete/tsp/

Note that the error with TwoPhase δ may exceed 2% since it focuses on a subset of the Pareto front.

Acknowledgements

This work was supported by the Paris Ile-de-France Region.


J J Benjamin Biemond
email: benjamin.biemond@tno.nl

Romain Postoyan

W P Maurice

H Heemels

Nathan Van De Wouw

J J B Biemond

On the graphical stability of hybrid solutions with non-matching jump times

We investigate stability of a solution of a hybrid system in the sense that the graphs of solutions from nearby initial conditions remain close and tend towards the graph of the given solution. In this manner, a small continuous-time mismatch is allowed between the jump times of neighbouring solutions and the 'peaking phenomenon' is avoided. We provide conditions such that this stability notion is implied by stability with respect to a specifically designed distance-like function. Hence, stability of solutions in the graphical sense can be analysed with existing Lyapunov techniques.

Introduction

Hybrid systems feature both continuous evolution in time and discrete events, and are valuable for the modelling, analysis and control of many engineering applications, see [START_REF] Goebel | Hybrid dynamical systems: Modeling, Stability and Robustness[END_REF][START_REF] Leine | Stability and convergence of mechanical systems with unilateral constraints[END_REF] and the references therein. While the stability of stationary points and sets for hybrid systems is relatively well understood, far less is known about the stability of a given time-varying and jumping solution. The stability of time-varying solutions to hybrid systems is challenging [START_REF] Leine | Stability and convergence of mechanical systems with unilateral constraints[END_REF][START_REF] Morȃrescu | Trajectory tracking control of multiconstraint complementarity Lagrangian systems[END_REF] as two nearby solutions typically show 'peaking behaviour', i.e. they experience jumps at close, but not identical, jump times and during this timemismatch interval, the state distance between both solutions will generally not be small. To analyse stability of a hybrid solution, set-stability techniques [START_REF] Forni | Follow the bouncing ball: Global results on tracking and state estimation with impacts[END_REF], ignoring the state difference in this interval [START_REF] Morȃrescu | Trajectory tracking control of multiconstraint complementarity Lagrangian systems[END_REF], or non-Euclidean distance-like functions [START_REF] Biemond | Distance function design and Lyapunov techniques for the stability of hybrid trajectories[END_REF][START_REF] Biemond | Tracking control for hybrid systems with state-triggered jumps[END_REF] have been proposed. The approaches in [START_REF] Morȃrescu | Trajectory tracking control of multiconstraint complementarity Lagrangian systems[END_REF][START_REF] Forni | Follow the bouncing ball: Global results on tracking and state estimation with impacts[END_REF] seem to be hard to generalize, while stability with respect to non-Euclidean distance-like functions is hard to interpret. Here, we investigate the stability of a single solution (in contrast to incremental stability [START_REF] Li | Basic properties and characterizations of incremental stability prioritizing flow time for a class of hybrid systems[END_REF][START_REF] Biemond | Incremental stability of hybrid systems[END_REF]) using closeness of the graphs of solutions. We define stability of a given solution in graphical sense, which implies that the Euclidean distance between the states of this solution and nearby solutions, when compared at close continuous-time instances, tends towards zero and, in addition, the difference between the continuous times used in this comparison also tends to zero when time evolves (a similar definition is used in [START_REF] Sanfelice | A hybrid systems approach to trajectory tracking control for juggling systems[END_REF]). We prove that graphical stability is implied by stability with respect to a well-designed non-Euclidean distancelike function, cf. [START_REF] Biemond | Distance function design and Lyapunov techniques for the stability of hybrid trajectories[END_REF][START_REF] Biemond | Tracking control for hybrid systems with state-triggered jumps[END_REF]. This implication allows to prove (asymptotic) stability in graphical sense using existing Lyapunov-based techniques. For the first time, a graphical and intuitive notion of stability is provided for a given solution that applies to a large class of hybrid systems and can be analysed with existing techniques.

2 Hybrid system and stability definitions

Let (x, y) = [x T , y T ] T for (x, y) ∈ R n × R m and for a set-valued mapping F , dom F := {x ∈ R n : F (x) = ∅}. Given S 1 , S 2 ⊂ R n , S 1 + S 2 denotes {y 1 + y 2 : y 1 ∈ S 1 , y 2 ∈ S 2 } and B s denotes {x ∈ R n : x ≤ s}, with • the Euclidean norm. Let ρ B1 (x, y) := inf (u,v)∈B1 (x - u, y -v) for B 1 ⊂ R n × R n . We study hybrid systems ẋ ∈ F (x), x ∈ C; x + ∈ G(x), ∈ D, (1)
as in [START_REF] Goebel | Hybrid dynamical systems: Modeling, Stability and Robustness[END_REF] where we impose the hybrid basic conditions (cf. [START_REF] Goebel | Hybrid dynamical systems: Modeling, Stability and Robustness[END_REF]) given by: Let hybrid time domain, maximal solutions to [START_REF] Goebel | Hybrid dynamical systems: Modeling, Stability and Robustness[END_REF] and tangent cone at a point x ∈ R n to a set B ⊆ R n (denoted T B (x)) be defined as in [START_REF] Goebel | Hybrid dynamical systems: Modeling, Stability and Robustness[END_REF]. We call solution φ t-complete if sup (t,j)∈dom φ t = ∞ and bounded when φ(t, j) ≤ R for all (t, j) ∈ dom φ and some R > 0.

Exploiting the Hausdorff distance between solution graphs (see [START_REF] Sanfelice | A hybrid systems approach to trajectory tracking control for juggling systems[END_REF][START_REF] Goebel | Solutions to hybrid inclusions via set and graphical convergence with stability theory applications[END_REF]), we define graphical stability of a solution as follows:

Definition 1 A t-complete solution φ to system (1) is stable in graphical sense if the following condition holds.

For any ε > 0, there exists δ( ) > 0 such that for any maximal solution φ with φ (0, 0) -φ(0, 0) < δ( ) it holds that (i) for all (t, j) ∈ dom φ , there exists (t , j

) ∈ dom φ with |t -t | < ε such that φ (t, j) -φ(t , j ) < ε, (ii) for all (t , j ) ∈ dom φ, there exists (t, j) ∈ dom φ with |t -t | < ε such that φ (t, j) -φ(t , j ) < ε.
The solution φ is asymptotically stable in graphical sense if, in addition, there exists r > 0 such that for any ε > 0 and any maximal solution φ with φ (0, 0)φ(0, 0) < r there exists T ≥ 0 for which the following statements hold: (iii) for all (t, j) ∈ dom φ with t ≥ T , there exists (t , j ) ∈ dom φ with |t -t | < ε such that φ (t, j)φ(t , j ) < ε, (iv) for all (t , j ) ∈ dom φ with t ≥ T , there exists (t, j ) ∈ dom φ with |t -t | < ε such that φ (t, j)φ(t , j ) < ε. 2

Definition 1 prioritises continuous time over the jump counter (in [START_REF] Biemond | Incremental stability of hybrid systems[END_REF] and references therein, this is used for incremental stability). Focussing on hybrid systems that cannot exhibit consecutive jumps without flow, [5, Definition 1] yields the distance-like function ρ A (x 1 , x 2 ) with

A := (x 1 , x 2 ) ∈ (C ∪ D ∪ G(D)) 2 : x 1 = x 2 or x 2 ∈ D ∧ x 1 ∈ G(x 2 ) or x 1 ∈ D ∧ x 2 ∈ G(x 1 )
and by [8, Lemma 1] the stability definition in [START_REF] Biemond | Tracking control for hybrid systems with state-triggered jumps[END_REF] is equivalent to:

Definition 2 A t-complete solution φ to system (1) is stable with respect to ρ A if the following conditions hold. For any ε w > 0, there exists δ w (ε w ) > 0 such that for any maximal solution φ with ρ A (φ (0, 0), φ(0, 0)) < δ w (ε w ) it holds that (i) for all (t, j) ∈ dom φ , there exists (t, j ) ∈ dom φ such that ρ A (φ (t, j), φ(t, j )) < ε w , and (ii) for all (t, j ) ∈ dom φ, there exists (t, j) ∈ dom φ such that ρ A (φ (t, j), φ(t, j )) < ε w . The solution φ is asymptotically stable with respect to ρ A if, in addition, there exists r w > 0 such that for any ε w > 0 and any maximal solutions φ with ρ A (φ (0, 0), φ(0, 0))<r w there exists T w ≥0 for which it holds that (iii) for all (t, j) ∈ dom φ with t ≥ T w , there exists (t, j ) ∈ dom φ such that ρ A (φ (t, j), φ(t, j )) < ε w , and (iv) for all (t, j ) ∈ dom φ with t ≥ T w , there exists (t, j ) ∈ dom φ such that ρ A (φ (t, j), φ(t, j )) < ε w . 2

Comparison of stability concepts

The following result allows to compare both definitions.

Theorem 1 Consider system (1) with a t-complete solution φ and let the following conditions hold:

(i) G(D) ∩ D = ∅, G(D) ⊂ C and G is single-valued and proper; (ii) ∀x ∈ C ∩ D, F (x) ∩ T C (x) = ∅; (iii) ∀x ∈ C ∩ G(D), -F (x) ∩ T C (x) = ∅; (iv) either D is bounded or φ is bounded.
Then, for all ε > 0 there exists s > 0 such that for any t-complete solution φ to (1) that satisfies the conditions: (v) φ (0, 0) -φ(0, 0) < s;

(vi) for all (t, j) ∈ dom φ , there exists (t, ) ∈ dom φ such that ρ A (φ (t, j), φ(t, )) < s, it holds that for any (t, j)∈dom φ , there exists (t , j )∈ dom φ with |t-t |< ε such that φ (t, j)-φ(t , j ) < ε. 2

The existence of a solution φ verifying (v) and (vi) is guaranteed by [1, Proposition 6.14] and that the conclusion of Theorem 1 coincides with (i) of Definition 1. Condition (i) in this theorem restricts the possibility of forward and backward jumps of the hybrid system (e.g., excluding Zeno-type solutions). With the definition of the set A, this condition is essential to draw conclusions on the difference between solutions φ and φ when condition (vi) holds. Condition (ii), combined with the fact that F (x) is nonempty for x ∈ C, guarantees that solutions that are close to the jump set, will indeed jump in the near future (with a uniform bound on the jump time mismatch). In particular, it ensures C ∩ D has zero Lebesgue measure and in case D is a submanifold with C located on one side of this manifold, (ii) enforces transversal intersection of solutions with D. Condition (iii) has a similar role for solutions backward in time and also guarantees that solutions cannot enter G(D) by flow. To infer compactness results for a bounded subset of the jump set D including those points of D explored by the solution φ , condition (iv) is imposed. We now formulate our main result below.

Theorem 2 Consider a t-complete solution φ to system (1), suppose that conditions (i)-(iv) of Theorem 1 hold. If the solution φ is (asymptotically) stable with respect to ρ A (as in Definition 2), then it is (asymptotically) stable in graphical sense as in Definition 1. 2

Analysis of stability in graphical sense is facilitated by Theorem 2, which states that this stability notion is im-plied by stability of the solution with respect to a specifically constructed distance-like function. Hence, existing Lyapunov-based approaches as in [START_REF] Biemond | Tracking control for hybrid systems with state-triggered jumps[END_REF] can be used to prove asymptotic stability in graphical sense. An example where Theorem 2 is used to prove stability in graphical sense for a bouncing ball tracking problem can be found in [START_REF] Biemond | On the graphical stability of hybrid solutions with non-matching jump times: Extended paper[END_REF]. An open question is if the stability definitions in graphical sense, or in terms of ρ A , are equivalent.

A Proofs of Theorems 1 and 2 Lemma 3 Consider system (1), suppose (ii) of Theorem 1 holds and let K > 0 be given. If D is bounded then for all 2 > 0, there exists 1 > 0 such that for any tcomplete solution φ to (1) and any (t, j) ∈ dom φ such that φ(t, j)

∈ (C ∩ D) + B 1 ,
γ > 0, F > 1 such that f < F for all f ∈ F (x) and x ∈ (C ∩ D)+B γ .
Given 2 > 0, we show that 1 > 0 can be selected such that for any solution φ and any hybrid time (t, j) ∈ dom φ such that φ(t, j) ∈ (C ∩ D) + B 1 , we have

φ(s, j)-φ(t, j) ≤ |s-t| F ≤ 2 , for all (s, j) ∈ dom φ with s ∈ [t, t + min( γ 1+ F , 2 F )). (A.1)
Namely, impose 0

< 1 < min( K, γ 1+
F ) and consider a t-complete solution φ (with φ(t, j) ≤ K if D is unbounded due to the hypothesis of the lemma) and (t, j) ∈ dom φ with φ(t, j)

∈ (C ∩ D) + B 1 . Exploit- ing 1 < K and y > 2 K for all y ∈ D \ D in case D is unbounded, we find φ(t, j) ∈ (C ∩ D) + B 1 . In- troducing δt := min( γ 1+ F , 2
F ), we let s ∈ [t, t + δt] be the maximum scalar such that for all s ∈ [t, s], (s, j) ∈ dom φ and φ(s, j) ∈ {φ(t, j)} + B γ-1 hold.We deduce φ(s, j)

∈ (C ∩ D) + B γ for s ∈ [t, s] and, hence, f < F for all f ∈ F (φ(s, j)) and s ∈ [t, s]. By definition of the solution, dφ(s,j) ds ∈ F (φ(s, j)) for almost all s ∈ [t, s], such that φ(s, j) -φ(t, j) ≤ |s -t| F , for s ∈ [t, s] holds. This directly implies that for s ∈ [t, s], φ(s, j) - φ(t, j) ≤ δt F ≤ γ 1+ F F = γ -γ 1+ F < γ -1 . Hence, we find s = min{t + δt, max(s ∈ R : (s, j) ∈ dom φ)} and [t, s] coincides with {s ∈ [t, t + δt] : (s, j) ∈ dom φ}.
With s as above and |s -t| ≤ δt ≤ 2 F , (A.1) is proven. For the sake of contradiction, we now suppose: S1: there exists 2 > 0 such that for all 1 ∈ (0, γ), there exists a t-complete solution φ to system (1), with φ(0, 0)

∈ (C ∩ D) + B 1 , and φ(t , 0) / ∈ C ∩ D for all t ∈ [0, min( γ 1+ F , 2 F )].
Let 2 > 0 be as in S1. We select 1 ∈ (0, γ) such that (A.1) holds for any solution φ and hybrid time instant (t, j) ∈ dom φ for which φ(t, j)

∈ C ∩ D + B 1 and introduce δt := min( γ 1+ F , 2 F , γ-1 F ). Select 0 1 = 1 , a sequence { i 1 } i∈{0,1,...} of strictly positive scalars such that i 1 < i-1 1
for all i ∈ {1, 2, . . .}, lim i→∞ i 1 = 0, and a sequence {φ i } i∈{1,2,...} of solutions to (1) with

φ i : [0, δt ] × {0} → R n , i ∈ {1, 2, . . .}, φ i (0, 0) / ∈ C ∩ D ; φ i (0, 0) ∈ (C ∩ D) + B i 1 and φ i (t , 0) ∈ ({φ i (0, 0)} + B 2 ) \ (C ∩ D) for all t ∈ [0, δt ] (cf. S1). Introducing the compact set K = (C ∪ D) ∩ (C ∩ D + B γ ) we find with (A.1) and δt ≤ γ-0 1 F that for each solution φ i , i ∈ N, it holds that φ i (t , 0) ∈ K for t ∈ [0, δt ],
where we exploited the bound dφ(t,j) dt ≤ F that holds in this time interval. Hence, the elements of the sequence {φ i } i∈N are contained in the bounded set of absolutely continuous functions [0, δt ] × {0} → K. Consequently, there exists a convergent subsequence within {φ i } i∈N that graphically converges to a function φ : [0, δt ] × {0} → K with φ(0, 0) ∈ C ∩ D according to Theorem 5.7 in [START_REF] Goebel | Hybrid dynamical systems: Modeling, Stability and Robustness[END_REF]. This function φ is a solution to the hybrid system by sequential compactness of solutions to hybrid systems satisfying A1)-A3), cf. [1, Theorem 6.8 and Definition 6.2(a)]. However, the existence of such a solution is excluded by item (ii) of Theorem 1 as no solutions to (1) can flow on C ∩ D. Hence, S1 is contradicted and we have proven the lemma for (t, j) = (0, 0). Time-invariance of (1) concludes this proof.

2

Lemma 4 Consider a hybrid system (1) satisfying items (i), (iii) of Theorem 1 and let K ∈ R be given. If D is bounded then for all 4 > 0, there exists 3 > 0 such that for any solution φ to (1) and any (t, j) ∈ dom φ such that φ(t, j) PROOF. If j > 0, we observe that the flowing solution segment of φ to (1) is characterised by the differential inclusion x ∈ -F (x), x ∈ C as long as x ∈ G(D) and the direction of continuous time is reversed. Hence, we deduce that the statement of Lemma 4 is proven by application of Lemma 3 after replacing D with G(D). 2

∈ (C ∩ G(D)) + B 3 , j ∈ {1,
PROOF of THEOREM 1. If the jump set D is unbounded (cf. (iv)) we construct K > 0 that verifies Lemmas 3 and 4 and is such that ρ A (φ (t, j), φ(t, )) in (vi) can be written as the distance from a compact set when both (iv) and (vi) hold for some s < K. For this purpose, we select K > 0 such that both

G -1 (G(D)∩B 2K )+B K ⊂ D ∩ B K and (G(D ∩ B 2K ) + B K ) ⊆ B K hold. We define D = D ∩ B K if D is unbounded and D = D otherwise.
The set D is closed by (A1), it is thus compact. In addition, G D := G( D) is compact by locally boundedness and outer semi-continuity of G.

We introduce

A 01 = {(z 1 , z 2 ) ∈ (C ∪ D) 2 : z 2 = G(z 1 )
, z 1 ∈ D} such that the set A 01 is compact as G is locally bounded and outer semi-continuous. Introducing the symmetrical set

A 10 = {(z 1 , z 2 ) ∈ (C ∪ D) 2 : z 1 = G(z 2 )
, z 2 ∈ D}, compactness of this set follows from the symmetry. From item (i) of Theorem 1, we conclude A 01 and A 10 are not intersecting. Furthermore, the intersection of A 01 and A 00 := {(z 1 , z 2 ) ∈ (C ∪ D) 2 : z 2 = z 1 } is empty, as, for points (z 1 , z 2 ) in this intersection, z 1 = z 2 = G(z 1 ) should hold, contradicting item (i); similarly, we find A 10 ∩ A 00 = ∅. Since A 00 , A 01 , A 10 are disconnected, closed and the latter two sets compact, there exists s > 0 such that {(x, y) ∈ (C ∪D) 2 : ρ A00 (x, y) ≤ s}, {(x, y) ∈ (C ∪ D) 2 : ρ A01 (x, y) ≤ s} and {(x, y) ∈ (C ∪ D) 2 : ρ A10 (x, y) ≤ s} are mutually disconnected. As D and G D are compact and F is locally bounded by (A2), there exist positive scalars r, F such that F > 1 and f ≤ F for all f ∈ F (x) and x ∈ ( D ∪ G D ) + B r . Now, fix ε > 0 as in Theorem 1 and take 1 as in Lemma 3 with 2 = min(ε,2r) 2(1+ F ) and 3 as in Lemma 4 with 4 = min(ε,2r) 2(1+ F ) , where, if D is unbounded, K is used. Selecting s > 0 such that s < min s, ε

2 √ 2 F , r F +1 , K, 1 2 min u∈ D,w∈ G D u -w , 1 , 3 (with K = ∞ if D is bounded)
, we prove that condition (i) in Definition 1 holds. Considering any pair (φ , φ) of t-complete solutions to (1) satisfying (v),(vi) and selecting (t, j) ∈ dom φ arbitrary, we find by (vi) that there exists (t, ) ∈ dom φ such that ρ A (φ (t, j), φ(t, )) < s. Exploiting strictness of this inequality and the infimum defining ρ A , this implies that there exists (z 1 , z 2 ) ∈ A such that

(φ (t, j) -z 1 , φ(t, ) -z 2 ) < s (A.2)
holds, with z 1 , z 2 satisfying one of the following three cases that are generated by the 'or' conditions in the definition of A. We now construct (t , j ) ∈ dom φ as in the theorem. Case 1:

z 1 = z 2 ∈ C ∪ D. We directly ob- serve (z 1 , z 2 ) ∈ A 00 and select (t , j ) = (t, ). From (φ (t, j) -z 1 , φ(t, ) -z 2 ) ≥ min z∈R n (φ (t, j) -z, φ(t, ) -z) = 1 √ 2 φ (t, j) -φ(t, ) and (A.2), we conclude φ (t, j) -φ(t, j ) < √ 2s ≤ ε since s < ε 2 √ 2 F < ε √ 2 .
Hence, (t , j ) satisfies the conditions imposed in the theorem. Case 2: z 1 ∈ D, z 2 = G(z 1 ). Since φ(t, ) is close to G(D), we will apply Lemma 4 to prove that φ experienced a jump shortly before the time instant (t, ) and select time (t , j ) before this jump and show item (i) of Definition 1 holds. First, observe that z 1 ∈ D holds also in the case where D is unbounded following (A.2) and (iv). Hence, (z 1 , z 2 ) ∈ A 01 holds. To prove  > 0 in (A.2), suppose the contrary, i.e.  = 0. Let t ≤ t denote the minimum continuous time such that inf z∈ D (φ (τ, j) -z, φ(τ, 0) -G(z)) < s for τ ∈ (t , t] and [t , j] ∈ dom φ . Continuity of hybrid arcs during flow either implies (t , j -1) ∈ dom φ or there exists

z 1 ∈ D such that (φ (t , j) -z 1 , φ(t , 0) -G(z 1 )) = s. The first option implies φ(t , j) ∈ G( D) such that φ (t , j) -z ≥ min u∈ D,w∈G( D) u -w > s, con- tradicting inf z∈ D (φ (τ, j) -z, φ(τ, 0) -G(z)) < s.
Otherwise, by design of s, we find inf z∈C∪D (φ (t , j) -z, φ(t , 0) -z) > s > s (A.3) and inf z∈D (φ (t , j) -G(z), φ(t , 0) -z) > s > s. Again exploiting continuity of hybrid arcs during flow, there cannot exist a hybrid time interval [τ , t ) × {j}, with τ < t and inf z∈C∪D (φ (τ, j) -z, φ(τ, 0) -z) < s or inf z∈D (φ (τ, j) -G(z), φ(τ, 0) -z) < s for τ ∈ [τ , t ) and, since (vi) holds, the only remaining option is (t , j) = (0, 0), in which (v) contradicts (A.3). A contradiction is found in every scenario and  > 0. Since (A.2) implies φ(t, ) -G(z 1 ) ≤ s, s < K and z 1 ∈ D has been obtained above, we find φ(t, ) ∈ G(D ∩ B 2K ) + B K , such that φ(t, ) < K follows from the construction of K. As, in addition, φ(t, ) ∈ G(D)+B s holds, s < 3 and 3 is selected as in Lemma 4 with 4 = min(ε,2r) 2(1+ F ) , there exists t ∈ [t -min(ε,2r) 2(1+ F ) , t] such that φ(t , ) ∈ G(D). Similarly, we infer that

inf z∈ D (φ (τ, j) -z, φ(τ, ) -G(z)) ≤ s for τ ∈ [t , t] and t = max(t , min{t ∈ R : (t, j) ∈ dom φ }). Hence, φ (t , j) ∈ D + B s (A.4)
is obtained, which implies φ (t , j) / ∈ G(D) and t = t . From φ(t , ) ∈ G(D),  ≥ 1 and item (iii), we find (t ,  -1) ∈ dom φ. Since |t -t | ≤ min(ε,2r) 2(1+ F ) < ε, we will conclude this case and show that (t , j ), with j =  -1, satisfies φ (t, j) -φ(t , j ) < ε. We first prove

ρ A00 (φ (t , j), φ(t ,  -1)) ≤ ε 2 (A.5)
holds by considering the case of t = 0 separately, followed by the case in which t > 0. If t = 0, we use φ (t , j) / ∈ G(D) obtained above to deduce j = 0 and φ(t ,  -1) ∈ D to deduce  -1 = 0 (since G(D) ∩ D = ∅ by item (i)), such that φ (t , j) -φ(t ,  -1) = φ (0, 0) -φ(0, 0) ≤ s by (v). As ρ A00 (x, y) ≤ x -y for all x, y ∈ C ∪ D and s ≤ ε 2 , we obtain (A.5). If t > 0, using items (i) and (iii) and the inclusion (t ,  -1) ∈ dom φ, we observe that there exists a time t < t such that for τ ∈ (t , t ), the equality

{ ∈ {0, 1, . . .} : (τ, ) ∈ dom φ} = { -1} (A.6)
holds, i.e., no jumps of φ occur in the open continuoustime interval (t , t ). From (A.4) we find that (τ, j) ∈ dom φ holds for all τ ∈ (t , t ), and some t ∈ [t , t ). For τ ∈ (t , t ), (vi) implies ρ A (φ (τ, j), φ(τ,  -1)) < s for τ ∈ (t , t ). Hence, for a sequence {τ k } k∈N with τ k ∈ (t , t ) and lim k→∞ τ k = t we find lim k→∞ φ (τ k , j) ∈ D+B s , lim k→∞ φ(τ k , -1) ∈ D and lim k→∞ ρ A (φ (τ k , j), φ(τ k , -1)) ≤ s. For each sufficiently large k, ρ

A (φ (τ k , j), φ(τ k ,  -1)) = ρ A00 (φ (τ k , j), φ(τ k ,  -1)). Namely, if x 1 ∈ D + B s and x 2 ∈ D, then ρ A (x 1 , x 2 ) ≤ s implies ρ A (x 1 , x 2 ) = ρ A00 (x 1 , x 2 )
. By continuity of ρ A00 and continuity of the hybrid arcs for fixed j, , we find lim k→∞ ρ A00 (φ (τ k , j), φ(τ k ,  -1)) = ρ A00 (φ (t , j), φ(t ,  -1)) ≤ s. With ρ A00 (φ (t , j), φ(t ,  -1))

≥ 1 √ 2 φ (t , j) -φ(t ,  -1) and s < ε 2 √ 2
F , we find (A.5). From ρ A01 (φ (t, j), φ(t, )) = inf z∈D (φ (t, j)z, φ(t, ) -G(z)) < s, we find φ (t, j) ∈ D + B s and, since s ≤ r F +1 and |t -t| ≤ min(ε,2r) 2(1+ F ) ≤ r F +1 , we obtain

φ (τ, j) ≤ F for τ ∈ [t , t]. Exploiting |t -t| ≤ ε 2(1+ F )
and F 1+ F < 1, we get φ (t , j)-φ (t, j) < |t -t| F ≤ ε 2 . With (A.5) and j =  -1, we find φ (t, j) -φ(t , j ) ≤ φ (t, j) -φ (t , j) + φ (t , j) -φ(t , j ) < ε, such that (t , j ) satisfy the theorem conditions. Case 3: z 2 ∈ D, z 1 = G(z 2 ). Since φ(t, ) is close to D, we will apply Lemma 3 to prove a jump of φ will occur soon, and select (t , j ) directly after this jump. Subsequently, we conclude this case by showing that item (i) of Definition 1 holds for this hybrid time instant. For this purpose, first, we observe that z 2 ∈ D holds also in the case where D is unbounded. Namely, as (A.2) implies φ (s) -G(z 2 ) < s and s < K, we find with (iv) that G(z 2 ) < 2K. Hence, z 2 ∈ G -1 (D ∩ B 2K ) ⊆ B K holds by construction of K and (z 1 , z 2 ) ∈ A 10 is verified. Since φ(t, ) -z 2 < s < K follows from (A.2), we find φ(t, ) < K and φ(t, ) ∈ D + B s . Since s < 1 , we can apply Lemma 3 and conclude there exists a time t ∈ [t, t + min(ε,2r) 2(1+ F ) ] such that φ(t , ) ∈ D. With items (i) and (ii), we find (τ,  + 1) ∈ dom φ for τ ∈ [t , t ), with some t > t . Reasoning analogously as in Case 2, we obtain (t , j) ∈ dom φ , φ (t , j) ∈ G(D) + B s , such that φ (t , j) / ∈ D and, choosing t sufficiently small, we find (τ, j) ∈ dom φ for τ ∈ [t , t ). Taking a sequence {τ k } k∈N with τ k > t and lim k→∞ τ k = t , we find lim k→∞ ρ A00 (φ (τ k , j), φ(τ

k ,  + 1)) = ρ A00 (φ (t , j), φ(t ,  + 1)) ≤ s. As s < ε 2 √ 2 F and F ≥ 1, we find φ (t , j) -φ(t ,  + 1) < ε 2 . From (A.2) and z 2 ∈ D, we find φ (t, j) ∈ G( D) + B s and, since s < r F +1 , we obtain φ (τ, j) ≤ F for τ ∈ [t, t ], since |t -t| < min(ε,2r) 2(1+ F ) < r F +1 . We deduce φ (t , j)-φ (t, j) < |t-t | F ≤ ε 2 from |t-t | ≤ ε 2(1+ F
) . Selecting j =  + 1, we obtain φ (t, j) -φ(t , j ) ≤ φ (t, j) -φ(t , j) + φ (t , j) -φ(t , j ) ≤ ε, and the conclusion of the theorem is verified. As (t , j ) has been constructed for all three cases and arbitrary (t, j) ∈ dom φ , the theorem is proven. 2 PROOF of THEOREM 2. We exploit items (i) and (ii) of Definition 2 and Theorem 1 to conclude stability in graphical sense as defined in Definition 1, and exploit the combination of (iii), (iv) of Definition 2 and Theorem 1 to conclude (asymptotic) stability in graphical sense. Given a pair of solutions (φ , φ), Theorem 1 only provides statements for all (t, j) ∈ dom φ , (see item (i) of Definition 1). Statement (ii) of Definition 1 is attained by another application of Theorem 1 for the solution pair (φ , φ ), which we select as (φ, φ ). (Asymptotic) stability of φ with respect to ρ A will be used to show that (vi) of Theorem 1 holds for both solution pairs (φ , φ) and (φ , φ ). Consider system (1) and solution φ satisfying the conditions of Theorem 2. We select K > 0 such that either φ (t, j) < K holds for all (t, j) ∈ dom φ or

x < K for x ∈ D, cf. (iv) of Theorem 1. Exploiting also item (i) in Theorem 1 and local boundedness of G if D is unbounded, we can select a scalar K such that the combination of the inequality ρ A (φ (t, j), y) < K for any y ∈ D and item (iv) of Theorem 1 implies y < K for any (t, j) ∈ dom φ . Let φ be stable with respect to ρ A . To prove existence of a scalar δ > 0 for any ε > 0 (see Definition 1) such that (i) and (ii) in Definition 1 hold, we first fix an arbitrary ε > 0. By application of Theorem 1, we find a scalar s > 0 such that for any solution φ that satisfies (v)-(vi) of Theorem 1, item (i) of Definition 1 holds. Considering the pair (φ , φ ) = (φ, φ ) of solutions, we apply Theorem 1 and find s > 0 such that if φ satisfies (v)-(vi), with s replaced by s , for every (t , j ) ∈ dom φ, there exists a hybrid time (t, j) ∈ dom φ , with |t -t | < ε, such that φ (t, j) -φ(t , j ) < ε. We now select δ = δ w (min(s, s , K)), with δ w (•) given in Definition 2, and consider an arbitary solution φ with φ (0, 0) -φ(0, 0) < δ (cf. Definition 1) and will show that conditions (i) and (ii) in Definition 1 hold. We note that ρ A (φ (0, 0), φ(0, 0)) ≤ φ (0, 0) -φ(0, 0) < δ ≤ δ w (s) implies that item (i) of Definition 2 is verified with ε w = s. Hence, for the solution pair (φ , φ) with scalar s, condition (v) of Theorem 1 holds and, by item (i) of Definition 2, we conclude that (vi) of Theorem 1 holds. As (iv) holds by assumption, we apply Theorem 1 and conclude item (i) of Definition 1. We note that ρ A (φ (0, 0), φ(0, 0)) ≤ φ (0, 0)φ(0, 0) < δ w (s ) implies that item (ii) of Definition 2 is verified with ε w = s . Hence, aiming to apply Theorem 1 with the solution pair (φ , φ ) = (φ, φ ) and scalar s , we observe that condition (v) holds and, by item (ii) of Definition 2, we conclude that (vi) of Theorem 1 holds (also for the pair (φ , φ )). As (iv) of the same theorem holds by assumption, we can apply Theorem 1 to conclude that item (i) of Definition 1 holds for the solution pair (φ , φ ). As a direct consequence, item (ii) of Definition 1 holds for the solution pair (φ , φ). Since we have proven (i) and (ii) for this solution pair, and φ is selected arbitrarily, the solution φ is stable in graphical sense. We now show asymptotic stability. Assume that φ is asymptotically stable with respect to ρ A , let r = r w > 0 be as in Definition 2, and select ε > 0 arbitrarily. Given ε, let s be as in Theorem 1 and let s be selected as above. In addition, consider any φ with φ (0, 0)-φ(0, 0) ≤ r. Applying Lemma 3 with 2 = ε, we find some positive scalar 1 . Let s > 0 be as in the proof of Theorem 1. We now consider Definition 2 with ε w = min( s

  PROOF. Let D be given by D if it is bounded and by D ∩ B 2 K otherwise. Since D is closed by (A1), D is compact. Since F is locally bounded by (A2) of the hybrid basic assumptions, we can select

there exists t ∈ [t, t + 2 ] such that (t , j) ∈ dom φ and φ(t , j) ∈ C ∩ D. For unbounded D, such 1 and t exist if φ(t, j) ∈ B K .

2