Laurent Schnebelen

M Jourdain

M Capitaine

Yves Rocher

Vincent Piquemal

Cédric Brun

Mickaël Cozic

Jérôme Lang

Enfin

Multichoice Mckp

Knapsack

BIMOP Binary variable (or 0 ¡ 1) MOP MOCO Multiobjective Combinatorial Optimization Problem GCP Generic Coupled Problem GUCP Generic Uncoupled

Keywords:

Je remercie tout d'abord Daniel Vanderpooten d'avoir accepté d'encadrer mes travaux de thèse, et je le remercie pour ses contributions, ses suggestions, sa disponibilité, sa minutie et sa patience pendant ces quatre années.

Je remercie Jean-Michel Vanpeperstraete de m'avoir proposé un sujet de thèse aussi passionnant qu'ouvert, me laissant explorer une variété d'approches, sanctionnées par des critères de performance clairs. Je remercie également Jean-Michel pour ses nombreuses suggestions, pour nos discussions philosophiques, littéraires et cinématographiques, ainsi que pour l'honnêteté de ses conseils tant relativement à la conduite des travaux qu'à mon orientation dans un monde postdoctoral dans lequel j'ai bien peu de repères.

Je remercie Kathrin Klamroth et Vincent T'Kindt d'avoir accepté d'être rapporteurs à ma soutenance de thèse, et de se consacrer à nouveau à une lecture constructive de mon manuscrit. J'ai bénéficié lors de la présoutenance de leurs remarques bienveillantes et de conseils dont j'ai taché de tirer parti pour terminer ce document. Je remercie également Dominique Quadri et Laetitia Vermeulen-Jourdan de l'intérêt qu'elles ont témoigné à ce sujet de thèse en acceptant d'être membres du jury de soutenance.

Je remercie ensuite Hélène et Marc Kerbérénès, mes parents, pour leur soutien tout au long d'un parcours d'études long et aux revirements parfois difficiles à justifier. Je remercie Dana Mukanova d'être à mes côtés chaque jour, de m'avoir convaincu de venir étudier à Paris, de m'avoir soutenu quand j'ai connu des échecs. Je la remercie pour sa bienveillance, son intransigeance et son bon goût en toutes choses.

Je remercie Satya Tamby pour son aide et le témoignage de son expérience. Je remercie Thomas Pontoizeau pour les mêmes raisons ainsi que pour sa complicité tant intellectuelle qu'artistique, et j'ai hâte de le retrouver sur scène. Je remercie Alex Savatovsky pour la série de mardi après-midi que nous avons passés à discuter, débattre et définir des sujets de recherche parfois sérieux, nous aventurant dans des régions situées de part et d'autre de la limite de nos domaines de compétence respectifs.

Notations

R, Z, N sets of real numbers, relative integers, and positive integers p N the number or objectives, or of decision criteria R p , Z p , N p sets of p-dimensional real-valued and relative or positive integer-valued vectors I t1, ..., nu N a set of integers X the feasible set of some optimization problem Y a set of points in the objective space

x px 1 , ..., x n q an n-dimensional vector with component x i , i t1, ..., nu.

x px i | i Iq an n-dimensional vector presented as a family indexed by I N. f pxq pf 1 pxq, ..., f p pxqq a vector valued function f pXq The set of all images of elements of X by f , i.e. f pXq tfpxq R | x Xu y y I if and only if y j ¥ y I j for any j t1, ..., pu (y weakly dominates y I ) y © y I if and only if y y I and y y I (y dominates y I ) y ¡ y I if and only if y ¡ y I for any j t1, ..., pu (y strongly dominates y I ) xR f x I dR t , ©, ¡u if and only if f pxqRfpx I q.

EpX, f q the set of efficient solutions of X according to objective functions f 1 , ..., f p ẼpX, f q a subset of EpX, f q which is isomorphic to Np f pXqq NpYq the set of non-dominated points in Y 0 for some contextual p N, a vector such that 0 i 0 for all i t1, ..., pu tpxqrx i Ð as a term t where x i (either free or already valued) is assigned value a

x ¡i the subvector of x composed of all components except i, x ¡i px 1 , ...x i¡1 , x i 1 , ..., x n q x I for I N, the subvector of x composed of components with index in I. xĪ for I N, the subvector of x composed of components with index in NzI. x px ¡i , x i q x is equal up to some permutation to px ¡i , x i q R p , R

Contents

Introduction

Complex systems and their optimization

A complex system is a collection of subsystems, which are independent enough to be identified, but still linked together in significant ways. A group of agents who share the ability to perform some tasks, or share access to some resources, may be thought of as a complex system. Other examples may be a firm producing on multiple sites with partially overlapping capabilities, or a manufactured product made of several interacting components. In all of these cases, multiple subsystems are to act or be acted upon, and this action has some consequences on the state of the other subsystems.

When attempting to operate or design complex systems, a trade-off thus appears naturally between, on the one hand, bottom-up reasoning -trying to infer the optimal state of the system from the optimal states of subsystems obtained independently, and on the other hand, top-down reasoning -assuming the position of a central planner. The first type of approach is usually computationally easier than the second type, since individual subsystems are smaller than the whole system. However, by not taking their interactions fully into account, one may propose either infeasible or suboptimal solutions.

The main challenge of complex systems optimization is, therefore, to take into account the interrelations between subsystems, while never considering the problem of the whole system at once. In bottom-up approaches, it can be done implicitly by adding, to the representation of the state of a subsystem, information that accounts for some scenario of interaction with the other systems. Top-down approaches need not reduce themselves to solving the original formulation of the problem using a generic method. If interactions between subsystems can be described as discrete, and are in finite number, a top-down approach may attempt to enumerate them, and explore them in an efficient way. Finally, some approaches may appear as hybrid between bottom-up and topdown. Examples of hybrid approaches include iterative methods to achieve coordination between the subsystems: the state of the interaction is fixed, subsystem problems are solved independently, and the state of the interaction is then adapted, so as to improve the results of the next subsystems optimization.

CONTENTS

Related work and contribution

The present study focuses on a representation of complex as coupled systems, in the sense that the interaction between subsystems is modeled by coupling constraints in the mathematical representation of the optimization problem. Manipulation of these constraints and of the variables which appear in them are key to the resolution of coupled system optimization problems. The goal of these manipulations is to reduce, as much as possible, the resolution of the original, coupled problem, to the resolution of a sequence of decoupled problems that can be solved independently.

In the literature on the optimization of complex systems, this is called decomposition and coordination. In such approaches, the state of the interaction can be represented explicitly, by an upper level problem dedicated to making strategic decision coordinating subproblems, as in multilevel optimization [START_REF] Colson | An overview of bilevel optimization[END_REF], [START_REF] Migdalas | Multilevel optimization: algorithms and applications[END_REF], [START_REF] Liu | A hierarchical decomposition approach for the optimal design of a district cooling system[END_REF]). It can also be represented more implicitly by the value of a penalty term in the objective function, measuring the violation of coupling constraints, as in Lagrangian relaxation (see e.g. [START_REF] Bertsekas | Dynamic programming and optimal control[END_REF]).

We consider multiobjective integer optimization problems (MOIP), which distinguishes it from most of the classical literature on complex system optimization. Recent work has used multiobjective optimization as a device for the optimization of complex systems, see, e.g., [START_REF] Yildiz | Decomposition branching for mixed integer programming[END_REF], [START_REF] Dietz | Introducing multiobjective complex systems[END_REF]. In these approaches, each subsystem problem can be either a single, or a multiobjective problem, and the objective space of the global problem is the product of the objective spaces of the subproblems. The interaction between subsystems is represented in a higher dimensional objective space in which the objective spaces of the subproblems are embedded. The goal of such approaches is to find trade-offs between competing subsystem optimizations, which is not what we pursue here. Rather, in our case, subsystems and the global system share the same objective space, and the objective functions of the global systems are additively separable into the objective functions of the subsystems, so that subsystems are not competing.

There exists, on the one hand, literature on the application of classical decomposition and coordination methods such as Lagrangian relaxation to integer programming, e.g., [START_REF] Fisher | The lagrangian relaxation method for solving integer programming problems[END_REF], [START_REF] Lemaréchal | Lagrangian relaxation[END_REF] or [START_REF] Geoffrion | Lagrangian relaxation for integer programming[END_REF]. On the other hand, the multiobjective optimization of complex systems is a well developed topic in the literature on metaheuristics, especially genetic algorithms applied to multidisciplinary design. However, our work is aimed at exact methods, and most existing work in exact methods using decomposition in the multiobjective case is dedicated to continuous optimization (e.g., [START_REF] Nakayama | Geometric consideration of duality in vector optimization[END_REF][START_REF] Nakayama | Duality theory in vector optimization: an overview[END_REF], [START_REF] Tenhuisen | Vector optimization and generalized lagrangian duality[END_REF]) and the characterization of the efficiency of single solutions. Thus, although work from these streams of research could be built upon for the computation of supported solution to the linear relaxation of a discrete multiobjective problem, and may be adapted to approach supported points, they cannot be applied to enumerate all efficient solutions of a discrete multiobjective problem. Some topics of interest in this study are related to work by [START_REF] Gardenghi | Multiobjective Optimization for Complex Systems[END_REF] and [START_REF] Gardenghi | Algebra of efficient sets for multiobjective complex systems[END_REF] on the algebra of efficient sets, and pertain to the combination of efficient solutions to multiobjective subsystem problems, as well as the conditions under which these operations can yield efficient solutions to the global problem. We explore this topic in the particular case of CONTENTS combinatorial problems, and we study the algorithmic aspects of these operations in addition to their algebraic aspects.

Layout of the thesis

In chapter C.1, we recall basic notions of multiobjective optimization, with a particular focus on dominance filtering algorithms. We then formally introduce coupled problems, and we define the key concept of decomposition, in the special case of uncoupled subproblems. Finally, we show that in the uncoupled case, a simple implementation of decomposition already yields considerable performance increase.

In chapter 2, we study the issue of combining solutions from subproblems of a decomposable MOCO problem. We propose several ways of increasing the performance of this operation, which is fundamental to all resolution approaches subsequently developed.

In chapter 3, we propose generic multiobjective upper and lower bounds (called bound sets) to the non-dominated set of a coupled problem, and we show that they can be obtained by solving uncoupled variants of that problem. We detail several of these variants and evaluate them experimentally.

In chapter 4, we present an application problem: a multiobjective 'multi-site assignment problem denoted by REF. After having introduced the formal framework of Dynamic Programming, we show that REF admits such a resolution method, and we show how decomposition can be used to improve this initial resolution method, in a "bottom-up" approach. We also show how notions of bound sets introduced in chapter 3 adapt to the resolution of REF, using dynamic programming. We then turn to a "top-down" approach to solving REF and show that in some cases, all possible ways to decouple the problem can be enumerated and solved using decomposition. We speed this "top-down" approach by using dynamic programming to solve subproblems, by the precomputation of parts of the problem shared by all decouplings, and by applying branch & bound tricks, using an original decoupling branching scheme.

Finally, our conclusions recall the contributions made in each chapter of the thesis, as well as limitations of our results.

Chapter 1

Fundamental notions Chapter Abstract

In this chapter, we begin with recalling basic notions of multiobjective optimization, and an overview of dominance filtering algorithms. We then introduce the key concepts in this thesis: complex system optimization problem, uncoupled problems and decomposition. We define a generic coupled problem and a generic uncoupled problem, and illustrate the power of decomposition as a means to solve generic uncoupled problems.

Notions of Multiobjective optimization

A multiobjective optimization problem (MOP) is the problem of optimizing a vector-valued function f : R n Ñ R p , for n, p N. We will call an element of the domain of f a solution, and its image by f a point. In order to optimize a vector-valued function, one must specify the order relation, the minima or maxima of which are to be computed. Usually, the component-wise order over points, called the dominance relation, is used. It is only a partial order, and the set of optimal elements of chains in this partial order defines the most common solution concept of multiobjective optimization: the Pareto set, or efficient set, and its image by f : the set of non-dominated points.

In applicative contexts, the optimization of a vector-valued function is meant to solve a multicriteria decision problem. Thus, to the p dimensions of co-domain of f correspond p objective functions which are usually given independently and explicitly. It will generally be assumed that no point is optimal for all p functions, but rather, that the objectives are somewhat conflicting and cannot be optimized simultaneously.

Notations

For N t1, ..., nu a set of variable indices, let f : R n Ñ R p be a vector-valued function and X R n be the feasible set. Thus, any feasible solution x px 1 , ..., x n q X is evaluated by f pxq pf 1 pxq, ..., f p pxqq. Accordingly, the multiobjective optimization problem (MOP) associated with function f and feasible set X can be written as: max pf 1 pxq, ..., f p pxqq

x X

Feasible solutions to a MOP can be considered from the perspective of the decision space R n , as elements in X, or from the perspective of the objective space R p , as feasible points, i.e.

elements of f pXq tfpxq | x Xu. We will write Y : f pXq only when objective functions and solution sets are fixed. For any x R p , we will write y : py 1 , ..., y p q f pxq for short. In the MOP literature, comparisons in the decision and objective spaces are expressed in the following terminology (in accordance with our application cases, we use maximization as the default case for dominance):

Definition 1.1. Given y, y I R p , y is said to dominate y I , written y © y I , if y is at least as good as y I according to all objectives, and is strictly better than y I according to at least one objective.

Formally:

y © y I ô 5 dj t1, ..., pu, y j ¥ y I j hk t1, ..., pu, y k ¡ y I k CHAPTER 1. FUNDAMENTAL NOTIONS Variants of dominance in the objective space include the case where points can be equal, y y I ô dj t1, ..., pu, y j ¥ y I j and the case where strict inequalities are required on all objectives y ¡ y I ô dj t1, ..., pu, y j ¡ y I j Definition 1.2. A point y Y is non-dominated if and only if there is no y I Y such that y I © y, and is weakly non-dominated if and only if there is no y I Y such that y I ¡ y. Definition 1.3. A solution x X is efficient with respect to f if and only if there is no x I X, x I such that f px I q © f pxq, and is weakly efficient if and only if there is no x I X such that f px I q ¡ f pxq.

Thus, the usual solution concept for multiobjective optimization is the enumeration of either the set of efficient solutions defined as EpX, f q : tx X | x is efficient with respect to f u or the set of non dominated points defined as NpYq : ty Y | y is non-dominatedu in the latter case, the decision maker resorting to this solution concept is assumed to be indifferent between solutions sharing the same value vector, while in the former case two solutions with the same value vector will be returned, and discrimination between those two, if needed, should be justified by other means. Finally, let us characterize a particular subset of efficient solutions.

Definition 1.4. An efficient solution x I EpX, f q is supported if and only if there exists some λ pλ j | j t1, ..., puq R p ¡ such that x I is an optimal solution for max xX p j1

λ j f j pxq where R p ¡ ta R p | a j ¡ 0, dj t1, ..., puu. λ will be called a weight vector. A set of supported solutions can be found fairly easily by defining some policy for exploring the weight space.

In the bi-objective case, a dichotomic method introduced by [START_REF] Aneja | Bicriteria transportation problem[END_REF] makes it possible to exhaustively and efficiently enumerate the set of supported solutions (see also [START_REF] Ulungu | The two phases method: An effcient procedure to solve bi-objective combinatorial optimization problems[END_REF]). As is well known (see e.g. [START_REF] Ehrgott | Multicriteria Optimization[END_REF]), all supported solutions are efficient, but not all efficient solutions are supported. It has furthermore been observed, e.g. by [START_REF] Visée | Two-phases method and branch and bound procedures to solve the bi-objective knapsack problem[END_REF] in the case of the bi-objective knapsack problem, that the number of unsupported solutions actually grows faster than the number of supported solutions as the instance size increases.

Finally, let us define the lexicographic order relation © Lex , which is often used in the context CHAPTER 1. FUNDAMENTAL NOTIONS of multiobjective optimization, and will appear in several of the methods we develop or employ.

Definition 1.5. For y py 1 , y 2 , ...y p q and y I py I 1 , y I 2 ..., y I p q R p , y © Lex y I if and only if y 1 ¡ y I 1 , or y 1 y I 1 and py 2 , ...y p q © Lex py I 2 , ..., y I p q.

1.1.2 The specific difficulties of 0 ¡ 1 integer MO

We take special interest in problems with binary decision variables, i.e. problems which admit a formulation of the form:

max pf 1 pxq, ..., f p pxqq x X x i t0, 1u d i t1, ..., nu Such problems may be called binary integer multiobjective problems (BIMOP). We focus on this type of problems because the application we study in chapter 4 is a generalized assignment problem, in which individual assignments are best described as binary variables.

Most problems studied in the BIMOP literature, which include multiobjective combinatorial optimization problems (MOCO), are NP-hard, even when their single-objective counterparts are of lower complexity (see, e.g., [START_REF] Serafini | Some considerations about computational complexity for multi objective combinatorial problems[END_REF]). Several aspects of this difficulty can be highlighted, which are of particular interest to our work.

In a multiobjective problem with sufficiently conflicting objectives, there are not one but many solutions to be outputted, according to the non-dominated set solution concept. It has been shown (see, e.g., [START_REF] Bazgan | On the number of non-dominated points of a multicriteria optimization problem[END_REF]) that the number of efficient solutions to an integer MOP (IMOP) grows exponentially with two factors. On the one hand, with the number of dimensions of the objective space, which impacts the filtering strength of the dominance relation: the higher the number of dimensions, the lower the likelihood that two random solutions will be comparable. Thus, algorithms requiring multiple dominance tests between feasible or partial solutions may be viable in the bi-objective case. But beyond, as dominance loses its ability to discriminate between solutions, one should be mindful of the size of sets one attemps to compute the efficient subset of. On the other hand in the case of IMOP, the discreteness of the decision space is reflected in the set of feasible points: the magnitude of the coefficients determines the number of possible vector values, and thus the sheer number of different points. This in turn affects the potential number of non-dominated points. When a MOP problem admits a family of instances for which the number of non-dominated points is exponential with respect to the number of variables, the problem is said to be intractable [START_REF] Ehrgott | Multicriteria Optimization[END_REF]). In this case, no algorithm can enumerate the non-dominated set in polynomial time with respect to the size of the instance. In particular, the multiobjective versions of combinatorial optimization problems like shortest path, minimum spanning tree, assignment, knapsack, traveling salesperson are intractable.

The current literature in IMOP algorithmics may roughly be divided between decision space CHAPTER 1. FUNDAMENTAL NOTIONS methods and objective space methods. The former focus on the constructive search of solutions, trying to shorten this search by eliminating as many as possible of the candidate partial solutions being constructed using bounding reasoning, comparisons between partial solutions, and between partial solutions and already complete and feasible solutions. Branch & Bound (see the survey by [START_REF] Przybylski | Multi-objective branch and bound[END_REF] and [START_REF] Sourd | A multiobjective branch-and-bound framework: Application to the biobjective spanning tree problem[END_REF]) and Dynamic Programing (see e.g. [START_REF] Klamroth | Dynamic programming approaches to the multiple criteria knapsack problem[END_REF], [START_REF] Bazgan | Solving efficiently the 0-1 multi-objective knapsack problem[END_REF], [START_REF] Delort | A hybrid dynamic programming approach to the biobjective binary knapsack problem[END_REF]) are examples of such decision space methods, and they tend to run into the first previously described difficulty: comparisons in higher dimension spaces are rarely possible.

Objective space methods, rather than constructing and comparing partial solutions, attempt to modify the original problem so that each non-dominated point can be obtained with as few calls to a single objective solver as possible. To this end, two main tools are used: scalarizations, and the addition of constraints. Whereas weighted sums can only yield supported points, more complex scalarizations, such as Tchebychev distance to a reference point, can yield any non-dominated point. This was the basis for some early generic integer MOP methods, such as [START_REF] Eswaran | Algorithms for nonlinear integer bicriterion problems[END_REF] and [START_REF] Sayın | The multiobjective discrete optimization problem: A weighted min-max two-stage optimization approach and a bicriteria algorithm[END_REF]. To avoid the optimization of a nonlinear function such as Tchebychev distance and still deliver non-supported solutions, one can restrict the scalarized problem to subregions of the criterion space so that the desired unsupported solution becomes locally supported, or a solution to a local lexicographic optimization, or even the local optimization of one objective only. In the bi-objective case, this gives rise to the most competitive generic integer MOP methods: the e-constraint method, first introduced by [START_REF] Haimes | On a bicriterion formulation of the problems of integrated system identification and system optimization[END_REF]. For more than two objectives, recent literature is seeing the rise of generic algorithms, e.g. generalizations of the ε-constraint method by [START_REF] Laumanns | An efficient, adaptive parameter variation scheme for metaheuristics based on the epsilon-constraint method[END_REF], [START_REF] Kirlik | A new algorithm for generating all nondominated solutions of multiobjective discrete optimization problems[END_REF] and improvements on the constrained exploration of the objective space by [START_REF] Lokman | Finding all nondominated points of multi-objective integer programs[END_REF], [START_REF] Boland | A new method for optimizing a linear function over the efficient set of a multiobjective integer program[END_REF] as well as [START_REF] Tamby | Enumeration of the nondominated set of multiobjective discrete optimization problems[END_REF], based on the notion of search region introduced by [START_REF] Klamroth | On the representation of the search region in multi-objective optimization[END_REF]. The main heuristic described by [START_REF] Bentley | Fast linear expected-time algorithms for computing maxima and convex hulls[END_REF] consists in moving the candidate point y I found to dominate some input y to the front of r

Dominance filtering algorithms

Y. We denote it as MCtF for Move Candidate To Front. This is based on the intuition that the next input point should be compared first to candidates CHAPTER 1. FUNDAMENTAL NOTIONS found to dominate previous inputs. Points that are often found to dominate other points have good chances to be non-dominated, thus, comparing them early to each input should be profitable. [START_REF] Cornu | Local Search, data structures and Monte Carlo Search for Multi-Objective Combinatorial Optimization Problems[END_REF] interpreted this heuristic as also prescribing an insertion of a non-dominated input at the front of r Y, which we also consider, and denote by MItF, for Move Incumbent to Front.

The basic algorithm can also be improved upon very significantly by taking advantage of a prior ordering of the input set, as has been described by [START_REF] Kung | On finding the maxima of a set of vectors[END_REF]. Indeed, if Y is sorted according to © Lex , a point cannot dominate any of its predecessors. Thus, half of the dominace tests required for filtering can be saved, as it is only necessary to check whether an input is dominated by no candidate. In the biobjective case, a crucial observation is that, if Y is sorted lexicographically, then input solution y is dominated by an element in r Y if and only if it is dominated by the last element of r Y. We call a dominance filtering algorithm based on this principle a unidirectional dominance filtering algorithm (UF). Using the unidirectional algorithm described as Algorithm 7, MCtF can still be used if p ¥ 3. However MItF becomes irrelevant, since an input point cannot dominate a candidate point. The first type of alternative to list-based dominance filtering algorithms is filtering by insertion into a KDTree in the sense of [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF], which these authors incorporate into a two phase method of dominance filtering, which is described as Algorithm 2. We present this algorithm in details, because we will later adapt it to dominance filtering for sets of more complex objects than mere points.

Algorithm 1: Unidirectional DominanceFilter input : Y R p output: NpYq 1 S ortpY, © Lex q 2 r Y Ð ty 1 u 3 /*
Assume the points set is structured into a list, and initialize a second list of candidate points. Considering points in order, check whether each incoming point is dominated by one of the previously considered candidates. If it is not, add it to the list of candidates. It is now guaranteed that no point is dominated by its predecessors. Reverse the list. Now, by repeating the same operation, we will make sure that no remaining point is dominated by any of its former successors, and thus that no remaining point is dominated at all. Algorithm 2: Two-phase dominance filtering algorithm using KD-Trees InsertpT ree 2 , yq 20 return ContentpT ree 2 q Chen et al. ( 2012) improve this method in three main ways. First, by embedding the lists built in each of the two phases into KDTrees. A KDTree is a binary tree, in which each node has therefore two children, denoted LChild and RChild, and is associated with a point y R p , and two bounds. u y , l y denote respectively the multiobjective upperbound and lower bound to the subtree CHAPTER 1. FUNDAMENTAL NOTIONS rooted in y, meaning that for all point y I in the subtree, u y © y I © l y . As is well known, binary trees allow one to speed up the verification that a point is non-dominated, by avoiding tests using bound reasoning, and [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF] use them in a classic fashion, described by Algorithm 3. If a point is found to be non-dominated, it is inserted in the KDTree using Algorithm 4. Second, Chen et al. (2012) use a sieve point to try and perform first dominance tests which are likely to succeed. This is essentially the same trick as in move-to-front heurisitc. In this case the dominance power of a point is heuristically evaluated by its l 2 norm. For each input point y considered in the first phase, first check whether it is dominated by the sieve point, and if it has greater l 2 norm than the current sieve point, make y the new sieve point. Finally, [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF] use pruning, a common concept in search algorithms over tree-like structure: within the first phase of their algorithm, they sometimes check whether some nodes of the tree, seen as points, dominate previously inserted nodes, thus performing in advance some operation belonging to phase two by taking advantage of the implicit dominance relations induced by the tree structure. According to the authors, pruning is to be performed after k |Y| δ points have been considered in the first phase, for some δ R. In the following, we use parameter values recommended by the authors. (2017). An NDTree is a tree which is used as an archive to store non-dominated points. Leaf-nodes in the tree are associated with sets of points, and internal nodes are associated with points which are respectively upper and lower bounds to the elements associated with the leaf-nodes of their subtrees. An NDTree is generally used in procedures where individual candidate points, which may be dominated, are presented successively. Usage of NDTrees can be broken down into two main procedures.

input : Y R p output: NpYq 1 /*Assume that Y ty 1 , ...y n u. */ 2 Initialize T ree 1 3 s Ð y 1 4 k 0 5 for
The first procedure tests whether an incoming candidate is dominated by some point in the leaf-nodes. Because of the tree-like structure of the archive, an incoming point needs not be compared to every point in the archive. At each node, the incoming point may be dominated by the lower bound of the internal node, then it is known to be dominated, and it is rejected. It may dominate the upper bound of the internal node, and then it can eliminate that node and its subtree, and be accepted. If the candidate point weakly dominates the lower bound of a node, and is weakly dominated by the upper bound of a node, then one of three cases occurs. 1. If the node is a leaf node, then the candidate point is compared to points in the associated subset, the points it dominates are removed from it, and if it is found to be non-dominated within that set, it is added to it. 2. If the node is an internal node, the candidate is passed down to the child nodes. 3. If the candidate point is neither within the bounds of the node, nor dominated by it, nor dominating is, we stop the testing of the candidate in this subtree.

If the candidate was not found to be dominated during the previous procedure, a second procedure inserts it into the leaf node which it is closest to, for some notion of distance between a point and a node. Since each leaf node has a maximum storage capacity, they have to be split when the capacity is exceeded by the entry of a new point. Then, upper and lower bounds of nodes in paths of the tree which end with modified leaf nodes must be updated.

We implement this method using parameter values recommended by the authors, and within the frame of the following, straightforward algorithm: for each point in Y, insert it in the NDTree. When all input points have been inserted, "harvest" the contents of the leaf-nodes of the NDTree, the union of which is NpYq.

Choosing a dominance filtering algorithm

Throughout this document, we propose to tackle optimization problems with specialized methods which, we intend to show, perform better than their generic counterparts. In doing so, we use dominance filtering as a basic algorithmic building block. In order for the performance improvements we exhibit to represent as faithfully as possible the contributions we put forth, we need to be assured that we use the best dominance filtering method available.

When p 2, we will use a unidirectional algorithm described by [START_REF] Kung | On finding the maxima of a set of vectors[END_REF] in which the set Y to be filtered is first sorted in decreasing order of value of the first component. Because points are in this case of dimension 2, if an input point is dominated by a previously considered point, then it is dominated by the last point found to be non-dominated. When p ¥ 4, we use the two phase method together with insertion into a KDTree, both proposed by [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF]. We conducted additional preliminary experiments comparing their approach to [START_REF] Jaszkiewicz | Nd-tree-based update: a fast algorithm for the dynamic nondominance problem[END_REF], and found that [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF] performed better. When p 3, the picture is less clear and we will have to experiment with both unidirectional dominance filtering algorithms and [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF], and always consider the one that allows the best performance for the method we are challenging.

Coupled problems and decomposition

Decomposing an optimization problem means replacing it with several simpler subproblems, and then combine their optimal solutions in some way, so as to obtain the set of optimal solutions to the original problem. This notion has frequently been investigated in the context of the optimization of continuous complex systems (e.g., by [START_REF] Gardenghi | Multiobjective Optimization for Complex Systems[END_REF]; [START_REF] Ishibuchi | Multiobjective optimization for interwoven systems[END_REF]) in which the intuition justifying decomposition is that the sum of the computational efforts spent to solve the subsystem problems and to combine their solutions should be significantly less than the effort required to solve the whole problem. The question is then of whether we can guarantee the optimality of the combined solutions for the original problem, or at least their approximate optimality.

In the formal context of optimization problems, decomposition has two aspects: one related to the functions to be optimized (we will call it separability), one related to the feasible set (we will call it structural decomposability). Separability has to do with how the optimization of the subsystems is reflected in that of the global system: the global evaluation of the system is obtained by aggregating the objective functions of the subsytems. Decomposition is possible when this aggregation operator satisfies a separability property. Structural decomposability has to do with how feasible solutions for the subsystems relate to feasible solutions for the original problem: feasible solutions to the global system problem should be obtained by combining feasible solutions to the subsystem problems as straightforwardly as possible.

The combination of these two aspects allows for decomposition, i.e. the computation of the desired solution by solving subproblems and combining their solutions. A complex system is, however, generally not straightforwardly decomposable. For example, it can be the case that some constraints induce such an interdependence between subsystem solutions, that the optimum for one subproblem cannot be obtained independently from those of other subproblems, or that combining solutions obtained independently in subproblems yields an infeasible solution.

Decomposition and coordination approaches, which are classical in the single objective continuous case, are usually not applicable in the multiobjective case. Informally the reason is that these approaches critically rest on the ability to define a sequence of optimizations which converges to a unique optimal value. This can be done either by finding, from a feasible solution to the problem, a direction of improvement for the evaluation of variables involved in the interaction between subsystem, or, from an infeasible solution to the problem, by a modification of the objective function that will reduce the violation of coupling constraints. In the multiobjective case, because there is no such thing as a unique optimal value, and because for each feasible solution there would be as many directions of improvement as there are objective functions, such a converging sequence cannot be defined.

Decomposable formulations

Let us set the stage for a generic formulation of coupled optimization problems. The subproblems into which a coupled optimization problem is to be decomposed should reflect such structural features of a complex system as: the division into technical components in the case of the design of a technical device, distinct geographic areas in a planning problem with a spatial aspect, etc. Formally, the vector of decision variables should be split into subsystem-wise subvectors. The models we study here are such that all decision variables are local to some subsystems, at least in the basic formulations. This is to restrict ourselves to the case where the existence of a coupling is the consequence of the presence of a coupling constraint.

Thus let N t1, ..., nu be the set of decision variable indices. Let a partition S of N represent the set of subsystems. A subsystem s S is thus identified with the subset of decision variable indices involved in the constraints defining this subsystem. Thus, x, the vector of all decision variables, can be seen as px 1 , ..., x n q, or px s |s S q, where x s px i | i sq px s 1 , ..., x s n s q regroups the decision variables associated with subsystem s. In this framework, interactions between sub-CHAPTER 1. FUNDAMENTAL NOTIONS systems are modeled by constraints involving local variables belonging to different coupled subsystems.

Example 1.1. Let S ts 1 , s 2 u, with |s i | n i , i t1, 2u, since N s 1 s 2 and s 1 s 2 r variable vector x can be described as px s 1 1 , ..., x s 1 n 1 , x s 2 1 , ..., x s 2 n 2 q, which we will rather denote by px s 1 , x s 2 q.

Let K PpNq be a set of couplings between subsystems, represented by the set of indices of the variables it couples. We assume that any coupling constraint involves variables from at least two subsystems. In the following, each coupling k K corresponds to a unique coupling constraint indexed by k.

Definition 1.6. Two subsystems s, s I , s s I will be said to be coupled by constraint k if there exist i s k and j s I k. An optimization problem will be said to be uncoupled if and only if K r.

Note that for some k K, the associated coupling constraint may not involve all the variables of the subsystems which are coupled by k. We denote by x k the vector of decision variables involved in constraint k, and by x s k the vector of decision variables from subsystem s involved in constraint k, i.e. x s k px i | i s kq. Finally, we will denote by X s the feasible set of the suproblem associated with each s S , defined over variable subvector x s , and for each k K, X k will denote the feasible set, defined by coupling constraint k over variable subvector x k . Using these notations, we can express the coupling structure in the following formulation of a coupled system optimization problem. max f pxq pf 1 pxq, ..., f p pxqq

s.t. x k X k d k K x s X s d s S (P)
where for all j t1, ..., pu, f j pxq °sS f s j px s q. Assume that |S | m and px s |s S q px 1 , ..., x m q. A well-known example of a coupled problem is the block-angular structure displayed below, with, for all s S , f s pxq : ± is X i Ñ R p , where X i is the domain of variable x i :

max f 1 px 1 q f 2 px 2 q . . . f m px m q s.t. A 1 k x 1 k A 2 k x 2 k . . . A m k x m k ¤ b k d k K C 1 x 1 ¤ d 1 C 2 x 2 ¤ d 2 . . . . . . C m x m ¤ d m x i X i d i N CHAPTER 1. FUNDAMENTAL NOTIONS x 1 k x 1 k x 2 k x 2 k k Figure 1
.1: Graph representation of a coupled system described in Example 1.2

In the following, for some decision variable vector x of which a variable subvector is x s , we denote by x s the variable subvector defined by all components of x other than those of x s . We slightly abuse notation to do this irrespective of the actual order of components in x, when it causes no confusion. In short, when we write x px s , x sq, we mean that px s , x sq is equal to x up to a permutation, so that we can isolate the part of it that is related to s.

Example 1.2. Figure 1.1 provides a toy example which will illustrate the modifications later applied to obtain either restrictions or relaxations of the problem. Consider a graph G pV, Eq where vertex v V is a subset of decision variables and each edge e pv i , v j q represents the existence of constraints binding variables belonging to vertices v i and v j V. In the following graph, x s k for s t1, 2u denotes the vector of variables of subsystem s coupled by constraint k, while x s k denotes the vector of variables of subsystem s which are not coupled by constraint k. We only label coupling constraints, i.e. constraints which link variables from different subsystems.

Observation 1.1. If problem pPq is uncoupled, then its feasible set X is such that X ¹ sS X s (1.1)
where ±

is the cartesian product of a family of sets.

Proof. Each feasible solution x X is decomposed according to S , i.e. x px s | s S q, with S a partition of N and x s X s for each s S . If K r, any x px s | s S q with x s X s for all s S is feasible, that is x ± sS X s is equivalent to x X.

Now that we have defined what it means to be uncoupled for a problem formulation (i.e.

K r) and established the link between formulation decomposability and feasible set decomposability pi.e. X ± sS X s q, we will look at how optimization can take advantage of this property.

Uncoupled optimization

Crucial to uncoupled optimization is the fact that the objective function can be computed straightforwardly from subterms that relate to the independent subsystems. Recall that f j is the j-th CHAPTER 1. FUNDAMENTAL NOTIONS objective function for the optimization of the global system over set X. Under the assumption that f j pxq is additively separable for all j t1, ..., pu, the objective function of the subproblem associated with each s S is straightforwardly given as f s px s q pf s 1 px s q, ... f s p px s qq, with f s j : ± is X i Ñ R for j t1, ..., pu. Thus each subsystem is associated with a multiobjective problem max x s X s f s px s q.

Definition 1.7. Given j t1, ..., pu, f j : ± iN X i Ñ R p , and for all s S f s j : ± is X i Ñ R p , f j is said to be additively separable along S if, for all x px s | s S q with x s px i | i sq, f j pxq şS f s j px s q

For an uncoupled problem, assuming that objective functions f j are additively separable along S for j t1, ..., pu, we wish to investigate the relation between the global non-dominated set Np f pXqq and the non-dominated sets of all subsystem problems Np f s pX s qq, for s S . In particular, a natural question is whether the following equality is valid:

EpX, f q ? ¹ sS pEpX s , f s qq (1.2)
that is, we want to be able to obtain the efficient set of the original problem by simply combining solutions obtained by solving the subsystems problems. In Equation (C.1.2.2), decomposability of optimization is expressed in the decision space. The corresponding property, expressed in the objective space, would be:

Np f pXqq ? sS N pf s pX s qq
where °sS Y s is the Minkowski or set sum. We recall that for p N, and A, B R p , the Minkowski sum of A and B is defined as

A B : ta b | a A, b Bu
then for S t1, ..., mu, we define the Minkowski sum (MS) of a finite family pY s | s S q of finite subsets of R p as sS Y s : p...pY 1 Y 2 q ... Y m¡1 q Y m Proposition 1.1. If ± sS X s X and for all j t1, ..., pu, f j is additively separable along S , then Np f pXqq sS Np f s pX s qq Proof. Let f pxq Np f pXqq, with x px s | s S q. We prove that f s px s q Np f s pX s qq for all s S . If there is s S and x s X s such that f s px s q Np f s pX s qq, then there is x I s X s such that f s px I s q © f s px s q. Since X ± sS X s and x I s X s , x ¦ px I s , x sq is feasible. Moreover, since for 1 ¤ j ¤ p, f j pxq °sS f s j px s q, we get that f px ¦ q © f pxq, contradicting f pxq Np f pxqq

When p 1, the converse is true, i.e. we have max xX °sS f s px s q °sS max x s X s f s px s q. However, this is no longer the case when p ¥ 2, and we have °sS Np f s pX s qq Np f pXqq.

In other words the combination of non-dominated solutions with cartesian product can generate dominated solutions. For a counter-example, let S t1, 2u, X s tx s , x Is u for s S . In the following we write y Is j for f s j px I q. Let x px 1 , x 2 q, x I px I1 , x I2 q X X 1 ¢ X 2 , with the following values: y 1 p0, 1q, y 2 p1, 0q y I1 p2, 0q, y I2 p0, 2q we have Np f 1 pX 1 qq ty 1 , y I1 u and Np f 2 pX 2 qq ty 2 , y I2 u. Therefore y 1 y 2 pNpf 1 pX 1 qq pNpf 2 pX 2 qqq. However y I1 y I2 © y 1 y 2 , and so y 1 y 2 Np f pXqq. Thus, since °sS Np f s pX s qq may contain dominated points, only the following result holds:

Corollary 1.1. If X ±
sS X s and for j t1, ..., pu, f j is additively separable, then

Np f pXqq Np sS Np f s pX s qqq
In other words, filtering by dominance the Minkowski sum of non-dominated sets of the subproblems yields the solution of the global, uncoupled problem.

Finally, an important remark. In further applications, the process of combining solutions of subproblems, and then filtering the set of combinations, which we call a pooling process, intervenes at different stages of complex computations. In somes cases, we want to produce only the non-dominated points which are the result of pooling. Sometimes, even if pooling is the final stage of the computation, which can be the case for a decomposable problem, we may want to exhibit the efficient solutions associated with non-dominated points. And in some other cases, when pooling is an intermediary step in the computation, we may need to consider solutions rather than their values, in order to perform operations which demand that structural constraints be satisfied.

When the objective function is additively separable along S , the operation of combining two solutions is equivalent, in the objective space, to computing the sum of their values. In the decision space, this operation is a concatenation. For simplicity, we will consider that we combine objects which have two attributes: a solution vector (or matrix), and a value vector, which we will note, for object e, solpeq and valpeq. This way we can also make it apparent that computing the value of a combination of solutions can be done quicker than by iterating over the solution vector: we can simply compute the sum of their value attributes. Thus, we can use the same algorithm to perform the Minkowski sum or the cartesian product of two sets of points or of solutions, and the dominance algorithm we use apply to these solution CHAPTER 1. FUNDAMENTAL NOTIONS objects irrespective of whether one wants to retrieve the solution or the value of a solution, but in mathematical notations, we will still write ± and E when operating explicitly over sets of solutions, and °and N when operating explicitly over points.

The Generic Uncoupled Problem and the efficiency of decomposition

The first question we will adress empirically is that of whether decomposing an uncoupled problem, i.e. solving the subproblems independently, then pooling the solutions to the subproblems, leads to faster resolution than solving the uncoupled problem at once and as a whole. We study the case where solutions are retrieved, rather than values only.

The number of decision variables in the original problem is the sum of the number of decision variables in each subproblem. Thus, if the complexity of the resolution method applied to the whole problem is higher than linear in the number of decision variables, and if the pooling algorithm applied to the sets of efficient solutions of the subproblems is sufficiently fast, one can expect that decomposition in itself will yield a decrease in the time required to solve the problem.

We test this hypothesis by comparing the running time for solving a generic uncoupled (GUCP) MOCO problem using a generic MOCO method, with the running time for solving each subproblem using the same generic method, and pooling the obtained local non-dominated sets. Throughout this document, the generic MOCO method used is the e-constraint method. For p 2 we use our own implementation of the classic algorithm by [START_REF] Changkong | Multiobjective decision making[END_REF], and for p ¥ 3, we use the generalized method by [START_REF] Tamby | Enumeration of the nondominated set of multiobjective discrete optimization problems[END_REF]. We will refer to them indistinctly as e-constraint methods, and the algorithm we use will be clear from the number of criteria in context. GUCP is to be regarded as a particular case of a generic coupled GCP MOCO problem, the latter being formulated as follows : max f pxq p şS is

π i 1 x i , ...,
şS is

π i p x i q şS isk a i x i ¤ b k d k K is c i x i ¤ d s d s S x i t0, 1u d i N (GCP)
where for all i N, a i , c i , b k , d s R ¥ and π j i R ¥ for all j t1, ..., pu, so that f pxq : t0, 1u n Ñ N p . Thus GUCP is formulated as

CHAPTER 1. FUNDAMENTAL NOTIONS max f pxq şS is π i x i is c i x i ¤ d s d s S x i t0, 1u d i N (GUCP)
In each of trial of the next experiment, an instance of the GUCP is generated randomly according to a set of parameters. Instances are generated as follows. Each variable is assigned to a random subsystem according to a uniform law, which results in the number of variables in each subsystems being, on average, n s . Then for every i t1, ..., nu, we randomly choose π i , a i , c i , r0, 1000s, and for each s S , d s °i

S c i 2 .
We solve GUCP instances using two methods. First, we solve the whole problem using the generic method. Second, we decompose the instance into m : |S | subproblems, solve the instances using the same generic algorithm, and perform the pooling of the m non-dominated sets using a naive pooling Algorithm (NA) described by Algorithm 5, where E s denotes the set of efficient solutions of subproblem s S , and E Ý Ñ s denotes the result of successive pooling operations up to subproblem s S .

We simply concatenate solutions subproblem after subproblem, and then we filter the set of all combinations by dominance, using the dominance filtering algorithm we found to perform best for each value of p. Note that if we wanted to retrieve only the non-dominated points of the problem, rather than the efficient solutions, we would not need to perform line 7 of Algorithm 5.

Algorithm 5: Naive sequential pooling of solutions sets (NA) All experiments throughout this document are performed on a computer equipped with a 2.11 GHz processor and 32,0 GB RAM. All implementations of algorithms were done in C++, except for the algorithm by [START_REF] Tamby | Enumeration of the nondominated set of multiobjective discrete optimization problems[END_REF] which was implemented and Java an interfaced with our C++ methods. Solver IBM Ilog Cplex(TM) 12.9.0 was called to solve any integer or continuous problem, in particular within e-constraint methods.

input : pE s | s S q output: Ep ± sS E s q 1 E Ý Ñ 1 Ð E 1 2 for s t2, ...,
The results of this first experiment suffice to show that decomposition yields a significant performance improvement for any value of p, even when pooling is done in the most naive way, applying dominance filtering only after having combined the solutions of all subproblems. This performance improvement seems to scale both with the number of variables increases, and with the number of subsystems increases. The latter is explained by the fact that as we keep the number of variables fixed and further divide the problem in smaller subsystem problems, those subproblems are solved much faster, so that the time spent solving the problem with decomposition exceeds greatly the sum of the durations required to solve the subproblems. However, even if smaller subproblems have fewer efficient solutions, the time required by the dominance filtering of the set of combinations of solutions also grows very fast with the number of subproblems, as can be observed in column T.Pool.

For larger values of p, as computing time with the generic algorithm grows very fast, we had CHAPTER 1. FUNDAMENTAL NOTIONS to consider smaller values of n. As a consequence, the size of sets to be pooled are on average much smaller, so that filtering time is negligible, and its increase with |S | is not apparent.

Conclusions

After recalling basics notions and challenges in integer multiobjective optimization, a swell as defining basic tools, such as dominance filtering algorithms, we have introduced the topic of this work: the optimization of coupled problems, which are models for complex systems.

We identified the conditions under which a problem that is composed of several subproblems can be solved by solving the subproblems independently. Our work ranges beyond the case of uncoupled problems, however our strategy for solving coupled problems requires us to be able to solves uncoupled problems as fast as possible. The reason for this requirement is that, for solving a coupled problem, we will reduce it to several uncoupled problems.

In chapter 2 we will see that multiobjective upper and lower bounds of a coupled problem can be obtained by solving uncoupled variants of that problem. Approaches developed in further chapters will attempt to reduce the main problem, or parts of it, to collections of uncoupled subproblems. In each of theses cases, the same methods for pooling sets of solutions of subproblems will be applied. In the next chapter we consider the algorithmic details of the pooling operation, and investigate ways to further increase the advantage yielded by decomposition.

Chapter 2

Computing efficiently the non-dominated subset of the Minkowski set sum1

Chapter Abstract

In this chapter, we study the computation of the non-dominated subset of the Minkowski sum (or set sum) of a finite family of finite sets of multidimensional, real-valued points. At first sight, this computation may be thought of as computing all combinations of points and applying dominance filtering to the resulting set of points. However, structural properties can be exploited algorithmically to perform this task more efficiently, by avoiding the comparison of some points. We propose two main approaches of doing so. First, we show that using lexicographic ordering over each set to be added to the sum, we can speed up the sorting phase required to efficiently apply dominance filtering to elements of the set sum. Second, we introduce box-based methods, defining dominance relations between sets of combinations of points. We adress the problem of computing Np °sS Y s q, for Y s R p for each s S , where °denotes the set sum as defined in Section C.1.2.2. It is a multiobjective problem in the sense that it requires us to output the optima of a partial order relation on a set. However it is not defined as the optimization of a vector valued function of some variables, rather, each set of points in the family is given explictly. For a set A R p given explicitly, one may consider that this optimization is the mere filtering out of dominated points in A. This is true for a set that has no particular structure that we can exploit. However, the set °sS Y s is the result of the MS operation over a family of sets, which induces some structure that can be exploited algorithmically. This operation consists in pairing and summing up points from sets of the family. Dominance filterings, pairings and summations, together with other algorithmic devices, should be organized in a way that makes computation faster than first computing the MS, and then filtering. To the best of our knowledge, this is a problem that has not yet been addressed.

Considering two finite sets Y, Z R p , computing NpY Zq requires tackling efficiently two computational challenges.

First, filtering Y Z takes up to |Y Z| 2 dominance tests using a naive algorithm. However we can improve on this, either by taking inspiration from algorithms which transfer part of the filtering effort to a prior lexicographical ordering, and then perform a unidirectional filtering as in [START_REF] Kung | On finding the maxima of a set of vectors[END_REF] (an approach further denoted as UA), or by using tree-like data structures which reduce the number of tests necessary to determine whether a point is dominated or not, as in [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF]. We can also impose a partial ordering on Y ¢ Z, based on the lexicographical orderings of Y and of Z, and use it to save dominance tests when filtering Y Z. This will be developed in Section C.2.1.3.

Second, Y Z cannot be obtained in fewer than |Y| ¢ |Z| steps. However, it may not be necessary to compute all of this set sum, since we are only interested in NpY Zq Y Z. To avoid the computation of dominated sums of points, we do not test dominance between sums of points, but between sets of sums of points. These sets of sums of points can be represented by combinations of boxes, where a box delimits a subset of either Y or Z, using an upper bound and a lower bound. One comparison between sums of box bounds can eliminate many dominated sums of points. This approach will be developed in Section C.2.2.

Algorithms implementing these methods outperform the naive approach of computing the set sum of two sets, and then filtering it. Because solving NDMSP for an arbitrary family of sets can be reduced to a sequence of summing pairs of sets and filtering the results, tricks that improve the computing time NpY Zq can have significant repercussions on the computing time for a family of arbitrary length. Next, Section C.2.1.1 details a simple, but crucial improvement on the most naive approach of first computing the whole of °sS Y s and then applying dominance filtering to the result, a pair of operations we refer to as pooling the family of sets pY s | s S q.

Intermediary filtering

Consider a finite family pY s | s S q of finite sets of points in R p . The first algorithm we consider for solving NDMSP, Algorithm 6, is a sequential algorithm in the sense that it considers sets Y s to be summed one after the other. At step s S , it generates all combinations between points from Y s and points obtained by combining points from subproblems 1 to s ¡1. The latter set is denoted by Y ÝÝÑ s ¡ 1 . The most naive sequential algorithm generates the whole of Y ÝÝÑ s ¡ 1 Y s at each step s S of the computation.

Sequential Pooling

The first algorithm we consider, Algorithm 6, is a sequential algorithm in the sense that it considers sets Y s to be summed one after the other. At step s S , it generates all combinations between points from Y s and points obtained by combining points from subproblems 1 to s ¡ 1. The latter set is denoted by Y ÝÝÑ s ¡ 1 . The most naive sequential algorithm, as shown by Algorithm 6 generates the whole of Y ÝÝÑ s ¡ 1 Y s at each step s S of the computation.

However, we can avoid the computation of the whole of °sS Y s by applying dominance filtering after adding each term of the set sum. The next result guarantees that this intermediary filtering removes no non-dominated element of °sS Y s .

Proposition 2.1. For any family pY s | s S q with Y s finite for any s S and S t1, ..., nu, for all s I ¤ n,

N £ 1¤s¤s I Y s N £ N £ 1¤s¤s I ¡1 Y s Y s I
Proof. We first prove the following lemma used twice in the proof of the proposition:

Lemma 1. Let A °1¤s¤s I Y s and B Np °1¤s¤s I ¡1
Y s q Y s I . Let y AzB. Then there is some y I B such that y I © y.

Proof of Lemma 1. If y AzB, then y y j y s I , with y j °1¤s¤s I ¡1

Y s zNp °1¤s¤s I ¡1
Y s q and y s I Y s I . Then there is some y I j °1¤s¤s I ¡1 Y s such that y I j © y j . Because Y s is finite for any s S , °1¤s¤s I ¡1 Y s is also finite, so we can assume that y I j Np °1¤s¤s I ¡1 Y s q. Then there is some y I y I j y s I such that y I © y. And because y I j Np °1¤s¤s I ¡1 Y s q, we have that y I B.

Using notations of Lemma 1, let us show that NpAq NpBq. Let y NpAq, and assume, by contradiction, that y NpBq. Either y B, or y B. If y B, then there is y I B, y I © y. Because B A, we have that y I A, and then y NpAq, which contradicts our assumption. If y B, by Lemma 1, there is some y I B such that y I © y, which yields the same contradiction. Now let us show that NpBq NpAq. Let y NpBq, and assume, by contradiction, that y NpAq. Because B A, there must be some y I AzB such that y I © y. This implies, because of Lemma 1, that there is some y P B such that y P © y I . By transitivity, we get y P © y, which means that y NpBq, in contradiction with our assumption.

Algorithm 6 uses this fact to apply an intermediary dominance filtering after combining points from each new subset considered, so that Y Ý Ñ s NpY ÝÝÑ s ¡ 1 Y s q. In pseudo-code, this is controlled by the boolean variable FiltEachS tep.

The filtering strength of the dominance relation is clearly impacted by the number of dimensions of the objective space: the higher the number of dimensions, the lower the likelihood that two randomly chosen points will be comparable. When the number of criteria is low, the cost of intermediary filtering is certainly worth paying, since the filtered set may be far smaller than the original set, profitably limiting the number of combinations to consider at the next step. When the number of criteria increases, the size of the filtered set tends to get closer to that of the original set. In spite of that, we observe that intermediary filtering remains quite profitable.

Algorithm 6: Sequential Pooling algorithm without (NA) or with (IF) intermediary filtering.

input : pY s | s S q, FiltEachS tep output: Np °sS Y s q 1 Y Ý Ñ 1 Ð Y 1 2 for s t2, ..., |S |u do 3 Y Ý Ñ s Ð r 4 for y Y ÝÝÑ s ¡ 1 do 5 for y I Y s do 6 Y Ý Ñ s Ð Y Ý Ñ s ty y I u 7 if FiltEachS tep then 8 Y Ý Ñ s Ð NpY Ý Ñ s q 9 if not(FiltEachStep) then 10 Y Ý Ñ |S| Ð NpY Ý Ñ |S| q 11 return Y Ý Ñ |S |

Experimental Results

The results of the next experiment measure the effectiveness of applying dominance filtering after each step of adding a set to the set sum, and whether it is impacted by the number of objectives.

We generate families of |S | sets of size |Y s | each. Every point is drawn randomly according to a uniform distribution over a p-dimensional hypercube of width 1000. We keep each newly drawn point only if it is not dominated by, and does not dominate, any previously kept point. The reason for this requirement is that, if some point y Y s is dominated by y I Y s , then any sum of points involving point y would be dominated by the same sum where y is replaced by y I . Here, the specific computational challenge we address is the elimination of dominated combinations of points in °sS Y s , rather than that of points already dominated in some Y s , which would necessarily be eliminated more efficiently by prior filterings of sets Y s , for s S . We generate 20 instances for each parameter set, and compare the computation time of the naive sequential pooling algorithm (NA) when dominance filtering is applied only after having computed the MS entirely, and when it is applied after each addition of an individual set to the set sum (IF).

IF Time

Step s denotes dominance filtering time at step s (so that, e.g., for 3 sets, there are 2 such steps, when filtering Y 1 Y 2 and then filtering Y Ý Ñ 2 Y 3 ). For s S , IF Size Reduction

Step s denotes the reduction in size from Y Ý Ñ s to NpY Ý Ñ s q. As regards to these measurements, experimental results suggest that the computational burden rests heavily on later stages of the computation. This is apparent in the fact that step wise size reduction increases from one step to the next. This imbalance increases significantly with p: for 5 objectives and 4 sets, step 2 takes approximately 29 times longer than step 1, and step 3 approximately 9 times longer than step 2. Meanwhile, for 4 objectives and 4 sets, step 2 takes approximately 16 times longer than step 1, and step 3 approximately 5 times longer than step 2.

Pairwise Pooling

As an alternative to the sequential structure of the sequential pooling algorithm with intermediary filtering (SqIF), we investigate whether solving NDMSP by summing and filtering separately pairs of sets, and then iterating the summing and filtering on the results of this operation, could lead to better results. A similar line of inquiry was developed by in the context of a particular dynamic Step-wise increase in performance due to filtering at each step, as a function of the number of objectives, the number of sets and the number of points in each set. Average values for trials.

programming algorithm by [START_REF] Stiglmayr | On the multicriteria allocation problem[END_REF], who named this alternative pooling scheme "cascadic". The implementation of this approach is described by Algorithm 8.

We iterate through the family of sets using an iterator incremented by two at each step. This iterator defines a pair of sets, of which we compute the MS, before to filter the latter by dominance. The result is appended to the family of sets that will be subjected to the same operation at the next step of the algorithm. If, at any step of the algorithm, the current family of sets to be pooled is of odd length, then there is one set remaining after iterating through the family with the twin iterator. This remaining set is added to the next family of sets to be pooled. If the current family of sets to be pooled has only one element, it is the result. Intuitively, pairwise pooling should perform very similarly to the SqIF for |S | 2, 3, since the sequence of pooling operations will be the same in both approaches.

Algorithm 8: Pairwise Pooling input : pY s | s S q output: Np °sS Y s q 1 Q Ð tY 1 , Y 2 , ..., Y |S| u 2 i Ð 1 3 while |Q| ¡ 1 do 4 Q I Ð r 5 while i 1 |Q| do 6 Y I Ð NpY i Y i 1 q 7 Q I Ð Q I tY I u 8 i Ð i 2 9 if i |Q| then 10 Q I Ð Q I tY i u 11 Q Ð Q I 12 return Y I Q
We test this approach against the SqIF. Instances are generated in the same way as in the previous experiment. Results are provided in Table 2.3. In this experiment, variability was high, and only for the case of p 2, |S | 3 and |Y s | 100 did we find average time gain relative to SqIF which remained positive within one standard deviation. It was not, however, part of any apparent trend. In the context of their application problem, [START_REF] Stiglmayr | On the multicriteria allocation problem[END_REF] had found pairwise, or cascadic pooling to perform worse than the sequential alternative. For these reasons, we doubt that pairwise pooling can constitute a reliable alternative to sequential pooling, and we do not consider it anymore.

In this section, we have seen that when computing Np °sS Y s q for |S | ¡ 2, the enumeration of many elements of °sS Y s can be avoided by filtering after each step. But it still requires us to compute all of Y ÝÝÑ s ¡ 1 ¢ Y s at each step s, which cannot be done in less than Op|Y ÝÝÑ s ¡ 1 | ¢ |Y s |q steps. However because we are aiming at NpY ÝÝÑ s ¡ 1 Y s q Y ÝÝÑ s ¡ 1 Y s , we can hope to avoid enumerating the combinations of points that we can prove in advance will lead to dominated points at step s. We present two types of approaches to achieve that end: a unidirectional pooling method, and box-based methods.

A unidirectional method for pooling 2.3.1 Definition and proof of correction

The application of a unidirectional dominance filtering algorithm (UF) requires that the input elements y Y be presented in a dominance preserving order, and one such order is the lexicographic order © Lex . Given two sets of points Y, Z R p , in order to obtain NpY Zq, we could compute all of Y Z and sort it according to © Lex , obviously producing a sequence over Y Z, which would allow the use of a UF algorithm. However, this section presents an alternative procedure to generate such a sequence over Y Z without requiring the sorting and prior computation of Y Z.

Rather, it requires only that Y and Z be sorted in a dominance preserving order.

Assume Y and Z are both sorted lexicographically, and indexed according to their rank in each respective lexicographic order. We define a structure in which y i z j Y Z can either have two children, y i 1 z j and y i z j 1 , if i 1 ¤ |Y| and j 1 ¤ |Z| respectively, or only one child if only one of these conditions is met, or no child if none of the conditions is met. Similarly, any element y i z j has at most two parents: y i¡1 z j and y i z j¡1 , if i ¡ 1 ¡ 0 and j 1 ¡ 0 respectively.

This parenthood relation defines a lattice over Y Z which is a subset of © Lex over Y Z, as illustrated in Figure C.2. We use a table of size |Y| ¢ |Z|, denoted by T , to record, at entry pi, jq, the number of parents of y i z j Y Z which have not yet been tested for insertion at any iteration of the procedure. By definition, at the beginning of the procedure, the entry for the top of the lattice reads 0, entries for elements of the leftmost and rightmost paths in the lattice read 1, and all others read 2.

Example 2.1. Let Y tp6, 3q, p3, 5q, p2, 7q, p1, 10qu, and Z tp7, 4q, p6, 5q, p5, 6q, p2, 9qu, sorted by decreasing lexicographic order and indexed accordingly. Let T stand for the matrix representation of the parenthood relation. p8,12q p2,3q p8,11q p3,2q p8,12q p4,1q p8,14q p2,4q p5,14q p3,3q p7,13q p4,2q p7,15q p3,4q p4,16q p4,3q p6,16q p4,4q p3, 1 2 3 4 1 0 1 1 1 2 1 2 2 2 3 1 2 2 2 4 1 2 2 Let P NpY Zq denote the set NpY Zq being constructed. Let H denote the subset of elements of Y Z whose associated entry in T is 0. H is structured as a heap ordered by decreasing lexicographic order. It will be used to determine the next point to be tested for insertion, so as to ensure that the sequence of insertions is the © Lex order over Y Z.

Algorithm 9 (UPool) describes the unidirectional pooling process. The general idea is to guarantee that sum points are tested in decreasing lexicographic order. First, Y 1 is sorted. As in Algorithm 6, each iteration of the procedure adds set Y s to set sum Y ÝÝÑ s ¡ 1 . In the case of UPool, this set sum is also being kept sorted, so we do not need to sort

Y ÝÝÑ s ¡ 1 before computing NpY ÝÝÑ s ¡ 1 Y s q.
At initialization, we know that point y 1 z 1 belongs to NpY ÝÝÑ s ¡ 1 Y s q, and can be inserted into P (line 6). The next two largest elements are its children, y 1 z 2 and y 2 z 1 , which now have no untested parents, and are inserted into heap H in lexicographic order (line 7). Parenthood table T is initialized accordingly (lines 8 to 15). We continue inserting elements h which appear at the root of H into P and update H and T until H r.

We check whether h is dominated in P (line 18), using a generic procedure notDominatedInpy, Aq which determines whether point y is non-dominated by an element of set A. In practice, if p 2, it is sufficient to compare the incoming point with the last inserted point, as in [START_REF] Kung | On finding the maxima of a set of vectors[END_REF].

If p ¥ 4, previously accepted points are stored into a KDTree and we check whether the incoming point is dominated in the KDTree, adapting the approach by [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF]. If p 3, we may use either a generalized unidirectional method of dominance testing with heuristics such as move-to-front and insertion at the front of P, or we may use our adaptation of [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF].

If h is not dominated, we add it to P (line 19). We then have to update H and T , as described in Algorithm 10. We remove h from H (line 2) and for each child of h, if it exists (line 3 or 7), we decrement its entry in T (lines 4 and 8). If any of these has its two parents already tested (lines 5 or 9), we insert it into H in the correct place according to © Lex (line 6 or 10).

Algorithm 9: Unidirectional Pooling Algorithm (UPool) p10, 9qu. We also have T p1, 1q 0, and since its children are in H, T p1, 2q 0 and T p2, 1q 0. 1. p12, 8q is at the top of the heap and is tested for insertion in P. As it is not dominated by any Algorithm 10: Update heap H and tree T (update) input : h, H, T, Y, Z

input : pY s | s S q output: Np °sS Y s q 1 sortpY 1 , © Lex q 2 Y Ý Ñ 1 Ð Y 1 3 for s t2, ..., |S |u do 4 Y Ð Y ÝÝÑ s ¡ 1 ; Z Ð Y s 5 sortpZ, © Lex q 6 P Ð ty 1 z 1 u 7 H Ð ty 1 z 2 , y 2 z 1 u /*inserted according to © Lex */ 8 T p1, 1q Ð 0; T p1, 2q Ð 0; T p2, 1q Ð 0 9 for i t3, ..., |Y|u do 10 T pi, 1q Ð 1 11 for j t3, ..., |Z|u do 12 T p1, jq Ð 1 13 for i t2, ...,
1 /*Assume that h y i z j */ 2 H Ð Hzthu 3 if i 1 ¤ |Y| then 4 T pi 1, jq Ð T pi 1, jq ¡ 1 5 if T pi 1, jq 0 then 6 H Ð H ty i 1 z j u 7 if j 1 ¤ |Z| then 8 T pi, j 1q Ð T pi, j 1q ¡ 1 9 if T pi, j 1q 0 then 10 H Ð H ty i z j 1 u element of P, P Ð P tp13, 7q, p12, 8qu. We decrement entries of T associated with both children of p1, 2q, i.e. T p2, 2q Ð 2 ¡ 1 1 and T p1, 3q Ð 1 ¡ 1 0. Thus only p1,3q p11, 9q is inserted into H, while p1,2q p12, 8q is removed. 2. H t p1,3q
p11, 9q, p2,1q p10, 9qu. p11, 9q is at the top of the heap and is tested for insertion in P. It is not dominated in P, and thus P Ð tp13, 7q, p12, 8q, p11, 9qu. T p2, 3q Ð 2 ¡ 1 and T p1, 4q Ð 1 ¡ 1 0.

H t p2,1q

p10, 9q, p1,4q p8, 12qu. p10, 9q is at the top of the heap and is tested for insertion in P. It is dominated in P by p11, 9q, and thus P tp13, 7q, p12, 8q, p11, 9qu. T p3, 1q Ð 1 ¡ 1 0, and T p2, 2q 1 ¡ 1 0.

Ð tp13, 7q, p12, 8q, p11, 9q, p9, 11qu. T p4, 1q Ð 1 ¡ 1 0, T p3, 2q Ð 2 ¡ 1 1. 5. H t p2,2q p9, 10q, p4,1q p8, 14q, p1,4q p8, 12qu. p9, 10q is dominated by p9, 11q and thus P is unchanged. T p3, 2q Ð 1 ¡ 1 0, T p2, 3q Ð 1 ¡ 1 0. 6. H t p4,1q p8, 14q, p1,4q p8, 12q, p3,2q p8, 12q, p2,3q
p8, 11qu. p8, 14q is not dominated in P, thus P Ð P tp8, 14qu. p4, 2q is the only child of p4, 1q, and thus T p4, 2q Ð 2 ¡ 1 1. Next, we prove that Algorithm 9 yields the correct result, by showing that it fulfills the conditions for being a unidirectional dominance filtering algorithm.

Proposition 2.2. Let Y and Z be two subsets of points. Algorithm 9 generates NpY Zq.

Proof. Algorithm 9 is a unidirectional dominance filtering algorithm. Therefore, correctness is established if all elements of Y Z are tested by decreasing lexicographic order © Lex . For this purpose, we need to show that, p1q except for y 1 z 1 , all elements of Y Z are inserted into H at some iteration, and will thus be tested, and p2q the root of H at any iteration is lexicographically larger than all untested points.

For p1q let us assume by contradiction that g Y Z never appeared in H. This means that at least one of the parents of g has not been tested either, which means it did not appear in H. Reiterating this, we exhibit a chain of elements that never appeared in H, leading to the root element of the parenthood relation. This means that either y 1 z 2 or y 2 z 1 was never in H, which contradicts line 7 of Algorithm 9.

For p2q, let us assume by contradiction that, at some iteration of the while loop (lines 15-20), h Y Z is the root of H as in line 17, while there exists an untested element h I Y Z such that h I © Lex h. By line 7, we know that this may not occur at the first iteration. Therefore, we consider the earliest iteration where such a point h I might exist. Clearly h I is an untested point that does not belong to H, since H is sorted according to © Lex . However, there must exist h P H which is an ancestor of h I . If it were not the case, h I could never itself be brought into H, because h P would not be tested, implying that none of its descendants, among which h I , would never join H either. This would contradict p1q. H being sorted lexicographically, we have h © Lex h P . Because h P is an ancestor of h I , we have h P © Lex h I , which, by transitivity of © Lex , implies that h © Lex h I , contradicting the assumption.

We now compare IF and UPool from a computational complexity point of view. For this purpose, we consider both algorithms at a given iteration of the global for loop, when combining two sets Y and Z. Assume that |Y| m, |Z| n, with m ¥ n. IF and UPool differ on the prior sorting. IF requires the sorting of Y Z, performed in Opmn log mnq Opmn log mq steps. UPool requires the sorting of Y and of Z, performed in Opm log m n log nq Opm log mq steps, followed by the insertion of the mn elements into the intermediary heap structure H. By Proposition 2.3, we prove that the size of H is always less than n. Therefore, this second operation requires, in total, Opmn log nq steps.

Proposition 2.3. Given two sets Y and Z to be combined by Algorithm 9, the heap H used in this algorithm never contains more than mint|Y|, |Z|u elements.

Proof. Let Y ty 1 , ..., y m u and Z tz 1 , ..., z n u, with m ¥ n. Assume by contradiction that |H| ¡ n. Then, H must contain at least two elements of the form y i z j and y k z j . Assume, wlog, that i k, meaning that y i z j is an ancestor of y k z j in the parenthood lattice. An element is introduced in H when its entry in T is decremented to 0. The entry of an element in T is decremented only when both its father nodes, and thus all its ancestor nodes, have been removed from H. Thus, if y k z j H, then y i z j H, contradicting the assumption.

Both algorithms require filtering the mn elements of Y Z, which involves mn operations in the biobjective case, and at least Opmn log mnq Opmn log mq operations for more than two objectives. In the latter case, this term will in all likelihood dominate the one associated with sorting. Thus we cannot guarantee that UPool will perform better than IF, but we showed that the structure of NDMSP can be taken advantage of to shorten significantly the sorting phase. In the biobjective case, sorting corresponds to the dominant term, with an overall complexity of Opmn log mq for IF, versus Opmn log nq for UPool. Therefore, we can expect UPool to perform better than IF in that case.

Experimental Results

We generate families pY s | s S q in the same way as in the experiment about intermediary filtering (cf. Section C.2.1.2), and we measure time gain (TG) yielded by UPool relative to IF, using the dominance filtering algorithm which allows IF to perform best. It can be either UF with move-to-front heuristic and insertion of new non-dominated elements at the front, denoted by UF, or the two-phase algorithm by [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF], denoted by CH.

For p 2 and p 3, we found that IF performs best using UF for intermediary filtering (cf. bold-titled column), and thus we consider TG when testing and inserting elements into P using UF as well. We find significant TG, up to 81, 65 % for p 2 with a monotonous increase as |Y s | increases, and low variability. For p 3, we find the same monotonous increase, but only 49.33% TG at best. In both cases we observe a very slight increase of TG as |S | increases, suggesting the performance of the algorithm scales well with the size of the sets being handled. These results are reported in For p 4 and p 5, IF benefits more from dominance filtering with CH, and we thus measure the performance improvement of UPool by considering column TG{CH. In this case, we found no positive TG, much higher variability and no monotonous trend in TG. Results using UF within UPool performed even more poorly as compared to IF using CH, and are not reported here. Thus, for these higher values of p, we cannot advise using UPool. To summarize, for p ¤ 3, UPoolUF is the clear winner, whereas for p ¥ 4, IFCH should be prefered.

Box-based methods

Box-based methods for speeding up the resolution of NDMSP rely on dominance relations involving combinations of bounded boxes. We present two dominance relations involving combinations of boxes. The first one is a relation between a point and a combination of boxes, and thus will be called point-to-box dominance (with a slight abuse since we compare the point to a combination of boxes, and not to a box). The second one is a relation between combinations of boxes, and will accordingly be called box-to-box dominance.

Box-based dominance relations

Given some set Y R p , let a box B be defined by a subset of Y, an upper bound u such that for all y in the box, u y and a lower bound l such that for all y in the box, y l. The tightest value for u is the ideal point of B and the tightest value for l is the anti-ideal point of B. For simplicity, we will write y B to say that y is contained in the subset of box

j Z . ȳ point-to-box dominates pB j Y , B j Z q if and only if ȳ u i Y u j Z . Observation 2.1. If ȳ point-to-box dominates pB i Y , B j Z q, then for all y B i Y B j Z , ȳ y (see Figure 2.3). Therefore, if ȳ Y Z, then no point of B i Y B j Z belongs to NpY Zq.
We turn to the box-to-box relation. If combining the upper bounds of two boxes yields a point that is dominated by the combinations of lower bounds of two others boxes, then all combinations of points from the latter two boxes dominate all combinations of points from the former two boxes. Formally, we have:

Definition 2.2. Let B i Y , B j Y Y and B k Z , B l
Z Z, with associated lower and upper bounds.

pB i Y , B k Z q box-to-box dominates pB j Y , B l Z q if and only if l i Y l k Z u j Y u l Z Observation 2.2. If pB i Y , B k Z q box-to-box dominates pB j Y , B l Z q, then for all y B i Y B k Z , for all y I B j Y B l Z , y y I (see Figure C.4) B 1 Y B 1 Z u 1 Y u 1 Z u 1 Y u 1 Z ȳ Figure 2.3: Bi-objective illustration of point-to-box dominance. B 1 Y B 2 Y B 1 Z B 2 Z l 2 Y l 1 Z u 1 Y u 2 Z u 1 Y u 2 Z l 2 Y l 1 Z Figure 2.4: Bi-objective illustration of box-to-box dominance.
For efficiency, it may be desirable that the boxes be disjoint, so that dominated combinations of points are eliminated only once.

Algorithms for box-based dominance filtering

In a sequential approach to pooling, as described by Algorithm 11 at each step s S , we apply two boxing operations: on the one hand on Y ÝÝÑ s ¡ 1 , the result of the previous steps, and on the other hand on Y s , the s-th set that we want to add to the MS under construction. After removing potential combinations of points by applying either one of the box dominance relations, we apply regular dominance filtering over sums of points associated with the remaining combinations.

First, a boxing algorithm is called to determine how to group the contents of Y ÝÝÑ s ¡ 1 and of Y s into two families of boxes, using a parameter ε as explained in the next subsection. Then we

Algorithm 11: PoolNextSetBox input : Y ÝÝÑ s ¡ 1 , Y s , BtB, PtB, n λ , ε output: NpY ÝÝÑ s ¡ 1 Y s q 1 Pooled Ð r 2 pB 1 ÝÝÑ s ¡ 1 , ..., B n ÝÝÑ s ¡ 1 q Ð BoxingpY ÝÝÑ s ¡ 1 , εq 3 pB 1 s , ..., B m s q Ð BoxingpY s , εq 4 BoxCombinations Ð tB 1 ÝÝÑ s ¡ 1 , .., B n ÝÝÑ s ¡ 1 u ¢ tB 1 s , ...,

Point-to-box dominance

To filter according to the point-to-box dominance relation, we need a pool P Y ÝÝÑ s ¡ 1 Y s of point combinations. As described in Algorithm 12, we propose to use supported points, by optimizing positively weighted sums for each sets, using weight set Λ R p ¡0 generated according to parameter n λ , to produce supported points. Points of the resulting pool will be compared to the combinations of upper bounds of boxes in Algorithm 13.

The points whose values maximize some positive weight sum are of course non-dominated in the sets in which they are found. Combining two points which are non-dominated in their respective sets Y and Z can yield a dominated point in the set sum Y Z. In this case, however, we use the fact that for any sets Y, Z and λ Λ,

maxtλpy zq | y z Y Zu maxtλy | y Yu maxtλz | z Zu
to guarantee that these points are not only feasible but also non-dominated, i.e. they belong to NpY ÝÝÑ s ¡ 1 Y s q and can be immediately included in the set sum under construction.

A key step in producing supported points is generating the set of weights defining the scalarizations for which supported points are optimal. A good set of weights would produce points which evenly cover Y ÝÝÑ s ¡ 1 Y s , and do not lead to too many redundant points. Because generating those supported points takes place within a box-based dominance filtering algorithm which must be quick, we must also bound the number of weight vectors we want to consider.

In the bi-objective case, we can easily generate the set of supported points exhaustively using the dichotomic method by [START_REF] Aneja | Bicriteria transportation problem[END_REF], and an ε-dominance based threshold to stop the exploration of a weight-space interval and thus control for the desired number of weights vectors. Beyond the bi-objective case, there is no systematic and scalable method enumerating the whole set of supported points. Following [START_REF] Borges | A study of global convexity for a multiple objective travelling salesman problem[END_REF], we could generate sets of maximally dispersed weights, so as to hopefully reduce the number of redundant points. However, experimentally, we found that this approach did not yield better performance than generating random weight vectors.

The number of weight vectors to be generated is denoted by n λ , which we tuned manually. We produce supported points according to the computed set of weights, and we remove duplicates, i.e. points that were supported by several different weight vectors. Then for each box combination, we iterate over these supported points and test for point-to-box dominance between the supported point and the box combination, as described by Algorithm 13. For filtering according to the box-to-box dominance relation, we adapt the unidirectional dominance filtering algorithm (UF) to the case of comparing box combinations. Correction of a UF requires the input to be ordered in a dominance preserving order. It is easily seen that the decreasing lexicographic order over sums of upper bounds of combinations of boxes is a dominance preserving order for box-to-box dominance:

Algorithm 12: getSupported input : Y ÝÝÑ s ¡ 1 , Y s , Λ output: S upported 1 S upported Ð r 2 for λ Λ do 3 y 1 Ð maxtxλ, yy | y Y ÝÝÑ s ¡ 1 u 4 y 2 Ð maxtxλ, yy | y Y s u 5 y ¦ Ð y 1 y 2 6 if y ¦ S
Observation 2.3. Let B 1 Y , B 2 Y Y, and B 1 Z , B 2 Z Z. If u 1 Y u 1 Z © Lex u 2 Y u 2 Z , then l 2 Y l 2 Z u 1 Y u 1 Z
Algorithm 13 describes the pseudo-code for the box-to-box filtering procedure. © u Lex denotes the decreasing lexicographic order over the sums of upper bounds of combinations of boxes, i.e.

for pB 1

Y , B 1 Z q, and pB 2 Y , B 2 Z q two pairs of boxes, by definition, pB 1 Y , B 1 Z q © u Lex pB 2 Y , B 2 Z q if and only if u 1 Y u 1 Z © Lex u 2 Y u 2 Z .
If p 2, we use UF with a unique dominance comparison with the last combination of boxes found to be non-dominated. For p ¥ 3, we adapt the approach by [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF] of insertion into a KDTree to check whether an input box combination is dominated or dominates previously considered box combinations. We will denote generically by CheckDtdInpy, Aq an algorithm which determines whether some box combination y is dominated by an element of set A.

Once box-to-box filtering has been performed, we are left with a set of non-dominated pairs representing combinations of boxes. We iterate over these pairs, and for each pair of boxes, we compute the MS of their contents, and append it to Pooled. Finally we apply a dominance filtering algorithm to Pooled before returning it.

An algorithm for creating a family of boxes

A variety of algorithms can be used to produce a family of boxes from an input set of points. As preliminary work, we have compared five approaches to achieve that end. The first one consists in a quantization of each set: for each dimension, a desired coarseness is defined, and the orthogonal projection of the set onto this dimension is divided into intervals of equal length. Each point then belongs to a unique box which is defined by an interval on each dimension. The second one consists in defining a desired number of boxes, and using either ascendant hierarchical or k-means clustering to define k boxes. The third one consists in defining a minimum and a maximum number of points per box, and inserting points into their closest box until the closest box is full. When a box is full, we split it according to a rule that ensures that the new boxes contain the minimum number of points, and makes the two new boxes as distinct as possible. The fourth approach consists in defining a threshold distance, considering each point in the input set in turn, and adding it to the closest already defined box which lies at a distance smaller than the threshold. If such a box does not exist, a new box is created for the current point and the process continues. A fifth approach would be similar to the fourth one, using not a distance threshold but ε-dominance to

Y , B 1 Z q, ..., pB m Y , B n Z qu */ 8 if l m Y l n Z © u i Y u j Z then 9 Dtd Ð T rue 10 else 11 Dtd Ð CheckDtdInppB i Y , B j Z q, NDBoxCombq 12 if notpDtdq then 13 NDBoxComb Ð NDBoxComb pB i Y , B j Z q
14 return NDBoxComb determine whether two points should belong to the same box.

We will only describe in details, and provide experimental results for the fifth approach to creating boxes, based on ε-dominance, and we begin with defining the notion.

Definition 2.3. For all y, y I Y, y is said to ε-dominate y I , written y © ε y I if and only if y j p1 εq ¥ y I j dj t1, ..., pu Moreover, y is said to be ε-indifferent to y I , written y ε y I if and only if y © ε y I and y I © ε y

The boxing algorithm operates as follows. Consider a first point y Y, and create a box containing y. Then for each remaining point y I Y, check whether y ε y I . If it is the case, add y I to the box associated with y, and remove y I from Y. If not, do nothing. When all of Yztyu has been examined, remove y from Y and repeat the process with a new head of list. This procedure is formally described as Algorithm 15.

Experimental Results

In the next experiment, we measure the time gain (TG) yielded by the use of box-based methods in solving NDMSP, and the average number of points eliminated by each box ¡ based test (av. We tuned ε so as to get the best possible results. Preliminary experiments suggest that decreasing the value of ε as |Y| increases allows better results. However, as we were not able to analyze this dependency compellingly, we resorted to the use a single value of ε for each value of p. In this case, for box-to-box dominance we used ε 0.001 for p 2, ε 0.02 for p 3, ε 0.06 for p 4 and ε 0.1 for p 5. For point-to-box dominance, we used ε 0.015 for p 2 and ε 0.03 for p 3. This dominance relation requires us to choose a number of weight vectors to compute. Intuitively, it seems that n λ should vary with the number of "feasible points" at one step of the computation, i.e. with |Y ÝÝÑ s ¡ 1 ¢ Y s |. However we were not able to find a relationship between |Y ÝÝÑ s ¡ 1 ¢ Y s | and n λ for which performance was significantly better than when using a fixed value of n λ .

Results for box-to-box dominance, in Tables 2.6 and 2.7, suggest that TG does increase as |Y s | increases. For p 2, we found maximum TGs ranging from about 40.51 to about 60.74%, suggesting that box-to-box dominance is effective at speeding up the resolution of NDMSP. However, an increase in |S | appears to "flatten" the TG profile, and limit performance especially for higher values of |Y s |. This is most likely due to the use of a single parameter for boxing, which becomes unfit to larger set sizes. For p 3, we find more moderate TG. We observe, however, a more pronounced increase of performance as |Y s | increases, which suggests that scaling is probably still good, and thus the approach viable, for this number of criteria. For p 4 however, we find small TG, high variability, and no clear trend as S or |Y s | increases, thus cannot claim that the method is viable in this case. For p 5, we found no positive time gain.

Note that for p 2 and p 3 ( Average values for 20 trials.

we found that it is also supplanted by box-to-box (column TG{UPCH). However, box-to-box itself performs very similarly to IF or worse, so that IF should still be considered the most reliable approach for p 4.

Turning to results for point-to-box dominance, in Tables 2.8 and 2.9, when p 2, we found the best results when successful box-based-dominance tests tend to eliminate many points, suggesting that the optimal box size is larger. We found point-to-box dominance filtering to yield significant TG, from about 49.1 up to about 67.42%. Performance did not appear to vary significantly with increases in |Y s | or in |S |. It remains to be investigated whether this is due to tuning issues or a need for adapting the parameter values with increasing |Y s | and |S |. For p 3, we found slightly positive TG, although variability remained too high to warrant that this TG will always be positive.

For p 4, TG was always negative is not reported.

Although one may expect the optimal number of boxes to be closely dependent on the number of supported points generated, we found that for p 2, it was nearly always beneficial to generate the complete set of supported points using the dichotomic method. For p ¥ 3, we were not able to find a very stable relationship between n λ and TG. We tuned n λ by hand and measured the ratio of unique supported points to the number of generated weight vectors, and we found that about 20% of successful generation (column Sup. Ef ) of supported points was associated with highest TG. Smaller n λ will yield higher Sup. In all reported cases, box-based dominance filtering using point-to-box dominance performed worse than UPool (column TG{UPUF). Thus as a conclusion to this subsection, we note that both box-to-box and point-to-box dominance filtering appear to be promising tools for speeding up the resolution of NDMSP when p 2, and that box-to-box dominance filtering remains reliable when p 3, but beyond, it is beaten by the simple IF approach.

UPool lost only against box-to-box for p 4, but in this context, both methods performed slightly worse than IF. In all other cases, UPool performed similarly or better than box-based dominance approaches, and should thus be preferred.

Conclusions and discussion

When computing the non-dominated subset of the Minkowski Sum of a family of set, a great proportion of the sums of points generated at intermediary steps of computing the MS end up being dominated in the final result. This is why, as we have seen, filtering after having added each individual set reduces computing time hugely. We have sought to further decrease the resolution time of NDMSP by reducing the number of comparisons made at each intermediary step of the pooling process.

To this end, we proposed a unidirectional pooling algorithm, UPool, based on the lexicographic sorting of the input sets. This approach seeks to improve computation time relative to unidirectional dominance filtering algorithms by significantly saving on the initial sorting of the set sum. We also proposed using bound reasoning over boxes, i.e. subsets of the sets to be summed together. We defined two dominance relations: point-to-box and box-to-box. For p 2, filtering methods based on these two dominance relations have shown promising results, and the methods appeared to scale with the size of the sets. For p 3 however, only box-to-box dominance showed reliable benefits, and beyond, no box-based approach performed better than UPool or IF. Although we investigated the combining of both dominance relations, because optimal parameter values for creating boxes seem to differ between the two, we could not, so far, draw conclusive results, and further research on this topic is needed.

The main advantage of UPool as opposed to box-based methods is that it requires no parameter tuning on which the scaling of performance improvement with regards to size or number of criteria would depend. Moreover, in all presented experiments, UPool appeared to be the preferable option, as it performed better than box-based approaches in all cases but for p 4, the basic approach IF performed slightly better. In the next chapters, we will use UPool (adapted to the combination of solutions, not just points) as our default pooling method.

Chapter 3

Decoupling a coupled problem to obtain bound sets1

Chapter Abstract

In this chapter, we explore the notion of decomposition in multiobjective optimization. Because we deal with complex system problems which are not straightforwardly decomposable, we have to modify their original formulation to decouple them. This results in restricted problems, which provide lower bound sets to the non-dominated set, and relaxations, which provide upper bound sets to the non-dominated set. Decomposition allows the resolution of these decomposable variants in a fraction of the time needed to solve the original problem, an advantage that increases significantly with the number of dimensions of the objective space.

Introduction

In this chapter, we show how decomposability can be exploited in solving coupled MOCO problems. This is not straightforward since, typically, a coupled system is not such that ± sS X s X, where S is the set of subsystems and X s the feasible set of the subproblem associated with subsystem s S . While the original formulation may not be uncoupled, we can obtain decoupled variants of it by generic modifications. This chapter describes a number of these variants. Each of them is based either on a relaxation X I X or a restriction X P X of the feasible set X. The following enumerates these variants:

• Ignoring coupling constraints yields a relaxation of the original problem.

• Restricting the original problem to non-coupled variables yields a restriction of the original problem.

• Restrict-splitting coupling constraints, by which we mean fixing, for each coupling, the variables associated with each subsystem except one, yields a restriction of the original problem.

• Copy-splitting coupling constraints, by which we mean replacing each coupling constraint with |S | restrictions of itself, one to each subsystem s S , yields a relaxation of the original problem.

• More surprisingly, ignoring coupled variables in subsystems local constraints also yields an uncoupled problem. In this variant, the subsystems problems become uncoupled, and the coupled variables, bound only by coupling constraints, form a new, independent subproblem. The result is a relaxation of the original problem.

Bound sets

While the original formulation of an MOP may not be uncoupled, we can obtain decoupled variants of it through generic modifications. The set of non-dominated points associated with each of these variants will bound the set of images of efficient solutions of the original problem from either above or below.

Definitions

An application of relaxations is to obtain upper bound sets, and an application of restrictions is to obtain lower bound sets, which we will now define. The image y of an efficient solution of a relaxation delineates a zone of the objective space which is "unreachable", in the sense that no feasible solution can be found to dominate y.

©

In other words, in addition to ensuring that a strong upper bound set delineates an unreachable region, we also require that every point of the set which we want to bound from above be dominated by or equal to some element of the strong upper bound set.

Multiobjective relaxations and restrictions

Upper bound sets can be obtained by solving relaxations of the original problem, and lower bound sets by solving restrictions of the same problem. Informally, a relaxation is defined by replacing a collection of constraints of the original problem with a less constraining collection of constraints.

f 1 f 2 n 1 n 2 n 3 n 4 n 5 l 1 l 2 l 3 l 4 l 5 u 1 u 2 u 3 u 4 u 5 Figure 3.1: Let N tn i | i 1, ..., 5u, L tl i | i 1, ..., 5u, U tu i | i 1, ..., 5u. L is a
weak lower bound on N, U is a strong upper bound on N Definition 3.5. [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF] For f and f I taking values in R, max x I X I f I px I q is a relaxation of max xX f pxq if and only if X X I and for any x X, f pxq ¤ f I pxq.

We propose the following extension of the notion to the multiobjective case:

Definition 3.6. max x I X I pf I 1 px I q, ..., f I p px I qq is a multiobjective relaxation of max xX pf 1 pxq, ..., f p pxqq if and only if for all j t1, ..., pu, max x I X I f I j px I q is a relaxation of max xX f j pxq

A multiobjective relaxation yields a feasible set that is a superset of the original feasible set, so that no efficient solution of the relaxation can be dominated by a solution in the original problem.

Proposition 3.1. Let f pxq pf 1 pxq, ..., f p pxqq where f j takes values in R for each j t1, ..., pu. If max x I X I f I px I q is a relaxation of max xX f pxq, then Np f I pX I qq is an upper bound set on Np f pXqq.

Proof. Assume by contradiction that there exists y I f I px I q Np f I pX I qq and y f pxq Np f pXqq, such that y © y I . Then because X X I , we have x X I , and thus f I pxq ¡ = f pxq © f I px I q, so that y I Np f I pX I qq, contradicting our assumption.

A restriction of a problem with feasible set X is simply a variant with feasible set X I X.

Observation 3.1. For any X I X, Np f pX I qq is a lower bound set to Np f pXqq.

Proof. Assume that there is y I f px I q Np f pX I qq, and y f pxq Np f pXqq such that y I © y. Because X I X, we have x I X, contradicting the assumption that y Np f pXqq.

In the next two sections, we formally describe methods for obtaining variants of the original problem, each of which being either a relaxation, with a feasible set X I X, or a restriction, with a feasible set X P X, of the original feasible set X, while leading to a decomposable problem.

Decomposable restrictions and lower bound sets

Admissible variable values, neutral variable values and restrictions

In order to obtain an uncoupled restriction of the original problem, we need to fix all, or, as we shall see, most coupling variables to some admissible or neutral values. Given a feasible solution x X, we denote by x ¡i the vector decision variables where variable x i is omitted. Therefore, up to some permutation, we have x px ¡i , x i q. When focusing on a subset of I N of variables, we write x pxĪ, x I q.

Definition 3.7. Value e i is an admissible value for variable x i if there exists x X such that x i e i , or equivalently px ¡i , e i q X. Put simply, a value is admissible for a decision variable if we know that there exists a feasible solution in which the decision variable takes this value. This will help us fix the coupling to one possible interaction between the subsystems, considering one of the possible ways in which this interaction could occur, while maintaining the feasibility of solutions obtained under this assumption. A stronger notion is that of a neutral value assignment, by which we denote a value assignment that is feasible for any other variables assignment: Definition 3.8. Value e i is a neutral value for variable x i if, for all x X, px ¡i , e i q X. Observation 3.2. Let e i be a neutral value for variable x i , then e N : pe i | i Nq X.

For example, in the knapsack problem (KP), it is easily seen that some variable assignments do not alter the feasibility of a solution, and are thus neutral. Consider a partial feasible solution to KP, and an item i for which no decision has been made. Not taking the item, i.e. setting the corresponding decision variable x i to 0, can never yield an infeasible solution. But in general, we cannot assume that assigning the value 0 to a variable preserves feasibility.

Obviously, as long as the feasible set is non-empty, any neutral value is also an admissible value. Admissible values always exist, but must be determined as the corresponding variable value of a feasible solution, which implies that we find such a feasible solution. On the contrary, neutral values may not always exist, but can be substituted into the original problem formulation without having to first find a feasible solution. This is an advantage when a neutral value is known for a whole class of problem, such as 0 for any variable in an instance of KP.

Restrict-splitting coupling constraints

Consider again the formulation of the original problem, with feasible set X:

max f pxq s.t. x k X k d k K x s X s d s S (P) Let x s K px s
k | k Kq denote the vector of variables in subsystem s which are coupled by some coupling constraint. Thus x s px s K , x s K q, with x s K the vector of uncoupled variables in s. Then we will denote by x K the vector of all uncoupled variables.

The simplest way to obtain a decomposable restriction to the original problem is to fix all variables which appear in coupling constraints to an admissible value, or possibly a neutral value. Then, because no free variable remains linked by some constraint to free variables from other subsystems, the problem is obviously decomposable. However, the larger the proportion of coupled variables in the original problem, the poorer the quality of approximation of the non-dominated set of this original problem provided by the non-dominated set of this modified problem.

Intuitively, by fixing all coupled variables, we prevent any interaction between subsystems, we will call this weak lower bound the local restriction or LocRes lower bound. We can however, obtain decomposition by fixing only a subset of the set of coupled variables. Thus, in restrict-split variants of the original problem, we choose, for each k K, one subsystem s coupled by constraint k. We fix all coupled variables appearing in coupling constraint k, except those which belong to subsystem s. In this case, even if there were no strictly local variables, the restriction would not be reduced to a single feasible solution.

Let K s K be the subset of coupling constraints which involve variables from subsystem s, and in which we have chosen to keep free variables x s , while all other variables are set to some admissible value. For k K s , we denote by x s k the variables of subsystem s which appear in constraint k, and by x s k the variables of all subsystems except s, which appear in constraint k. Let x K f be the subvector of coupled variables which are left free. The restrict-split variant of problem pPq is:

max f px K , x K f , e K f q s.t. px s k , e s s k q X k ds S , d k K s x s X s d s S
which can be reformulated as k , e s s k q X k . In this restriction of the problem the constraints indexed by k are not coupling constraints anymore, because there is a single subsystem to which variables remaining free in the coupling constraint belong. Thus, because that S is a partition of N, we have that X I ± sS X Is . Finally, assum- ing that f pxq is separable along S , Corollary C.1 implies that this restriction of the problem is decomposable.

max f pxs K , x K f , e K f q s.t. px s s K , x s K s q X I s ds S
Example 3.2. The following illustrates the restrict split uncoupled variant of the original problem, obtained by fixing variables from subsystems other than s 1 which also appear in coupling constraint k.

Decomposable relaxations and upper bound sets

As specified by Definition 3.6, a multiobjective relaxation can be obtained either by replacing some objective functions f j by "relaxed" functions f I j or by considering a superset X I of the original feasible set X. We focus here on the latter possibility, which is to be achieved by relaxing the constraints defining X. The simplest way to achieve this is to simply ignore these constraints. The relaxations we wish to consider must also allow us to reach decomposability. For this purpose we need to ignore, in these constraints, only terms which are associated with some coupling variables, while ensuring that the resulting feasible set X I is a superset of X.

Consider N t1, ..., nu a set of variable indices, b R, gpx i | i Nq the function defining the left-handside of a coupling constraint, and feasible set X tx R n | gpx i | i Nq ¤ bu. Further consider I N. We denote px i | i Nq by x, px i | i Iq by x I and px i | i NzIq by xĪ. We say that g is relaxable for I if and only if there exists some function g I px I q such that g I px I q ¤ gpxq for all x in the domain of g. If g I relaxes g for I N, then for any

b I R such that b I ¥ b, X I tx R n | g I px I q ¤ b I u X.
Linear budget constraints (with positive coefficients) are common relaxable constraints: omitting any variable yields a relaxation of the original constraint.

Example 3.3. In the case of the knapsack problem, where N is a set of items, w i is the weight of item i N and b is the capacity, the budget constraint gpxq °iN w i x i ¤ b is relaxable for any subset I of objects. The relaxing constraint is that obtained by omitting all variables xĪ, so g I px I q °iI w i x i .

Copy-splitting coupling constraints

For each k K, let γpkq be the set of subsystems which have variables appearing in constraint k.

We consider a variant of the original problem obtained by producing as many relaxations of each coupling constraint k, as there are subsystems in γpkq. Each relaxation involves only variables from one of these subsystems.

Formally, assuming that coupling constraint k is relaxable for each s γpkq, we denote by g Is k px s k q ¤ b k the relaxation of g k px k q ¤ b k associated with s |γpkq|. Then the following problem is a relaxation of problem pPq, in the sense that its feasible set X I is such that X X I . max f pxq

s.t. g Is k px s k q ¤ b k d k K, s γpkq x s X s d s S
Setting X Is tx X s | g s k px s q ¤ b k dk K, s γpkqu the relaxation can be rewritten as: max f pxq s.t.

x s X Is Thus, we have X I ± sS X Is . Then, under the assumption that f is separable along S , Corollary C.1 implies that this relaxation is decomposable along S .

Example 3.4. The following illustrates the uncoupled variant of the original problem obtained by replacing coupling constraint k by two constraints k I and k P , which are copies of k which have been restricted to variables of each of the subsystems.

Local relaxation of coupled variables in subproblem constraints

In the next decouplable variant of the original problem pPq, we assume that each constraint x s X s is of the form g s px s q ¤ b s and is relaxable for s kK k. In other words, for K the set of noncoupled variables, g s px s q ¤ b s is relaxed into g Is px s K q ¤ b s . Applying this relaxation to all s S , we obtain the following relaxation:

s.t. x k X k d k K g Is px s K q ¤ b s d s S
Although the constraints for the subproblems now involve only strictly local variables, coupled variables still appear free in the objective function, and their assignments are still constrained by the coupling constraint, so the following formulation is a relaxation of the original problem.

C

defines a partition of N, with each set of constraints involving variables from one element of the partition, and thus

X I X K ¢ ¹ sS X s

K

Assume that for all j t1, ..., pu, f j is separable along tK, Ku, and that the term associated with K can be further separated along S , i.e.

f j pxq f K j px K q şS f s j, K px s K q
Then, Corollary C.1 implies that this variant is decomposable along S and along tK, Ku.

Example 3.5.

Experimental results

In this section, we report the results of an experiment serving two purposes. First, to report the computation times for our lower bound (LB) and upper bound (UB) concepts, and compare theses results with the computation time of the actual non-dominated set of the problem. Second, to evaluate how well they approximate this non-dominated set. We make these measurements as functions of the number of objectives (column p), the number of variables in the instance (column n), the number of subsystems in which the instance is subdivided (column (|S |), the number of coupling constraints in the instance (column |K|), and the proportion of variables in each subproblem that are coupled by each coupling constraint (column DV). We assume that when a constraint is coupling, it couples variables from all systems, at the rate defined by DV.

Experiments are conducted solving instances of the following generic coupled problem (GCP), which is formulated as: max f pxq p iN

π 1 i x i , ..., iN π p i x i q p
şS is

π 1 i x i , ...,
şS is

π p i x i q s.t. şS isk a i x i ¤ b k d k K is c i x i ¤ d s d s S x i t0, 1u d i N (GCP)
With f pxq : t0, 1u n Ñ N p . Because all constraints in the problem are budget constraints, for all i N, 0 is a neutral value for x i , which can then be used in all restrictions.

For every i N, we randomly choose π j i , a i , c i , r1, 1000s. When p 2, non decomposed original problems and subsystem problems are solved using a straightforward implementation of the bi-objective e-constraint method: while the problem is feasible we optimize lexicographically the first, then the second objective, under the constraint that each new solution is better than the previous one on the second objective. This is achieved by simply updating the lower bound of the same constraint. For p ¥ 3, we use the generic MOCO algorithm by [START_REF] Tamby | Enumeration of the nondominated set of multiobjective discrete optimization problems[END_REF]. Once solutions to subsystem problems have been obtained, they are pooled together according to Algorithm 6. At each step of the pooling process, some explicit dominance filtering needs to be performed. Figure C.9 exhibits the enframing of the original non-dominated set between the various notions of bounds presented in this article.

First, we study the increase in performance yielded by decomposition alone, in the case where the original optimization problem is uncoupled. In Table 3.1, we report the results of an experiment where we solved uncoupled instances without decomposition, and then using decomposition. Computing time for each of these operations is respectively denoted as T.NoDec, and T.Dec. Val- ues reported between squared brackets are standard deviations. We found that decomposition yields an improvement in computing time which increases with the number of criteria and with the number of variables when p 2, 3. For p 4, n 60 and |S | 3, resolution using decomposition took, on average, less than 1% of the time needed to solve the problem without decomposition.

We then consider uncoupled restrictions and relaxations of coupled instances. Solving a restrict split restriction of (GCP) requires us to produce some assignment of couplings to subsystems. For any k K, we need to select a subsystem s γpkq. The selection is made according to the following heuristic:

s k arg max sγpkq iks °p j1 π j i c i
In other words, we choose the subsystem coupled by k such that the variables of this subsystem which are coupled by coupling k maximize the average sum of profit-cost ratios over criteria.

Relaxed constraints for the copy split variant are obtained as follows. For any k K, if °sS °isk a i x i ¤ b k , then because coefficients a i are positive, for any s S , we have °isk a i x i ¤ b k , so that for any k K and s S , the latter constraint is a relaxation of the former. The copy split relaxation will thus be:

max f pxq s.t. isk a i x i ¤ b k d s S , dk K is c i x i ¤ d s d s S x i t0, 1u d i N (CS GCP)
Relaxed constraints for the locally relaxed variant are obtained by considering that, because coefficients c i are positive, we have °is s K c i x i ¤ °is c i x i , where s s K stands for sz kK k.

Thus °is s K c i x i ¤ d s is a relaxation of °is c i x i ¤ d s , which will be substituted to the latter in the locally relaxed variant. The locally relaxed variant will thus be:

max f pxq s.t. şS isk a i x i ¤ b k d k K is K c i x i ¤ d s d s S x i t0, 1u d i N (LR GCP)
To measure the quality of the approximation provided by a bound set, we use a posteriori ε-dominance, defined as follows : CHAPTER 3. DECOUPLING A COUPLED PROBLEM TO OBTAIN BOUND SETS Definition 3.9. Given y, y I Y, y is said to ε-dominate y I , written y © ε y I if and only if y j p1 εq ¥ y I j dj t1, ..., pu Definition 3.10. A is an ε-approximation of B for © ε if and only if dy B hy I A, y I © ε y

Here, we use this notion to evaluate the quality of approximation of a set A by a set B a posteriori. This is achieved by computing the smallest ε such that, for all y A, there is a point in the approximation B set which ε-dominates it. Formally, we compute

ε ¦ pA, Bq : mintε R | dy I A hy B, y I © ε yu
To measure how well a lower bound set L on a set Y approximates Y, we compute ε ¦ pY, Lq. To measure how well an upper bound set U to a set Y approximates Y, we compute ε ¦ pU, Yq. Another way to measure approximation quality is to measure the actual proportion of the non-dominated set of the original problem which is already contained in the bound set, which we will denote, e.g. in the case of an LB set, by |NDLB| |ND| .

In the following tables, RS denotes the restrict split restriction (Section 3.2.2) of the original problem, which provides a lower bound set. CS denotes the copy split relaxation (Section C. Let us begin by observing that in the bi-objective case, computing the bounds may not be relevant, because the whole non-dominated set is computed quickly. For more than two objectives however, as T.ND increases strongly, computing these bounds becomes worthwhile.

In Table 3.2, we report the average computing time T.ND and size of the non-dominated sets of the instances for which we computed bound sets, as well as their number of variables. This will serve as a reference to measure the advantages of computing bound sets using decomposition, rather than computing the original non-dominated sets.

From Table 3.3, we can observe the impact of varying instance parameter values on the time required to compute the RS lower bound. Increasing the number of subsytems while keeping n fixed makes the computation of the bound significantly faster, as the problem is divided into smaller subproblems, with an RS lower bound obtained in less than 2% of the time required to obtain the original non-dominated set when |S | 4. An increase in the number of coupling constraints |K| or in the proportion of coupling variables also speeds up the computation of the lower bound, because the greater the proportion of coupled variables in a problem, the smaller the number of variables in the restrict split variant. Approximation quality appears to be very good when |K| 1 and DV 10, with between 17 and 26% of non-dominated points contained in the RS lower bound when p 3 and between 27 and 51% when p 4. A slight deterioration in approximation quality is observed as |K| and DV increase, which is more readily explained by the neutralization of a higher number of variables. Now we turn to experiments with upper bounds. Results for the CS upper bound are reported in Table 3.4. The computation of this bound appears to require more time than does the RS lower bound. However, it benefits very similarly from an increase in n and in |S |, with a CS upper bound obtained in 1.2% of the time required to compute the original non-dominated set, for p 3 and |S | 4. An increase in |K| or in DV does not appear to have a significant impact, however, on computation time for the CS upper bound.

We observe that this upper bound consistently provides an excellent approximation quality of 1% when |K| 1 and DV 10, with on average between 34% and 47% of the non-dominated points already contained in the CS upper bound for p 3, and between 41% and 54% when p 4. Approximation quality decreases but remains good when |K| or DV increase.

Finally, turning to Table 3.5, we can observe that the LR upper bound is obtained strikingly faster than the set of non-dominated points: it required less than 1 % of the time in all our experiments. This most likely illustrates the power of dividing the original problem into a larger number of subsystems. Indeed in this case, instead of being associated with one of the |S | subsystems, the coupled variables constitute a separate, small subsystem. The impact of increasing n is similar to what was previously observed. An increase in |K| appears to have a detrimental effect on computing time for the LR upper bound, while an increase in DV has an unclear effect.

Surprisingly, an increase in |S | does not decrease the ratio between the computing time of the LR upper bound and that of the original non-dominated set. Indeed, it seems that increasing |S | also makes the global resolution of the original problem faster, albeit not to the same proportion as decomposition. However in the case of the LR upper bound, it appears that this effect was not as beneficial to computing the LR upper bound using decomposition.

On the downside, the LR upper bound provides a way poorer upper approximation of the nondominated set than the CS upper bound does. In all our experiments, the proportion of points from the original non-dominated set in the LR upper bound was 0, and the a posteriori value of ε was more than twice and up to ten times larger than for the CS upper bound.

In Appendix A, we propose the definition of a solution concept base on computing a lower bound and an upper bound on the non-dominated set of a MO problem. This double bound can either be used on it's own as an approximation of the non-dominated set of the original problem. The accuracy of this approximation can be defined either by ε-approximation, or by the difference in hypervolume dominated by the bounds. Alternatively, this double bound can be used as an element in a wider decision making process, providing information regarding the levels of performance that can be reached on each criterion. In the next chapter, we introduce an application problem that is a coupled problem, as well as the theoretical notions which are relevant to the study of its structure.

Conclusions and discussion

In this chapter, we have defined a collection of restrictions and relaxations of a coupled problem, which provide upper and lower bound sets to the non-dominated set of a coupled problem. The restrictions and relaxations were designed to be decomposable, and solved by solving independent subproblems and pooling the results.

On the one hand, we found that the restrict split and copy split bound provide very good approximations of the non-dominated set. The availability of decomposition to compute them provides a sizable advantage in terms of computation time, which increases as the number of variables and the number of subsystems increases. However, it may be argued that these bounds are still too expensive to compute.

On the other hand, we found that the locally relaxed upper bound can be obtained dramatically faster than the original non-dominated set, because it not only decomposes the problem along subsystems lines but also separates the coupled part of the problem from the rest. However, it provides a poorer approximation of the non-dominated set. A natural development of the work presented in this paper is to find ways to improve the quality of this upper bound concept while retaining the structural property that makes its computation very fast.

Part of the cost of computing bound sets comes from the sheer number of non-dominated points to be found, even for subproblems. When the goal is to produce an approximation of the non-dominated set of a problem, generating all solutions of the restriction or relaxation may not be useful. A priori ε-approximation (see e.g. [START_REF] Vassilvitskii | Efficiently computing succinct trade-off curves[END_REF], [START_REF] Bazgan | Approximate Pareto sets of minimal size for multi-objective optimization problems[END_REF]) could generate approximations of smaller size, but because it may contain dominated points, it cannot be used to approximate upper bound sets. Methods developed by [START_REF] Kaddani | Weighted sum model with partial preference information: application to multi-objective optimization[END_REF] can generate a subset of the non-dominated set of any problem, but to do so, it requires additional preference information from the decision maker. Knowing these challenges, further work should pursue the generation of subsets of decomposable upper bound sets.

Presentation of the REF problem

Original concepts developed in this work are intended to be applied to the resolution of a particular application problem, which we will refer to as REF. REF is an assignment problem where, to use classical terminology, tasks have to be assigned to machines in order to maximize some profit value associated with performing each task on some machine. Furthermore, machines are to be thought of as located on sites, with each site having limited capacity, which is consumed by assignments of tasks to machines located on it. The decision maker assumes the position of a central planner, who has to allocate tasks to machines on the different sites.

Formulation

REF can be dressed up as a multi-site production problem. Let T be a set of tasks, S a set of sites and M a set of machines. We will say that S ptq S is the set of sites on which task t can be performed. Symmetrically T psq T will be the set of tasks that can be performed on site s S . We assume that t T can in general be assigned to any machine of a site it can be performed on. Therefore, given Mpsq M the set of machines located on site s, we will denote by MpS ptqq the set of machines on which task t can be performed, i.e. for all t T, MpS ptqq sS ptq Mpsq.

Finally, N : tpt, mq T ¢ M | m MpS ptqqu will denote the set of all variable indices, with |N| : n.

Each assignment of a task t to a machine m yields profits g j tm , for j t1, ..., pu, and consumes an amount w tm of resources available on the unique site s S such that m Mpsq. This resource should be thought of as time or energy available to a manufacturing site. Each of the p objectives is to be maximized under the constraints that each task should be assigned at most once, and that for each site, resource consumption does not exceed some bound. Example C.5 provides a graph illustration of the assignment problem and examples of feasible and unfeasible solutions.

Example 4.1. Consider and instance of REF defined by T tt 1 , t 2 , t 3 u, M tm 1 , m 2 , m 3 , m 4 u and S ts 1 , s 2 u, such that S pt 1 q ts 1 , s 2 u, S pt 2 q ts 1 u, S pt 3 q ts 2 u, Mps 1 q tm 1 , m 2 u, Mps 2 q tm 3 , m 4 u. The following table records the weights associated with each assignment of a task to a machine. Assume that the capacity of each site is 1. A feature of this problem which captures our interest is that for some t T, |S ptq| 1, i.e. these tasks can be performed on one site only. We refer to the set of such tasks as T L , and to T zT L as T C . We will often call T L the set of local tasks, and T C the set of complex tasks. We will denote by T L psq the set T L T psq of local tasks which can be performed on site s. Thus we will say that subsystems or sites s and s I are coupled if there is a task t T C such that ps, s I q S ptq 2 . Sets of inequalities p2q and p3q represent local constraints, which involve only assignments of tasks which can be performed on one site only, and limit the consumption of resources on each site.

w tm t 1 t 2 t 3 s 1 m 1 0.6 0.5 m 2 0.5 0.5 s 2 m 3 0.6 0.5 m 4 0.5 0.6 t 1 t 2 t 3 m 1 m 2 p1q m 3 m 4 s 1 s 2 s 1 s 2 t 1 t 2 t 3 m 1 m 2 p2q m 3 m 4 t 1 t 2 t 3 m 1 m 2 p3q m 3 m 4 s 1 s 2 s 1 s 2 t 1 t 2 t 3 m 1 m 2 p4q

Relation to other problems and degenerate variants 4.1.2.1 REF and GAP

At first sight, REF bears a resemblance to the generalized assignment problem (GAP). Both problems combine assignment constraints and knapsack constraints. We temporarily focus on the single objective case, so as to use the complexity of single-objective GAP to characterize the complexity of REF. Let us recall the classical formulation of single objective GAP (cf. e.g. [START_REF] Kellerer | Knapsack Problems[END_REF]):

max n ķ1 m i1 g ki x ki m i1 x ki ¤ 1 d k t1, ..., nu n ķ1 w ki x ki ¤ b i d i t1, ..., mu x ki t0, 1u d i t1, ..., mu, d k t1, ..., nu
Closer comparison between GAP and REF reveals the following difference. Consider, as illustrated in Figure 4.2, a bipartite graph with set T , i.e. vertices corresponding to tasks on the lefthandside, and set M, i.e. vertices corresponding to machines on the right hand side. Each x ts t0, 1u d t T, s S which defines an instance of GAP. Thus, GAP reduces to REF, and is therefore NP-hard.

t 1 t 2 t 3 m 1 m 2 m 3 m 4 REF x 1,1 x 1,2 x 2,1 x 2,2 ¤ b 1 x 1,3 x 1,4 x 3,3 x 3,4 ¤ b 2 x 1,1 x 2,1 ¤ b 1 x 1,2 x 2,2 ¤ b 2 x 1,3 x 3,3 ¤ b 3 x 1,4 x 3,4 ¤ b 4 t 1 t 2 t 3 m 1 m 2 m 3 m 4 GAP

Multiple knapsack variant

Another particular case of REF is obtained when the assignments of some task yields the same profits and consume the same amount of resources for all machines on which the task can be performed. The problem then becomes that of selecting objects endowed with intrinsic value and weight, as in KP. Formally, we state that for any t T , m, m I MpS ptqq, we have g tm g tm I : g t and w tm w tm I : w t . We then get the following formulation: max ţT ş

A basic dynamic programming resolution method

Dynamic programming (DP), as introduced by [START_REF] Bellman | On the theory of dynamic programming[END_REF], is an algorithmic approach to solving decision problems which can be represented by multi-stage models of some decision process. The key feature of a DP algorithm is to build partial solutions incrementally at each decision stage. Stages can be organized in a sequence or a directed states graph [START_REF] Escoffier | Programmation dynamique[END_REF]). At each stage, all previous stages having been performed, a small optimization problem is solved, taking previous stages as input to generate new partial solutions. This optimization eliminates partial solutions which will not be built upon in later stages. On the one hand, this implies that some feasible solutions of the original problem cannot be reached anymore. On the other hand, this is a way to avoid considering all solutions to the original problem, and allows for a quicker resolution. The main challenge of DP is thus to define a sequential decision structure and intermediary optimization problems which allow the elimination of as many partial solutions as possible, while guaranteeing that efficient solutions remain accessible and are indeed outputted by the algorithm.

For simplicity, we will not define partial solutions on projections of the domain of complete solutions over subsets of variables. Rather we use the fact that 0 is an admissible value for any decision variable of REF. Recall that N : tpt, mq T ¢ M | m MpS ptqqu, and let 0 N be the solution such that for any i N, 0 i 0. By definition of an admissible value, 0 N is feasible. We set Q 0 : t0 N u to be the initial stage of the decision process. At each subsequent stage, building solutions according to decisions available at this step will amount to making an assignment of variables which previously took value 0. Thus, solutions generated during the sequential decision process are of equal dimensionality, and we can evaluate them without having to define one objective function per decision stage.

Let T t1, ..., tu, the set of tasks in REF, be used as the set of stages of the decision process, where t will denote the last decision stage. To each step t T , a DP algorithm must associate a set of decisions. In our case, decisions available for task t correspond to the machines m MpS ptqq to which t can be assigned. Additionally, the decision to assign task t to no machine must be available, and we denote it by e t . At step t T , taking decision m MpS ptqq boils down to assigning value 1 to variable x tm , and if decision e t is taken, no assignment is done, as x tm is set to 0 by default. Assigning value 1 to decision variable x tm in partial solution x will be denoted by xrx tm Ð 1s. Let Q t t0, 1u n be a subset of solutions built by taking decisions in periods 1 to t included, where only variables with indices pt I , m I q, for t I t1, ..., tu and m MpS pt I qq have been fixed according to decisions taken at these steps, while other variables keep value 0.

In general, for t T , consider x Q t and for ∆ t t : ± t t I t 1 MpS pt I qq te t I u, let δ ∆ t t denote a sequence of decisions from step t 1 to step t. We denote by τpx, δq the solution obtained by taking decisions prescribed by decision sequence δ. Finally Comppxq and Extpxq denote respectively the set of completions and the set of extensions of partial solution x:

Comppxq : 3 δ ∆ t t | τpx, δq X A Extpxq : tx X | hδ Comppxq, x τpx, δqu
To the final step t we associate the multiobjective function f pf 1 , ..., f p q, and to each step t t1, ..., t ¡ 1u we associate objective function pf, wq, where f is the multiobjective function in the original formulation of REF, and for any x X, wpxq b ¡ wpxq R |S| , is the function evaluating the residual capacity associated with solution x, where with b pb s |s S q is the bound vector of subsystem capacities.

Observation 4.1. For any t ¤ t, x, x I Q t , any δ Comppxq Comppx I q, f j pxq ¥ f j px I q ñ f j pτpx, δqq ¥ f j pτpx I , δqq dj t1, ..., pu and ws pxq ¥ ws px I q ñ ws pτpx, δqq ¥ ws pτpx I , δqq ds S Since several solutions can have the same value vector, a solution is not defined by its value vector. Iin a DP decision sequence however, the only information relevant to the sequential construction of a non-dominated set is the information defining which assignments it can still receive, and which value it can attain, and these will be inferred from additional criteria, as in [START_REF] Bazgan | Solving efficiently the 0-1 multi-objective knapsack problem[END_REF]. Thus two solutions with the same evaluation by the objective functions associated with the step at which they are generated are interchangeable from the point of view of DP. For this reason, Ẽ denotes a function that associates, to a each point of a non-dominated set, a single pre-image of this point.

Example 4.2 illustrates the general form and execution of a dynamic programming decision process.

Example 4.2. The following example illustrates a dynamic programming decision process, through the state-space graphs. At decision stages t 1 and t 2 , states Q 1 and Q 2 are generated by making feasible decisions from previously generated partial solutions, and then keeping only the nondominated newly completed partial solutions.

e N p0, 0q Q 0 Q 1 x 3 τpe N , δ 3 q p3, 3q
x 2 τpe N , δ 2 q p2, 2q

x 1 τpe N , δ 1 q p4, 1q δ3 p3, 3q δ2 p2, 2q δ1 p4, 1q © Q 2
x 6 τpx 3 , δ 6 q p6, 8q

x 5 τpx 1 , δ 5 q p5, 5q

x 4 τpx 1 , δ 4 q p7, 6q δ4 is not feasible in x3 δ5 p3, 5q δ5 p1, 4q δ4 p3, 5q © Figure 4
.3: Visualization of state-space graph of a dynamic programming decision process. T tt 1 , t 2 u. At stage t 1 , x 3 dominates x 2 , meaning that all the extensions of x 3 will dominate extensions of x 2 , so we avoid considering extensions of x 2 . We assume that decision δ 4 cannot be made from

x 3 , so that Q 2 contains no element τpx 3 , δ 4 q. At stage 2, x 5 is dominated by x 4 . Since t 2 is, in this case, the final stage, tx 4 , x 6 u is the set of efficient solutions.

Algorithm 16 describes the resolution of REF using dynamic programming. First, we initialize Q 0 with the zero partial solution te N u. Then for each task t T , we initialize Q t with value Q t¡1 , since in REF, the decision of not assigning task t is always feasible. Then for each partial solution

x Q t¡1 , we generate all assignments of tasks t which extend x and do not violate the resource constraint, i.e. have non-negative residual capacity w. Q t must then be filtered by dominance. If t is not the last step, we need to preserve state-elements which are not dominated for residual capacities. Only if t t do we eliminate elements dominated according to f .

: Q 0 Ð te N u output: ẼpX, f q 1 for t t1, ..., tu do 2 Q t Ð Q t¡1 3 for x Q t¡1 do 4 for m MpS ptqq do 5 x Ð xrx tm Ð 1s 6 /*Recall that for m I MpS ptqqztmu we already have x tm 0 */ 7 if wp xq 0 then 8 Q t Ð Q t txu 9 if t t then 10 Q t Ð ẼpQ t , pf, wqq 11 else 12 Q t Ð ẼpQ t , f q 13 return Q t
Note that because 0 is the default value of any variable, we do not have to explicitly generate partial solutions obtained by taking the neutral decision, as long as we begin each step t T by setting Q t Ð Q t¡1 . Next we show the correctness of this algorithm, i.e. that its produces the desired result: Proposition 4.2. Algorithm 16 returns the set Q t such that Q t ẼpX, f q Proof. : Let x ẼpX, f q, assume x Q t. Then there is some t I such that t I is the last step of the decision process at which some solution x I such that x Extpx I q also belongs to Q t I . If x I does not belong to Q t I 1 , then there is some x P Q t I such that pf, wqpx P q © pf, wqpx I q. Consider δ Comppx I q such that τpx I , δq x. Because wpx P q © wpx I q, we also have δ Comppx I q. Then because f px P q © f px I q, by Observation C.4, we also have, for x ¦ τpx P , δq, that f px ¦ q © f pxq. And thus, x ẼpX, f q, contradicting our assumption.

: Let x Q t, and assume x ẼpX, f q. This means that there exists x I X such that f px I q © f pxq. We can assume that x I ẼpX, f q, otherwise we replace it with the maximal element of a © chain it belongs to. Because ẼpX, f q Q t, x I Q t, and thus x Q t, being eliminated at line 12, contradicting our assumption. Algorithm 16 will be used as a subprocedure within many algorithms that are presented in this work. Thus, we use the following notation to describe an instance of REF and an application Table 4.1: Computing time (in seconds) using e-constraint and using Direct DP, for various instance parameter values and p 2, 3, 4. Average values for 10 trials. of Algorithm 16. We denote an instance by the tuple pQ 0 , T C , T L , M, S q, and we write B Ð DPpQ 0 , T C , T L , M, S q when we obtain a pre-image ẼpX, f q of Np f pXqq by solving the instance using Algorithm 16.

Experimental results

In all our experiments on problem REF, instances are generated as follows. We ensure that if |S ptq| ¡ 1, then S ptq S for any t T . In other words, if a task is complex, it can be performed on all sites. This is a way to both avoid considering cases where intermediary decomposition into several coupled subproblems would be possible, and to make sure that our performance measurements pertain to hard instances. For each pt, mq N and j t1, ..., pu, g j tm and w tm are both drawn uniformly randomly from interval r0, 1000s. Finally, we compute b s the bound of the capacity of site s S so as to adapt the notion of "hard" knapsack constraints, setting Unless otherwise stated, we will assume that |Mpsq| has the same value for all s S , and that |Mpsq| |T L psq|, i.e. that the subsystem assignment matrix is square. T.e ¡ const denotes the computing time to solve the instance using either biobjective e-constraint method for p 2 or the generic algorithm by [START_REF] Tamby | Enumeration of the nondominated set of multiobjective discrete optimization problems[END_REF] for p ¥ 3. Finally, T.DP Basic denotes computing time for the dynamic programming algorithm.

Experiments reveal the execution of this method to yield contrasted results. It performs well relatively to e ¡ const when p increases, but way more poorly as |T C | or |S | increase. The main reason for the increased computational cost, in the latter case, is the p |S | dimensionality of DP states, which results in the procedure generating a high number of mostly incomparable states until the very last stage of the process, where dimensionality is reduced to p. Several improvements can be made upon it to gain some speed, as will be demonstrated later by an experimental protocol comparing this basic approach to resolution methods that benefit from some improvements using decomposition. But as these tricks use upper and lower bound sets, we need first to define those for REF, and describe algorithms to obtain them.

Upper and lower bound sets for REF

Weak upper and lower bound notions developped in Chapter 3 can be straightforwardly applied to REF. Furthermore, these bound sets can be computed using dynamic programming algorithms.

Weak upper bound sets and a strong singleton upper bound 4.2.1.1 Copy-split relaxation

We recall that in the copy-split relaxation, each subsytem problem inherits a copy of each coupling constraints, which is restricted to variables associated with complex tasks assignments than can be made in it. In the notation introduced in Section C.3.3.1, we assume that coupling constraints imposed by coupling k K can be reformulated as g k px k q ¤ b k for some real valued function g. We further assume that for j t1, ..., pu, f j is additively separable along S and that so is g k for any k K. Recall that for any k K, γpkq ts S | s k ru. Recall that assuming that coupling constraint k is relaxable for each s |γpkq|, we denote by g Is k px s k q ¤ b k the relaxation of g k px k q ¤ b k associated with s |γpkq|. Then the copy split relaxation can be expressed as:

max f pxq s.t. g Is k px s k q ¤ b k d k K, s γpkq x s X s d s S
In REF, K tpt, mq | t T C , m MpS ptqqu and g s k px s k q °mMpSptqq x tm for t T C . For any t T , b t 1 . Then the copy-split relaxation of REF can be formulated as follows: max ţT mMpSptqq x tm ¤ 1 d t T, s S ptq ţTpsq mMpSptqq x tm ¤ 1 d t T psq ţTpsq mMpSptqq w tm x tm ¤ b s

w tm x tm ¤ b s d s S x tm t0, 1u d t T, m MpS ptqq t 1 t 2 t 3 m 1 m 2 m 3 m 4 t 1 t I 1 t 2 t 3 m 1 m 2 m 3 m 4 s 1 s 2 s 3 s 4
x tm t0, 1u d t T psq, m MpS ptqq

As this subproblem is itself an instance of REF, it can be solved using Algorithm 16. Then, we compute copy-split relaxation using Algorithm 17.

Algorithm 17: Dynamic Programming Algorithm for solving the copy-split relaxation of REF with decomposition. input : pte N u, T C , T L , M, S q 1 S ubRes Ð pr 1 , ..., r |S| q 2 for s S do As shown in Section C.3.3.2, a less obvious way to relax the problem is to relax, from local constraints, decision variables which appear in coupling constraints, a relaxation which we have called the LR relaxation. We recall that for K the set of non-coupled variables, g s px s q ¤ b s is relaxed into g Is px s K q ¤ b s . Applying this relaxation to all s S , this variant is formulated as follows, max f pxq

s.t. x k X k d k K g Is px s K q ¤ b s d s S
mMpSptqq x tm ¤ 1 d t T C mMpSptqq x tm ¤ 1 d s S , t T L psq ţT L psq mMpSptqq w tm x tm ¤ b s d s S x tm t0, 1u d t T, m MpS ptqq
This problem is separated in |S | 1 subproblems, the first |S | of which involve only tasks local to system s S . Each of these tasks must be assigned to at most one machine of site s under the budged constraints of s : max ţT L psq mMpsq

g j tm x tm d j t1, ..., pu s.t. mMpsq x tm ¤ 1 d t T L psq ţT L psq mMpSptqq w tm x tm ¤ b s x tm t0, 1u d t T L psq, m MpS ptqq
Each of these problems can be solved using DP, computing DPpte N u, r, T L psq, Mpsq, tsuq, for s S . Then, there remains the now-independent problem of assigning each complex task to at most one machine, under no budget constraint:

CHAPTER 4. APPLICATION PROBLEM max ţ T C mMpSptqq g j tm x tm d j t1, ..., pu s.t. mMpSptqq x tm ¤ 1 d t T C x tm t0, 1u d t T C , m MpS ptqq
This problem, the complex task relaxation problem, can also be solved using a DP algorithm, albeit a much simpler one that the one used for solving subproblems, because there is no need to check for the feasibility of partial solutions with regards to budget constraints. This simplified procedure is presented as Algorithm 18, the application of which is denoted by DPpte N u, T C , r, M, S q.

Algorithm 18: Dynamic Programming Algorithm for the complex tasks assignemnts in the LR relaxation of REF

input : pte T C u, T C , r, M, S q 1 Q 0 te T C u 2 /*Assume |T C | t */ 3 for t t1, ..., tu do 4 for x Q t¡1 do 5 for m MpS ptqq do 6 x Ð xrx tm Ð 1s 7 Q t Ð Q t txu 8 Q t Ð ẼpQ t , f q 9 return Q t
Thus, Algorithm 19 solves the LR relaxation by pooling results of the locally restricted subproblems and the complex tasks problem, using any pooling algorithm, here denoted generically by Pool.

Algorithm 19: Dynamic Programming Algorithm for solving the LR relaxation of REF

input : pte N u, T C , T L , M, S q 1 Complex Ð DPpte N u, T C , r, M, S q /*

Singleton strong upper bound

Improvements on the basic resolution method for REF can be made by eliminating partial solutions when a strong upper bound on the value of their extensions is found to be dominated by an incumbent point. Because these must be computed quite often, we look for a greedy algorithm to obtain them.

In the case of the single objective knapsack problem, it is well known (see e.g. [START_REF] Kellerer | Knapsack Problems[END_REF]) that the linear relaxation upper bound can be obtained using the following greedy algorithm. Objects t t1, ..., tu are ranked in decreasing order of ratio g t w t . Objects are added in order until residual capacity becomes negative, and the resulting sum of profits is an upper bound on the value of the optimal knapsack. This can be improved further using [START_REF] Martello | Lower bounds and reduction procedures for the bin packing problem[END_REF]'s concept of upper bound. This concept of singleton upper bound can be applied to any partial solution where objects t1, ..., tu have been previously decided upon, for t t. Objects upon which no decision has yet been taken are then added, in an order determined by the previously defined ratio. Finally, point pu G 1 , ..., u G 1 q obtained by computing the objective-wise greedy upper bounds is a strong upper bound on all non-dominated knapsacks, or extensions of a partial knapsack.

However none of these bound concepts straightforwardly applies to our case, even when assuming that |S | 1, as is evident from the following counter example.

Example 4.3. Let T t1, 2u, S tsu, Mpsq ta, bu, assume that b s 6 and consider the following matrix, where an entry reads as g tm {w tm , for t T and m Mpsq. The following sequence of pairs describes the ranking of assignments in decreasing order of ratio g tm w tm : p1, aq, p1, bq, p2, aq, p2, bq. The greedy algorithm takes p1, aq, after which the knapsack has value of 20 and a residual capacity of 2. p1, bq cannot be added to the knapsack because task 1 has already been assigned to machine b. Then, all remaining assignments would violate the knapsack contraints, so the knapsack is full. The candidate upper bound value of 22 is obtained by further adding the value of p2, aq, the split assignment. However, the optimal solution to this instance of the problem is to select p1, bq, for a profit of 24, which is higher than 22. Thus, the greedy algorithm failed to yield an upper bound solution.

We propose a greedy heuristic which gets around the issue of the assignment constraint. First, we copy complex tasks, so that we can decouple the problem as in the copy-split relaxation. This yields T I such that, for each t T C , for each s S ptq, and m Mpsq, we have t I T I L psq with w t I m w tm and g j t I m g j tm for each j t1, ..., pu. Then we replace T I with a set of artificial tasks For each j t1, ..., pu, we compute the j-th component of the objective-wise optimum as follows. For each s S , we rank T ¦ psq in decreasing order of g j¦ t w ¦ t , and we compute a candidate objective-wise upper bound by making assignments in order, until the assignment which exceeds capacity of site s has been made. We denote by u ¦ j the objective value of this greedy solution on component j. Proposition 4.3. pu ¦ 1 , ..., u ¦ p q computed by the previously described procedure is a singleton strong upper bound to the set of extensions of a partial solution.

Proof. Consider, for any instance of REF, its copy-split relaxation, with feasible set X I and objec- tive function f I . Then consider the variant of the copy-split relaxation where for any pt, mq N, g j tm is replaced by g j¦ t and w tm by w ¦ t . Since for any pt, mq N, we have w ¦ t ¤ w tm and g j¦ t ¥ g j tm , this variant with feasible set X P and objective function f P is a relaxation of the instance with feasi- ble set X I . max xX I pf I 1 pxq, ..., f I j pxqq is itself a relaxation of the instance of REF, and by transitivity max xX P pf P 1 pxq, ..., f P j pxqq is a relaxation of X. Now, let us show that pu ¦ 1 , ..., u ¦ p q is a singleton strong upper bound over the set Np f P pX P qq. Indeed, for s S , the subproblem associated with s of the problem with feasible set X P is equivalent to an instance of the knapsack problem, since the assignments of a task t T now all yield the same profit and have the weight value for any accessible machines in Mptq. For each j t1, ..., pu, component u ¦ j of greedy assignments according to decreasing order of g j¦ t w ¦ t (where the limit item is taken in the knapsack) is an upper bound to the linear relaxation of the same knapsack problem. Therefore, u ¦ j ¥ max xX I f P j pxq, and by transitivity applied twice, u ¦ j ¥ max xX f j pxq for all j t1, ..., pu.

Strong and weak lower bound sets 4.2.2.1 A strong lower bound set

We consider the feasible set X L of the restriction of the original problem to local task variables, i.e. X L tx X | x tm 0 dt T C , m MpS ptqqu As a particular case of the notions we 92

CHAPTER 4. APPLICATION PROBLEM t 1 t 2 t 3 m 1 m 2 m 3 m 4 s 1 s 2 s 1 s 2 t 1 t 2 t 3 m 1 m 2 m 3 m 4 Figure 4
.6: Visualization of the local restriction as applied to REF: assignments of complex task t 1 are removed from all subsystems studied in Chapter 3, it is clear that Np f pX L qq is a weak lower bound on Np f pXqq. Furthermore, the computation of ẼpX L , pf, wqq is separable, since X L Π sS X Lpsq for X Lpsq tx s X s | x s tm 0 dt T C , m Mpsqu. When computing the efficient subset of X L for f , we eliminate all solutions dominated for the p original criteria, even if these solutions are non-dominated when residual capacity is taken into account. If, to the contrary, we keep the latter solutions, but still evaluate them according to the p original criteria, we obtain set f p ẼpX L , pf, wqqq R p , which lies in the same objective space as Np f pXqq and has the following property: Proposition 4.4. f p ẼpX L , pf, wqqq is a strong lower bound set to Np f pXqq.

Proof. Assuming T L t1, ..., lu, ẼpX L , pf, wqq corresponds to Q l in Algorithm 16. This means that for any x ẼpX, f q, there exists x I Q l such that x Extpx I q. Therefore, we have f pxq © f px I q. Thus for any f pxq Np f pXqq, there exists f px I q f p ẼpX L , pf, wqq such that f pxq © f px I q, establishing the result.

Further on, we will call this strong lower bound Strong LocRes, as opposed to the weak variant of the same restriction, which has the same set of feasible solutions, but is solved only for the main p criteria. Let us recall that a strong lower bound set is not necessarily better than a weak lower bound set when used as an incumbent set. Rather, it is a guarantee that all non-dominated points lie in the union of the upper orthants associated with its elements, rather than in the incomparable regions around it. In practice, because EpX L q contains close-to-empty-solutions which maximize residual capacity, its image Np f pX L qq is not a very useful strong lower bound set. However, using singleton upper bounds, computed by taking into account all remaining complex tasks, i.e. all t T C , we can improve on it significantly. Namely, we can remove from ẼpX L , pf, wqq any solution, such that a strong upper bound on its set of extensions is dominated by another solution from the same set. Formally: Proposition 4.5. Let x EpX L , pf, wqq, u x a singleton strong upper bound on Extpxq, and © u : tpx I , xq X 2 | f px I q © u x u. Then f pEp ẼpX L , pf, wqq, © u qq is a strong lower bound set on Np f pXqq. In other words, filtering EpX L , pf, wqq by © u yields a strong lower bound set to Np f pXqq.

Proof. By Proposition 4.4, for each x ẼpX, f q, there exists x I ẼpX L , pf, wqq such that

x t 1 t 2 t 3 m 1 m 2 m 3 m 4 t 1 t 2 t 3 m 1 m 2 m 3 m 4 s 1 s 2 s 1 s 2 Figure 4
.7: Visualization of the restrict-split restriction as applied to REF: assignments of complex task t 1 are restricted to site s 1 . In this case, we have T ¦ s 1 tt 1 u and T ¦ s 2 r

.

Extpx I q. For x EpX, f q, let Restpx, Lq : tx I X L | x Extpx I qu, and RestpEpX, f q, Lq xEpX, f q Restpx, Lq. As long as filtering by © u removes no element from RestpEpX, f q, Lq, it yields a strong lower bound set. For any x I EpX L , pf, wqq, u x I is a strong upper bound on Extpx I q, meaning that for any x Extpx I q, u x I © x. Thus for any x I , x P EpX L , pf, wqq if f px P q © u x I , then there cannot be any x EpX, f q such that x Extpx I q, and thus x I RestpEpX, f q, Lq.

Weak lower bound obtained from solving the Restrict-split variant

The adaptation of the restrict-split variant is fairly straightforward. For each t T C , we choose one s S , and thus for each s S we obtain a subset of T ¦ s T C of complex tasks which become local tasks for s and are thus added to the set of tasks T L psq, as seen line 4 in Algorithm 20 and illustrated by Figure 4.7. Profits and weights of assignments which are kept from the original instance by this operation remain unchanged.

Algorithm 20 describes the computation of the restrict-split variant of REF using dynamic programming. Once we have computed the heuristic assignment of complex tasks to subsystems, we modify each subsystem instance, solve them independently, and pool the resulting solution sets together.

Algorithm 20: Dynamic Programming Algorithm for solving the copy-split relaxation of REF with decomposition. input : pte N u, T C , T L , M, S q output: Np f pX I qq 1 S ubRes Ð pr 1 , ..., r |S| q 2 Let pT ¦ s | s S q be the list of subsets of complex tasks restriced to each s S 3 for s S do perform poorly, but it can easily be improved. Indeed, it may contain solutions which still have enough residual capacity to receive assignments of complex tasks. Even more such assignments could be made on solutions of some subset of strong LocRes, i.e. some subset of ẼpX L , pf, wqq, but doing all possible non-dominated assignments on the whole strong LocRes would amount to solving the original problem, and our goal here is only to obtain a weak LB. We need to strike the right balance between the richness of the set from which we complete solutions, the speed of computing extensions, and the tightness of the obtained weak LB.

Let κpBq denote the set of solutions obtained by completing partial solutions from some set B with further feasible assignments. In order to obtain an incumbent set quickly, κ should be computable by a greedy algorithm. We use a ranking of complex tasks assignments variables to complete, when possible, solutions obtained after combining solutions of the local problems.

Following [START_REF] Bazgan | Solving efficiently the 0-1 multi-objective knapsack problem[END_REF], we sort assignments by decreasing order of g j tm w tm , for each objective j t1, ..., pu. We then globally sort the assignments in increasing order of the sum of their ranks in objective-wise orders, as described by Algorithm 21.

For each solution x B, we make feasible assignments of tasks to machine according to the previously defined ranking. Because a task can be assigned only once and within site budget constraint, when a task has been assigned, we can clear the ranked list of its other assignments. We remove all assignments to other machines of the task having just been assigned, and all assignments to machines on a site that would have been saturated at that step.

It appears, from preliminary experiments, that κp ẼpX L , pf, wqqq yields a far richer incumbent set than κp ẼpX L , f qq, because the solutions it contains admit more completions from additional variable assignments, since they belong to this set because of their low budgets consumption. However, computing κp ẼpX L , pf, wqqq without a preliminary filtering proves very expensive, while many quasi-empty solutions cannot get good enough completions. We can modulate the richness of the set on which to apply greedy completions in three ways.

First, remembering that X L ± sS X L , we may compute, for each s S , ẼpX s L , f q rather than ẼpX s L , pf, wqq. Second, for B s the set of solutions computed for decomposable steps associated with s S , rather than computing Npp f, wqp ± sS B s qq, we may compute Np f p ± sS B s qq.

Note that even if we have computed ẼpX s L , f q as suproblem solutions, we can still compute Ẽpp ± sS B s q, pf, wqq, which will be richer than Ẽpp ± sS B s q, f q.

Third, we can require that the f -dominated solutions to be completed by the greedy algorithm still score better than some minimal f values. Indeed, many solutions from a sufficiently rich set, be it either Epp ± sS ẼpX s L , f qq, pf, wqq or Epp ± sS ẼpX s L , pf, wqq, pf, wqq must still have too low values on f for their completions to yield non-dominated solutions. We can filter them out using an approach akin to ε ¡ dominance. We will say that x α¡dominates x I if for all j t1, ..., pu, f j pxq ¥ f j px I q α ¦ M, where M could be the maximum possible profit value, or could be f j px I q. We denote this filtering operation by ẼpB ¦ , αp f qq. Finally one has to choose an α. A relevant value is one such that the obtained incumbent, κpEpB ¦ , αp f qqq, is nearly as rich as κp ẼpB ¦ , pf, wqq, but is obtained in a sufficiently short time.

Experimental Results

The following experiments compare the time necessary to compute two weak lower bounds: the solution of the restrict-split variant of REF, and the greedy lower bounds described in the previous section. We will also assess the quality of these bounds as approximations of the original REF instance. |Mpsq| will denote the number of machines on each site s, assuming that all s, s I S , |Mpsq| |Mps I q|. Once again, the matrix of assignments of local tasks to machines of each site is square, with, for all s S , |T L psq| |Mpsq|. |S | denotes the number of subsystems, hence of subproblems. |N| is the number de decision variables in the instance.

In Table 4.2, ApQ.(ε) denotes, as in Section C.3.4, the approximation quality of a set as measured by a posteriori ε-dominance. We did not perform experiments using the greedy completions of ẼpX L , pf, wqq, because computing time for this weak lower bound proved to be prohibitive for the instance sizes considered. Rather, we compared weak lower bounds obtained by greedily completings three types of base solutions sets.

First, GrS LocRP denotes Ep ± sS ẼpX s L , f q, pf, wqq, which is obtained by solving the subproblems for criteria f , but keeping the combinations of subsystem solutions which are non-dominated for pf, wq. Second, GrS LocRPα denotes EpEp ± sS ẼpX s L , f q, pf, wqq, αp f qq, which is the same as the previous one, where the result of the pooling is further filtered according to α-dominance as described in Section 4.2.2.3. We use α |T C | |T L | , and relation x © α x I ô f pxq © f px I q α ¢ gMax where, for j t1, ..., pu, gMax j 1000 is the upper bound of the interval in which g j tm is randomly drawn, for any pt, mq N. Third, GrWLocR denotes the greedy completions of ẼpX L , f q, which is simply the weak LocRes lower bound. From Table 4.2, we can observe, first, the expected: the richer the set from which the completions are computed, the better the approximation quality (ApQ.(ε), to be minimized), but the more expensive the computation becomes. Because the tradeoff between these two measurement appears to be rather balanced, it is difficult to decide in advance which of these variants of greedy weak lower bound is the best. Comparing results for the greedy lower bounds with results for the restrict-split variant in Table 4.3, it appears that the latter dominates the former, in the sense that it provided a better approximation quality, and was computed quicker for all considered instance parameters. It is particularily interesting that it was computed even faster than the GrWLocR greedy lower bound, which was obtained from greedily completing the solution of a uncoupled restriction with fewer variables, and therefore solved faster than the restrict-split variant. This gives us good reasons to believe that in applications, the latter will perform better than the former. However, this remains to be tested empirically.

In Appendix B, we present additionnal lower bound concepts, which are based on correcting upper bound solutions obtained by solving either the CS or the LR relaxations. These corrections are performed according to greedy heuristics, removing in priority assignments which minimize the aggregated profit to weight ratio, until constraints of REF that were violated by the solutions to the relaxation become satisfied again.

In the next sections, we investigate ways in which the DP approach to solving REF can be improved, using, among other devices, decomposition in independent decision sequences, incumbent sets and singleton upper bounds. In Section C.4.3, we propose improvements over the basic DP method which we characterize as "Bottom-up" approaches, based on using decomposition in the treatment of local tasks first, and delaying as much as possible the consideration of all residual capacity criteria at once. In Section C.4.4, we tackle choices associated with complex tasks first, so as to decouple the problem entirely before solving instances defined by these decoupling choices. Hence we characterize those approaches as "top-down".

"Bottom-up" approaches based on Dynamic Programming

In this section, we present improvements of the DP algorithm which take advantage of the fact that some subsequences of steps in the decision process are independent. In our case, subsequences of steps associated with local tasks from each subsystem. This can take several forms. First, the precomputation of an incumbent set from a decomposable restriction, which can be used in bound reasoning to eliminate partial solutions. Second, a preliminary filtering of each state of the decision sequence can be done considering only p 1 criteria, by partitioning the set of partial solutions according to the site in which the last assignment has been made. Third, we can initialize the DP process with a strong lower bound set obtained from the decomposable restriction of the problem to local tasks, effectively ignoring |S | ¡ 1 criteria in all decision steps associated with local tasks.

Filtering with a set of incumbent solutions

At each step of the decision sequence, partial solutions can be eliminated by comparing strong upper bounds to the values of their extensions to the values of already known solutions. This is similar to fathoming a node in branch and bound, where a subtree can be avoided if it is shown that an upper bound over all solutions in this subtree is dominated by an incumbent solution. If an incumbent solution is found to dominate the strong upper bound on extensions of a partial solution, then all the extensions of this partial solution will be dominated by the known incumbent solution. It is thus useless to pursue extensions of this partial solution, and we filter it out. This procedure, akin to what was proposed by [START_REF] Figueira | Algorithmic improvements on dynamic programming for the bi-objective t0, 1u knapsack problem[END_REF], is described in its application to REF by Algorithm 22. The better the precomputed incumbent solutions, the more powerful this filtering mechanism, in terms of the number of extensions eliminated by a successful test. However, the cost of computing the incumbent has to be taken into consideration, as it weighs on the overall performance.

Partial filtering

When generating new feasible solutions during the DP process, it is possible to filter subsets of these new solutions without considering all |S | residual capacity criteria. In particular this is the case when two solutions which did not previously dominate each other are modified by assignments consuming resources in the same subsystem. Consider, for t T , s S ptq, Q t psq the set of states where task t was assigned to a machine of site s. Proposition 4.6. for x, x I ẼpQ t¡1 , pf, wqq such that τpx, δq, τpx I , δ I q Q t psq, pf, wqpτpx, δqq ¨pf, wqpτpx I , δ I qq if and only if pf, ws qpτpx, δqq ¨pf, ws qpτpx I , δ I qq

Proof. Left to right is trivial. For right to left, consider that if τpx, δq, τpx I , δ I q Q t psq, then for all s I s, ws I pτpx, δqq ws I pxq and ws I pτpx I , δ I qq ws I px I q. In other words, on the one hand, if ws pτpx, δqq ¨w s pτpx I , δ I qq, then wpτpx, δqq ¨wpτpx I , δ I qq, and on the other hand, if pfpτpx, δqq, ws pτpx, δqq ¨pfpτpx I , δ I qq, ws pτpx I , δ I qqq, then in particular f pτpx, δqq ¨fpτpx I , δ I qq.

The conjunction of these two facts yields the result.

This observation allows the following modification of the DP procedure, described as Algorithm 23, to produce the correct result. In essence this modification applies a sort of divide and conquer heuristic to filtering set of states Q t , the latter being partitioned into |S | sets of states, to which preliminary filterings are to be applied.

= 0 then 7 Q t psq Ð Q t psq txu 8 Q t psq Ð ẼpQ t psq, pf, ws qq 9 Q t Ð sS Q t 10 if t t then 11 Q t Ð ẼpQ t , pf, wqq 12 else 13 Q t Ð ẼpQ t , f q 14 return Q t 4.3.3 "Kickstarting"
Assume that T , which is both the set of tasks and the set of steps in the DP decision process, is ordered in such a way that all t T L come before all t T C , and that t C is the first complex task. Q t C can be computed by solving the subproblems associated with each s S for criteria pf, ws q and pooling the results together. Let l denote the last element of the set of local tasks, then according to Algorithm 16,

Q t c Ẽptτpx, mq X t c | x Q l , m MpS pt c qq e t C u, pf, wqq
What we propose is a way to compute Q l using decomposition. The subsequence of the DP decision process associated with T L can be seen as a DP process where both the intermediary and target objective function are pf, wq. Thus Proposition C.3 implies that Npp f, wqpX L qq pf, wqpQ l q X L is exactly the feasible set of the local restriction of REF, and we know from section 3.2.2 that X L ± sS X s L . Again by Proposition C.3, for any s S , T L psq t1 s , ..., l s u the set of local tasks associated with site s, we have that Npp f, wqpX s L qq pf, wqpQ l s q Algorithm 24: "Kickstarted" dynamic programming for solving REF

input : pT C , T L , M, S q output: ẼpX, f q 1 for s S do 2 Q l s Ð DPpr, T L psq, Mpsq, tsuq 3 Q t C Ð Ep ± sS Q l s , pf
Q t Ð Q t txu 11 if t |T C | ¡ 1 then Q t Ð ẼpQ t , pf, wqq 12 else Q t Ð ẼpQ t , f q 13 return Q |T C | 4.3.

Experimental Results

In experiments reported in Tables 4.4 to 4.6, we measure the performance improvement yielded by implementing the modifications of the DP algorithm which we described in the previous subsections. The implementation of partial filtering (denoted by Pt.Filt) in these experiments does not require additional clarifications. However, to apply the filtering of partial solutions using an incumbent set (denoted LB.Filt), as described in Section C.4.3.1, we first need to choose among of the available notions of strong upper bounds and weak lower bounds to constitute our incumbent. Our chosen notion of strong upper bound over the extensions of a partial solution will be the greedy singleton upper bound which we described in Section C.4.2.2 and proven to be an strong upper bound in Proposition C.4. We used the two types of lower bound which we have already experimented with in Section 4.2.3. Namely, the GrS LocRα lower bound, and the RS lower bound, which is the solution of the restrict-split restriction. Further work is required for filtering using an incumbent set in the kickstarted approach, using an upper bound that can eliminate partial solutions in earlier steps of the decision process. We have observed that using the objective-wise optimum of the LP relaxation of the completion problem tends to eliminate partial solutions earlier (and fewer in later stages), but it appeared too costly in preliminary experiments. Yet, it hints at a possible way of improving this method.

To summarize, all dynamic programming-based approaches suffer major impairments when |S | increases. This reflects the fact that an increasing |S | generates more incomparabilities between partial solutions, which increases the number of partial solutions to be handled until the very last stage. This remains the main flaw of the "bottom-up" approach to solving REF, which we have not yet found a way to overcome. As a consequence, we will now explore methods which use DP when only one residual capacity criterion is involved, i.e. when the problem is reduced to a collection of restrict-split variants, as will be investigated in the last section of this work.

In the next section, we shift the perspective on solving the REF problem. We started by solving it from the bottom up, starting from local task assignments and keeping sufficient residual capacity to assign complex tasks. Next, we will investigate how preliminary coordination of the assignment of complex task may allow for a faster resolution of REF than that available with generic methods.

"Top-down" approaches for solving REF

A central hypothesis in the decomposition approach to the optimization of complex systems is that a collection of independent subsystem problems is significantly easier to solve than the whole problem. If this discrepancy is sufficiently large, in particular if subproblems can be solved using methods which specifically exploit their structural properties to reach higher efficiency, it may be possible to solve not just the collection of subproblems associated with one variant, but several of such collections, and in less time than a generic algorithm requires to solve the whole original problem. Under some conditions, solving such collections of variants can lead to an exact solution to the original MOIP problem.

Enumerating all ways of decoupling REF

Recall from Section 3.2.2 that the restrict-split variant of any instance of the GCP is defined by assigning each coupling constraint to some subsystem s S , yielding subset K s K of coupling constraints involving variables from subsystem s, and in which we keep only variables x s free. The restriction is thus formulated as

max f px K , x K f , e K f q s.t. px s k , e s s k q X k ds S , d k K s x s X s d s S or equivalently, max f pxs K , x K f , e K f q s.t. px s s K , x s K s q X I s ds S where X I s X s kK s X I k , and X I k is a set such that x s k X I k if and only if px s k , e s s k q X k .
Note that the set S K of applications of K into S defines the set of all possible assignments of couplings constraints to subsystems. Thus to each σ S K corresponds one restrict-split variant, and for X σ the feasible set of any such variant, we have that X σ X. If we can prove, for a particular coupled problem, that σS K X σ X, then the following fact ensures that we can solve the original problem by solving all possible decoupled variants of the problem. To summarize, decoupling a problem amounts to fixing some variables to admissible values, effectively turning coupling constraints into non-coupling, local constraints. Using notation previously introduced for restrict-split variants of coupled problems, we saw that the set S K of functions σ associating each coupling to a particular subsystem describes an exhaustive exploration of all ways of decoupling the original problem. Tables 4.9 and 4.10 present results of an experiment where an instance of REF is solved by enumerating all of S T C (since K T C in this case). Instances are such that each site has the same number |Mpsq| of machines. For each σ S T C is solved as a whole (NoDec), and then enumerating all restrict-split variants (Enum), which are solved by decomposition, subproblems being solved using e-constraint method for p 2, and the generic algorithm by Tamby and [START_REF] Tamby | Enumeration of the nondominated set of multiobjective discrete optimization problems[END_REF] for p ¥ 3. This experiment shows that even in the biobjective case, enumeration of S T C (column Enum&Dec) can be competitive and even faster than generic resolution when the number of subsystems is relatively low compared to the size of subsystems. However, for smaller subsystems but a larger value of S T C , this approach performs worse than the generic method.

Pre-computation of independent dynamic programming states

The restrict-split variant can be solved by decomposition, and each of the subsystem problems of this variant can be solved by DP. Thus, solving all possible restrict-split restrictions of the original problem boils down to solving only DP subproblems with p 1 criteria. By definition, these subproblems are restricted to local tasks of some subsystem of the original formulation, in addition to the new local tasks obtained by assigning complex tasks to a particular subsystem according to σ S T C .

Since most of the tasks in the original problem are not complex but local tasks, their treatment will be the same in each node. It is well known that the output of a DP algorithm does not depend on the order in which the decisions are taken. Without loss of generality, we can assume that local tasks are always taken care of first, so that the |T L | first steps of the resolution of decoupled problems are the same in each leaf-node of the decoupling branching scheme. Thus, these steps can be precomputed: each of the subsystem specific subproblem in a leaf-node problem could be kickstarted to the first task step corresponding to a complex tasks having been assigned to this subsystem in the restrict-split variant. But contrary to the kickstarted method considered in Section C.4.3.2, we now solve problems with only one residual capacity criterion.

We recall that T L psq T psq T L , and we consider the feasible set X s L of subproblem s restricted to local task variables, defined as :

X s L : tx t0, 1u |T L psq| | ţT L psq mMpsq x tm ¤ 1 & ţT L psq mMpsq w tm x tm ¤ b s u (4.1)
Algorithm 26 describes the computation of Np f pXqq. First, we compute Q s L ẼpX s L , pf s , ws qq for all s S , i.e. the solutions of the restriction of subproblems to local tasks which are nondominated for f and for ws , the residual capacity criterion local to their subproblem. Then for each σ S T C , we complete the resolution of the |S | independent subproblem with the tasks made local to each of them by decoupling σ, and we compute the non-dominated subset of the cartesian product of these solutions. Since the restrict-split problem has no coupling, this yields all the efficient solutions of this decoupled version of the problem. We then compute the union of all the non-dominated sets obtained in restrict-split variants, and the non-dominated subset of this union is Np f pXqq. As is apparent from results in Table 4.11, the joint use of dynamic programming and the precomputation of the DP subsequence which is common to all restrict-split variants leads to very significant performance improvements. In the bi-objective case, for weakly coupled instances, we solved some instances of the problem 4 to 8 times faster than the generic algorithm, and for p 3, we find 18 to 52-fold speed-ups relative to the generic method.

Finally, we compare the performance of the best methods obtained in Section C.4.3 with the best "top-down" approach. Considering first instances solved in the previous section, we observe, from results presented in Table 4.12 that in all cases, the top-down approach denoted Dec DP PreC performs significantly better than the bottom-up approach denoted DP LB.Filt.

CHAPTER 4. APPLICATION PROBLEM

We further considered "harder" instances, where we increased either the number of local tasks and machines |Mpsq|, or the degree of coupling between the subproblems, with more subsystems and more complex tasks. The superiority of the top-down approach, relative to the bottom-up approach, was confirmed by results obtained on these instances and reported in Table 4.13. When |Mpsq| increases, the top-down approach maintains an advantage over the e ¡ constraint method. However, for p 2, when the degree of coupling increases, it is overcome, and e ¡ constraint remains the quicker method. In the next and last section, we pursue final improvements of the top-down approach aiming at staying competitive with e ¡ constraint. We have seen that in some cases, the full enumeration of S T C is already competitive with a generic method. Thus, we can expect to reach good results by embedding it in a branch and bound approach, where branching is done not over the binary decision variables of the initial formulation, but over additional decision variables restricting a coupling to a particular subsystem. That being said, we remain wary of the two following facts:

1. a tree defined by branching over all possible restrictions of some coupling grows faster in width than a binary tree.

2. a leaf of the tree enumerating the values of decoupling variables does not define a solution because, it merely restrict the possible values of the variables of the original problem. Rather, it defines a decomposable instance of the problem, that is assumed to be solvable quickly.

Single objective approaches to the optimization of complex systems rely on the notion of coordination, i.e. of modifications of the subproblems which allow the independent resolution of these subproblems to yield a feasible solution to the global problem. Decouplings, which correspond to restrict-split variants in the generic form of the integer coupled optimization problem, consist in such coordinations. When solving single objective REF, we may consider the optimal value of an instance, and deduce the decoupled variant from which it is obtained. If for each instance, there is a best decoupling, that which is associated with the optimal value, then we may look for ways to either make a guess about which decoupling to pursue first, or to use branch and bound techniques to eliminate decouplings which we could prove would not lead to the optimal solution.

In the multiobjective case, what would a best decoupling mean? Conversely, would we be able to eliminate some decouplings in hopes to shorten the enumeration of S T C ? Because we will compare sets of non-dominated points obtained by solving various restrict-split instances, let us define the notion of set dominance. Definition 4.1. Let p N, A, B R p . A dominates B if and only if for all y B, there exists some y I A such that y I © y.

Arborescent search and branching

Let T PpNq be a set of subsets of variables indices. These will define, for each level of the BB tree, the subset of variables which will be fixed to some values. In general, T needs not define a partition of N, however it must be the case that tT t r, so that we do not assign values to variables more than once. Let G pV, Eq be a tree of depth |T|, where ν V is called a node. To each node is associated a restriction of the original optimization problem, with X ν the feasible set at node ν. E is such that for all ν, ν I V, if pν, ν I q E then X ν I X ν . For each ν in the tree such that depthpνq t, for each ν I such that pν, ν I q E, and u an evaluation of all decisions variables in t,

X ν I tX PpX ν q | dx X ν , di t, x i u i u
In other words, each child of a given node of depth t has the subvector of variables corresponding to t fixed to some value. Let us now explain how these valuations are defined.

When formalizing a branch and bound resolution method with branching over subsets of variables, one can make use of additional variables subvectors x t for t T representing sets of variables. At each level t T , a branching decision must be made, i.e. an assignment of the variables x i for i t. Borrowing notation from DP (cf. Section C.4.1.2), the set of branching decisions that can be made at level t will be written ∆ t , and τpx, δq will be the solution obtained by modifying

x according to decision δ, i.e. to each δ ∆ t corresponds one assignments of the whole variable subvector x t px i | i tq. Thus, children node of node ν of depth t, as defined by the feasible sets of their associated problems, are The decoupling branching scheme thus defined is not, in the case of REF, such that all variables are fixed in the problems associated with its leaf nodes. This is no problem per se, as long as we have a method to solve efficiently the leaf node instances. Because, in our cases, these instances are restrict-split variants of REF, we can solve them using decomposition, DP and precomputation of DP subsequences for subproblems, as previously described.

Two properties are of particular interest for BB resolution methods, exhaustivity and exclusiv- 3. E is exclusive if and only if for all ν I , ν P such that pν, ν I q E, pν, ν P q E X ν I X ν P r (4.3)

A BB algorithm is correct if the branching scheme is exhaustive. It needs not be exclusive, but exclusivity generally allows for better performance as it avoids redundancy. However, in cases where a particular branching scheme can yield a significantly smaller tree than usual exclusive branching schemes (e.g. branching over binary variable evaluations), with still sufficiently easy node-problems, it may be worth considering a non-exclusive branching scheme. This challenge is undertaken by designing sets of decisions which perform the adequate coordinations.

Fathoming a node in the branch and bound search

Let us introduce the implicit enumeration trick used in BB (without loss of generality we consider maximization). Consider LB R p a set of values of feasible solutions, such that for all y LB, there does not exist any y I LB, y I © y. At any given moment in the development of the method, LB (a weak lower bound) will be the set of the best solution values discovered so far, called the incumbent set. We note S R the subset of R p in which there may still be solutions of interest to be found, i.e. solutions which are not dominated by elements of LB. Formally S R ty R p 2hy I L, y ¨yI u pLB ¡ R p ¡ q, where R p ¡ tx R p | x ¡ 0u. This notion was introduced as the search region by [START_REF] Klamroth | On the representation of the search region in multi-objective optimization[END_REF]. If p 1, |LB| 1, so S R is simply the real half-line anchored in some l R such that LB tlu. If p ¥ 2, S R is defined as the union of subsets of R p called zones. These are anchored in knee points defined by p-uples of incumbent points or partial solutions. They are the smallest set Z of points in R p such that zZ ptzu R p ¡ q S R.

The branch and bound principle is then the following: if one can prove that f pX ν q S R r, one shows that searching Gpνq, the subtree rooted in node ν, cannot lead to the discovery of new non-dominated points. Therefore, ignoring this subtree will not prevent the enumeration of Np f pXqq. This is called a separation problem. [START_REF] Ehrgott | Bounds and bound sets for biobjective combinatorial optimization problems[END_REF], following [START_REF] Sourd | A multiobjective branch-and-bound framework: Application to the biobjective spanning tree problem[END_REF], proposed to prove the absence of interesting solution in the subtree by enclosing f pX ν q in an upper bound set U, thus proving that pUB R p ¡ q S R r. This is illustrated in both single and multiobjective cases by Figure C.14.

For each node ν, let U B ν be an upper bound set to Np f pX ν qq. Then by transitivity of dominance, it is an upper bound to f pX ν q. Since for all ν I with pν, ν I q E, it is the case that X ν I X ν , we also have f pX ν I q f pX ν q. Thus we know that no solution x found in the subtree rooted in ν can be such that f pxq ty R | hy I U B ν , y ¡ y I u : U B R p ¡ . In other words, for Gpνq the subtree of G rooted in ν,

dν I Gpνq, f pX I ν q U B ν R p ¡ r (4.4)
The separation problem has been presented by [START_REF] Sourd | A multiobjective branch-and-bound framework: Application to the biobjective spanning tree problem[END_REF] as that of finding a separating function h : R p Ñ R separating the feasible set as bounded by the upper bound set from the search region, in the sense that hpyq ¥ 0 dy f pX ν q hpyq 0 dy S R (4.5)

A separating function is in general non-convex, reflecting the complexity of an hypersurface separating f pX ν q from the current incumbent set may, when such a surface exists. In practice, the separating function will be computed piece by piece. Each of these pieces defines a cone -which in extreme cases can be a half-space, such that there can be no feasible solutions within that cone.

One way to ensure this condition is to generate these cones together with upper bound points. If these upper bound points are supported by hyperplanes, the pieces will be hyperplanes. Thus the following is a separating function [START_REF] Sourd | A multiobjective branch-and-bound framework: Application to the biobjective spanning tree problem[END_REF]):

h ν Λ pyq min λΛ pmax xX ν xλfpxqy ¡ xλyyq (4.6)
where Λ R p and xλ, yy °p j1 λ j y j . The image of a supported solution is an upper bound point, so any y R p such that h ν Λ pyq 0 is above the hyperplane supporting the optimum of max xX ν xλfpxqy for some λ Λ. In other words it is unfeasible at node ν. Note that it would also be the case for supported solutions of the linear programming relaxation of the problem under consideration at node ν. Non-supported upper bound points still provide a set of cones that can be used for separation, using the following separating function:

h ν pyq 5 1 if hy I U B, y I ¨y ¡1 otherwise (4.7)
The reasoning is even clearer in this case. Given that no solution down in the subtree rooted at ν can dominate an upper bound point y I computed at ν, any point that would dominate y I is to be considered unattainable at ν, in the sense that there exists no solution x X ν such that y I ¨fpxq, or in other words, such that y I tfpxqu R p ¥ . In any case, to state that a point y in unattainable at ν, all we need is one piece of the separating function, either a λ such that max xX ν xλfpxqy ¤ xλyy, or max xX R ν xλfpxqy ¤ xλyy, or one upperbound solution y I such that y © y I , where X R ν is the feasible set of the linear programming relaxation the problem associated with node ν.

The approach to fathoming by bound that is exemplified by both [START_REF] Sourd | A multiobjective branch-and-bound framework: Application to the biobjective spanning tree problem[END_REF] and [START_REF] Parragh | Branch-and-bound for bi-objective integer programming[END_REF] is the following: if a separation function is found between S R and U B, the node is fathomed. The emptiness of S R f pX ν q is guaranteed by the fact that f pX ν q pUB ¡ R p ¥ q. If no separating function is found, branching must occur. [START_REF] Parragh | Branch-and-bound for bi-objective integer programming[END_REF] note that even when separation is not possible, the hyperplanes generated while searching for a separating function provide information restricting the subspace of the objective space where non-dominated points may lie. This leads, in their line of work, to objective space branching.

Our own approach also aims at using information gained by generating segments of a separat-ing function, but it differs from the other approaches in how it defines fathoming. We rely on the notion of search zone elimination that aims at providing two advantages. First, search zone elimination is valid for a whole subtree, and thus any child node will inherit only the zones which were not found to be unattainable at its parent node. Second, the set of upper bound points necessary to show that some zone is unattainable needs not constitute a strong upper bound set.

Given a node ν, consider Z ν the set of zones to be investigated at node ν. We call these zones active at node ν, and will see that oftentimes, S R ν : zZ ν ptzu R p ¡ q is a lot smaller than S R.

As hinted at earlier, in practice, we do not compute the whole separating function at every node. Rather, we generate pieces of the separating function, and for each piece, we try to eliminate as many search zones as possible. We do this until either separation is achieved, or the procedure for generating pieces of the separating function terminates.

For the purpose of our implementation, let us define the function hpc, yq : PpR p q ¢ R p . Where c R p is a piece of the separating function that is a cone: either the upper-right quadrant of a non-supported upper-bound solution, or the halfspace above the supporting hyperplane of an upper bound solution. Note that even when we fail to separate X ν and S R, we may still be able to delete some zones from Z ν . This information is useful, because for any ν I such that pν, ν I q E, X ν I X ν . So, if for some z Z ν , tzu R p ¡ X ν r, then tzu R p ¡ X ν I r. Thus there is no point in further investigating, in Gpνq, a zone that has been proven unattainable at node ν. The zones which remain active after the procedure has been applied are passed on to the children nodes ν I of ν.

in our abilities to fathom nodes by the initial incumbent, and because of the experiment performed in Section 4.4.3.1, we have reasons to believe that we, in fact, only managed to eliminate a small portion of dominated leaves. Further work is thus required, both beyond the biobjective case, and thus using other upper bound notions, and to integrate quicker or more selective update of the search region when feasible solutions are obtained.

Conclusions and perspectives

In this chapter, we have explored a wide array of methods aimed at efficiently solving a discrete, multiobjective complex systems problem we called REF. Starting from the definition of a dynamic programming algorithm to solve REF, we observed that this approach suffered from an increase in the number of sites, corresponding to additional criteria introducing incomparabilities between states of the DP process. We proposed to tackle this issue in two ways.

First "bottom-up", trying to eliminate as much partial solutions as possible while relying only on local information. We used filtering using lower bound computed from the local restriction of the problem, and a kickstarted method in which we obtained the first step of the DP process associated with a complex task by aggregating results of subproblems which were solved using only local information.

Second, we proposed a "top-down" approach, the basic version of which is an exploration of all possible assignments of complex tasks to sites, which are sufficient to decouple an instance of REF. The number of instances to solve in this approach obviously grows quickly with the number of complex tasks and of sites. To compensate for this, we used DP, which efficiently solves instances with p 1 criteria, and we precomputed the final stage of the local part of the problem, used to solve as many very small p 1 criteria DP problems as there are decouplings.

In both cases, we precompute the local part of REF, and we pass it on the a second phase of the computation. The duality between the two approaches is clear: the number of criteria, and thus of tests between incomparable DP states in the bottom-up approach is traded for the number of instances to solve in top-down approach. Our experiments revealed the best top-down approach to perform significantly better than the best bottom-up approach.

We attempted to improve the top-down approach further by presenting the enumeration of decouplings in the form of a branch and bound search, where branching is done on the decision of assigning a complex task to a site. We obtained encouraging but contrasted and limited results in the biobjective case, which invite further research.

General conclusions and perspectives

This work has explored the multiobjective optimization of complex systems. We have broken down the subject into three main topics. First the combination of solutions obtained from the resolution of independent subproblems, then the use of decomposition to provide lower and upper bounds to the non-dominated set of a complex system MO problem, and finally we have developed applications of decompositions for the resolution of a particular complex system optimization problem, named REF.

In Chapter 2, we achieve significant performance improvements in combining and filtering sets of solutions to optimization problems. This operation is a cornerstone of methods we developed further on. We put forth two main approaches to solving this problem. The first attempted to avoid the sorting of the whole set sum, and deduce, from the sorting of individual sets, a dominance preserving order over combinations of solutions, which allows for faster dominance filtering of the set of combinations. Another approach, which did not perform as well as the previous one, involved grouping points in the sets to be combined into boxes. Although we have considered multiple algorithms to produce these boxes -some of which allowed for good results in applying box-based dominance relations, there is room for improvement in this matter. We have managed to provide a box-building algorithm, the parameters of which had to be tuned so as to provide good overall performance, albeit within to confines of a particular number of criteria. Further work should investigate whether such a box-based method could be devised without parameters, or with parameters, the values of which would remain relevant across any number of criteria.

Decomposable relaxations and restrictions of the original problem, presented in Chapter 3, proved to provide a very good approximation of the set of non-dominated solutions of the generic coupled problem. We also showed that the speed up of the computation of the bounds which is allowed by decomposition scales well with both the size of the instances and the number of criteria. This strand of research should be pursued by putting forth more refined notions of lower and upper bounds, that could be adapted to problems in which the coupling between subsystems is more complex than the one we assumed. These would include using partially relaxed or mixed constraints such as surrogate relaxations, or supported solutions which could be computed using existing methods in single objective Lagrangean relaxations. We left untouched the subject of multiobjective Lagrangean duality, because it has been developed mostly in the continuous case, which seemed more foreign to our application cases that methods we could readily develop. However, intersections between this line of inquiry and our own work are to be pursued.

Chapter 4 provided several approaches to integrating decomposable optimization and dynamic programming in a "bottom-up" fashion, both breaking down the DP process into several independent subsequence, and to quickly obtain lower and upper bounds that allow the elimination of partial solutions in large number. For p ¥ 3, the relevance of DP was supported by the current limitations of generic algorithms. Although in most cases we achieved competitiveness for p 2 as well, we observed that as the number of subsystems, and thus the size of DP states, increased, only approaches which used DP for fully decoupled problems could beat the e-constraint method. This motivated so called "Top-down" approaches, where we enumerated fully decoupled variants of the original problem so as to exhaustively cover its feasible set. We initiated a strand of work combining this enumeration of decouplings with multi-objective branch & bound, which gave promising but limited results in the biobjective case, and which we hope to pursue in future research. Within this section, we do not present any strong lower bound concept for the GCP. However, for application problems presented in later sections of this work, we will.

RS) weak lower bounds in purple and orange respectively, and the copy-split (CS) and local omission of coupled variables (LR) upper bounds in light blue and green respectively. Secondly, we plot the randomly sampled points used to quantify the decision maker's uncertainty, conditioned on the computation of the RS and LR bounds only. Dark blue points are those which are either dominated by a point in the RS wLB, or dominate a point in the LR UB, while yellow points are those which fall in the area of uncertainty. Computing RS and CS yields much lower uncertainty, but plotting the resulting area of uncertainty would produce a less visible an thus unclear example. Next, we give further details on the computation which provides this result. avec, pour tout j t1, ..., pu, f j pxq °sS f s j px s q. Si |S | m et px s |s S q px 1 , ..., x m q, on dénote par x s les sous-vecteurs de variables n'appartenant pas à s.

Observation C.1. Si un problème pPq est non couplés, son ensemble réalisable X est tel que X ¹ sS X s (C.1)

C.1.2.2 Optimisation découplée

Si f j pxq est additivement séparable pour tout j t1, ..., pu, la fonction objectif du sous-problème associé à s S est donné par f s px s q pf s 1 px s q, ... f s p px s qq, avec f s j : ± is X i Ñ R pour j t1, ..., pu. Chaque sous-système est associé avec le problème max x s X s f s px s q.

Definition C.7. Pour j t1, ..., pu, f j : ± iN X i Ñ R p , et pour tout s S f s j : ± is X i Ñ R p , f j est additivement separable selon S si, pour tout x px s | s S q avec x s px i | i sq, f j pxq şS f s j px s q Une question d'intérêt pour la généralisation multiobjectif de la décomposition est la suivante, où °sS Y s dénote la somme de Minkowski. Si p 1, on a aussi max xX °sS f s px s q °sS max x s X s f s px s q. Mais pour p ¥ 2, °sS Np f s pX s qq Np f pXqq. Un autre filtrage par dominance est donc nécessaire.

C.1.3 Efficacité de la décomposition

L'efficacité de la décomposition est testée en résolvant un problème générique non couplés, d'abord en utilisant un algorithme générique sur le problème entier, puis en résolvant les sous-problèmes séparément, en combinant les solutions et en filtrant le résultat par dominance. Les résultats de cette expérience suffisent à montrer que la décomposition permet une amélioration significative de la performance pour n'importe quel nombre de critères, et pour une augmentation du nombre de sous-systèmes et du nombre de variables. On propose une heuristique gloutonne pour générer un point bornant fortement toutes les completions d'une solution partielle. Comme dans la relaxation copy-split, les tâches sont copiées, donnant un ensemble T I de tâches. Puis pour T I est remplacé par un ensemble de tâches "idéales"

C.2.2 Approches par boîtes

T ¦ tel que pour tous t T I , m Mpt I q, est l'ensemble des points non dominés par les solutions déjà connues, et est constituée de "zones de recherches" définies par des points charnières situés au minimum des paires de points réalisables déjà connus. En monocritère, cette région de recherche est simplement la demidroite réelle supérieure à la valeur de la meilleure solution connue.

Montrer que les solutions réalisables connues sont meilleures que les solutions d'un sous arbre revient à "séparer" la région de recherche de l'ensemble réalisable de ce sous arbre, ensemble réalisable borné par des ensembles bornant supérieurs, ou des hyperplans définis par des scalarisations, comme illustré dans la Figure C.14.

La séparation est réalisée en éliminant des zones de recherche à l'aide d'hyperplans définissant un demi-espace inaccessible. Si une zone de recherche est contenue dans ce demi espace, elle peut être supprimée du noeud actuel et de ses descendants. Ainsi la séparation n'a pas à être réalisée en un coup.

C.4.4.4 Application à la résolution de REF

On se limite dans cette expérience au cas biobjectif. L'arbre est exploré en profondeur jusqu'aux noeuds feuille, auxquels la variante restrict-split correspondante est résolue par programmation dynamique et décomposition. Aux noeuds intermédiaire, on génère une série d'hyerplans par la méthode dichotomique, pour essayer d'éliminer autant de zones de recherches actives que possible, en stoppant l'exploration à un certain seuil de précision défini par ε-dominance. La région de

  of Multiobjective optimization . . . . . . . . . . . . . . . . . . . . . . . 1.2 Coupled problems and decomposition . . . . . . . . . . . . . . . . . . . . . . . 1.3 The Generic Uncoupled Problem and the efficiency of decomposition . . . . . . 1.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Computing efficiently the non-dominated subset of the Minkowski set sum 2.1 The non-dominated subset of the Minkowski sum problem (NDMSP) . . . . . . 2.2 Intermediary filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 A unidirectional method for pooling . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Box-based methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Conclusions and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Decoupling a coupled problem to obtain bound sets 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Decomposable restrictions and lower bound sets . . . . . . . . . . . . . . . . . . 3.3 Decomposable relaxations and upper bound sets . . . . . . . . . . . . . . . . . . 3.4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Conclusions and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  of Multiobjective optimization . . . . . . . . . . . . . . . . . . . . 1.1.1 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 The specific difficulties of 0 ¡ 1 integer MO . . . . . . . . . . . . . . 1.1.3 Dominance filtering algorithms . . . . . . . . . . . . . . . . . . . . . 1.1.4 Choosing a dominance filtering algorithm . . . . . . . . . . . . . . . . 1.2 Coupled problems and decomposition . . . . . . . . . . . . . . . . . . . . . 1.2.1 Decomposable formulations . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Uncoupled optimization . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 The Generic Uncoupled Problem and the efficiency of decomposition . . . 1.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 1. FUNDAMENTAL NOTIONS

  1.1.3.1 Algorithms over linear list solution archivesThe most basic of dominance filtering algorithms is one which iterates through a list of input points Y, comparing each y Y to every point in r Y, a set of candidate points, i.e. points which have not yet been found to be dominated by any point in Y. For each y I r Y such that y © y I , y I is removed from r Y. If there is no y I r Y such that y I © y, then y is added to r Y. In the worst case, i.e. if Y is an independent set with respect to ¨, this algorithm requires |Y| 2 dominance tests. In practice, move to front-type heuristics as introduced by Bentley et al. (1993) can be used to bring about conclusive dominace tests sooner in the loop over r Y.

Algorithm 7 :

 7 Unidirectional Dominance Filtering for p 2 input : Y R p output: NpYq 1 S ortpY, ¥ Lex q 2 /*The for loop iterates according to ¥ Lex order */ 3 for y Y do 4 /*Assume r Y ty 1 , ..., y m u these two values, as well as the time gain (TG), defined as 100 T.NA¡T.IF T.NA, where T.NA denotes the computing time when filtering once at the end, and T.IF when filtering at each step. SR denotes, in percentage, the reduction in size from the MS to its non-dominated subset. For example, in a trial with p 2, |S | 3, and |Y s | 100, the MS has a size | °sS Y s | 10 6 , whereas the average size of the non-dominated set, denoted |ND|, is 561.9, corresponding to a size reduction of 99.94%. Values between squared brackets correspond to standard deviations.

  Figure C.3 exhibits table T and Figure C.2 the parenthood

19q Figure 2

 19q2 Figure 2.1: Lattice-representation of the parenthood relation over Y Z in Example C.1

2 Figure 2

 22 Figure 2.2: Table representation of the parenthood relation over Y Z in Example C.1

  Assume that pi, jq BoxCombinations is the combination of B Y a box of elements of Y, and B Z a box of elements of elements of Z, where Y Ð Y ÝÝÑ s ¡ 1 and Z Ð Y s . that NDBoxComb tpB 1

Algorithm 15 :

 15 Boxing using ε-dominance input : Y R p , ε R output: pB 1 , ..., B m q 1 Boxes Ð r 2 while Y r do 3 The goal of the latter measurement is to contrast the ways in which box-to-box and point-tobox dominance filtering yield time gain. We consider box-to-box and point-to-box independently. Instances are identical to those used in the previous experiment.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.1.1 Bound sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.1.2 Multiobjective relaxations and restrictions . . . . . . . . . . . . . . . . 61 3.2 Decomposable restrictions and lower bound sets . . . . . . . . . . . . . . . 63 3.2.1 Admissible variable values, neutral variable values and restrictions . . 63 3.2.2 Restrict-splitting coupling constraints . . . . . . . . . . . . . . . . . . 64 3.3 Decomposable relaxations and upper bound sets . . . . . . . . . . . . . . . 65 3.3.1 Copy-splitting coupling constraints . . . . . . . . . . . . . . . . . . . 66 3.3.2 Local relaxation of coupled variables in subproblem constraints . . . . 66 3.4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.5 Conclusions and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

  Definition 3.4.[START_REF] Ehrgott | Bounds and bound sets for biobjective combinatorial optimization problems[END_REF]) For any Y R p , LB is a strong lower bound set to Y (in the maximization case) if and only if 5 dy Y hy I LB, y © y I dy LB, y is a lower bound point to Y (3.2) or equivalently if Y LB R p © . Example 3.1. Figure C.5 provides an example of a set N tn i | i 1, ..., 5u which is weakly bounded from below by set L tl i | i 1, ..., 5u, because n 3 dominates no element in L. However U tu i | i 1, ..., 5u is a strong upper bound on N.

Figure 3 . 2 :

 32 Figure 3.2: Graph representation of the modification yielding the restrict split variant.

Figure 3 . 3 :

 33 Figure 3.3: Graph representation of the modification yielding the copy split variant.

Figure 3 . 4 :

 34 Figure 3.4: Graph representation of the modification yielding the local relaxation variant.

  in coupling constraints) Ks ds S (uncoupled variables local to subsystem s)

  Figure C.8 illustrates the uncoupled variant of the original problem obtained by removing terms associated with coupled variables from non-coupling constraints in subsystem local constraints.

  constraint k K, where |s| is the number of variables involved in subsystem s.

Figure 3 . 5 :

 35 Figure 3.5: Plot representation of the original non-dominated set and bound sets. Instance parameters are p 3, n 75, |S | 2, |K| 1, DV = 10

  3.3.1) of the problem, which provides an upper bound set. LR denotes the relaxation obtained by omitting coupled variables in local subproblem constraints (Section C.3.3.2), which also provides an upper bound set. T. denotes computing time, |.| denotes the size of a set of points, and ApQ.pεq denotes the quality of approximation in terms of a posteriori ε dominance.

CHAPTER

Figure 4 . 2 :

 42 Figure 4.2: Assignment graphs with knapsack constraints for REF (left) and GAP (right). REF has one knapsack constraint per site, while GAP has one knapsack constraint per machine.

CHAPTER 4 .

 4 APPLICATION PROBLEM Algorithm 16: Dynamic Programming Algorithm for REF input

  reported in Table 4.1, we solve instances of REF for various values of p and instance parameters |Mpsq|, |T L psq|, |T C |, and |S |.

  mMpSptqq

Figure 4 . 4 :

 44 Figure 4.4: Visualization of the copy-split relaxation as applied to REF: assignments of complex task t 1 are duplicated and made simple tasks assignments in each subsystem

  mMpsq

4 S

 4 ubRes s Ð DPpte N u, r, T I L psq, Mpsq, tsuq 5 return PoolpS ubResq 4.2.1.2 Relaxation of complex task assignment variables from non-coupling constraints

Figure 4 . 5 :

 45 Figure 4.5: Profits and weights associated with T ¦ when applied to Example 4.3

CHAPTER 4 .

 4 APPLICATION PROBLEM Algorithm 21: Greedy wLB from completions of solutions to LocRes input : B X L , Asg 1 /*Asg denotes the list of ranked assignment variable indices */ output: κpBq 2 OutS et tu 3 for x B do 4 while x X & Asg r do 5 for pt, mq Asg do 6 x Ð xr xtm Ð 1s 7 Asg Ð Asgztpt I , mq Asg | t I tu /*Remove all assignments of task t Asgztpt, mq Asg | m Mpsqu /*If capacity of site s is exceeded, remove all assignments of remaining tasks to machines of site s.

Algorithm 23 :

 23 Partial filtering in the DP formulation of REF input : pT C , T L , M, S q output: ẼpX, f q 1 for t t1, ...,

4. 4 . 3

 43 Towards a branch and bound method for solving coupled problems4.4.3.1 Dominance between decouplings

Figure 4 .

 4 Figure 4.8: Binary branching (top) versus decoupling branching (bottom)

Figure 4 . 9 :

 49 Figure 4.9: The notions of search region, search zones and separation illustrated in single and biobjective cases. In case (b), search zones are defined by knee points n 1 to n 5 , deduced from the lower bound tl 1 , ..., l 4 u and some objective-wise upper bounds

  the process of generating cone-like pieces of the separating function and using them to prove that S R X ν r.

Algorithm 27 :

 27 Separation Algorithm (Separate) input : ν a node, C ν a set of cones generated as node ν output: 1 if the separation succeeds, 0 otherwise1 Z ν Ð Z ν ¦ where ν ¦ is the father node of ν in G 2

Figure A. 1 :

 1 Figure A.1: An example of a double bound and a sample of points drawn randomly to measure area of uncertainty after computing the RS wLB and the LRUB

Figure A. 1

 1 Figure A.1 features a graphical representation of the local restriction (LocR) and restrict-split (RS) weak lower bounds in purple and orange respectively, and the copy-split (CS) and local omission of coupled variables (LR) upper bounds in light blue and green respectively. Secondly, we plot the randomly sampled points used to quantify the decision maker's uncertainty, conditioned on the computation of the RS and LR bounds only. Dark blue points are those which are either dominated by a point in the RS wLB, or dominate a point in the LR UB, while yellow points are those which fall in the area of uncertainty. Computing RS and CS yields much lower uncertainty, but plotting the resulting area of uncertainty would produce a less visible an thus unclear example. Next, we give further details on the computation which provides this result.

  Figure C.1: Graph representation of a coupled system

  EpX, f q ? ¹ sS pEpX s , f s qq ou, dans l'espace des objectifs Np f pXqq ? sS N pf s pX s qq Proposition C.1. Si ± sS X s X et si pour tout j t1, ..., pu, f j est séparable selon S , alors Np f pXqq sS Np f s pX s qq

C. 2

 2 Calcul efficace du sous ensemble non-dominé d'une somme d'ensembles C.2.1 Définition du problème Il s'agit de calculer efficacement Np °sS Y s q, avec Y s R p pour tout s S . Le défi principal de ce calcul est d'organiser correctement la sommation et le filtrage afin de générer le point de combinaisons dominées possible. C.2.1.1 Filtrage intermédiaire Un algorithme séquentiel considère les ensembles Y s l'un après l'autre, et à chaque étape s S , génère les combinaisons entre points de Y s et points de Y ÝÝÑ s ¡ 1 . Plutôt que de générer la somme pour les tous les s S puis filtrer, le résultat suivant assure qu'il est correct d'appliquer un filtrage par dominance à chaque étape. Proposition C.2. Pour toute famille pY s | s S q avec Y s fini pour chaque s S et S t1, ..., nu, pour tout s I ¤ n, peuvent être comparés du point de vue de la complexité. Supposons que |Y| m, |Z| n, with m ¥ n. IF et UPool diffèrent dans la phase de tri préalable au filtrage par algorithme unidirectionnel. IF nécessite de trier Y Z, en temps Opmn log mnq Opmn log mq. UPool nécessite de trier Y et Z indépendamment, en temps Opm log m n log nq Opm log mq, puis d'insérer les mn elements dans H. Il est possible de prouver que H est toujours de taille inférieure à n. L'opération est donc réalisée en temps Opmn log nq. Comme dans le cas biobjectif, le tri représente la majorité de temps de calcul (l'insertion étant en temps constant) UPool a de bonnes chances d'être plus performant que IF dans ce cas, où UPool a une complexité de Opmn log nq et IF de Opmn log mq. C.2.1.5 Résultats expérimentaux Pour p 2 on trouve expérimentalement que UPool offre des performances jusqu'à 81, 65 % supérieures à IF, augmentant avec une augmentation de |Y s |. De même pour p 3, mais seulement jusqu'à 49.33%. Pour p 4, 5, UPool n'offre pas de meilleures performances que IF.

  Figure C.10: Le graphe p1q représentes toutes les affecations possibles ainsi que les groupements de machines en sites. Le graphe p2q montre une solution réalisable, et les graphes p3q et p4q des solutions irréalisables

  mMpSptqq

  Figure C.13: Binary branching (top) versus decoupling branching (bottom)

  The for loop iterates according to © Lex order

	1.1.3.2 Algorithms over binary trees with p-dimensional intermediary node data (KDTrees)
			*/
	4 for y Y do 5 Dominated Ð False
	6	if p 2 then
	10 11 12 13 14 15	else /*p ¥ 3 for y I r Y do if y I © y then Y Ð moveT oFrontpY, y I q Dominated Ð T rue	*/
	16 17 18	if Dominated False then r Y Ð Y tyu /*If p 2, this means an insertion at the back. If p ¥ 2,
		depending on the heuristic used, this could mean an insertion at the front, or at the back of r Y.	*/

7 /*Assume r Y ty 1 , ..., y m u */ 8 if y m © y then 9 Dominated Ð T rue

  y Yzty 1 u do

	6 7 8	if s y then if ||s|| 2 ||y|| 2 then s Ð y
	11	if k	|Y| δ	then
	12		/*when	|Y| δ	input points have been considered	*/

9

if DominatedpT ree 1 , yq False then 10 InsertpT ree 1 , yq 13 PruningpT ree 1 q 14 k Ð k 1 15 Initialize T ree 2 16 Assume ContentpT ree 1 q ty 1 , ..., y m u 17 for y ty m , ..., y 1 u do 18 if DominatedpT ree 2 , yq False then 19

mu do 3

  Each line of Table1.1 reports parameter values and average measurements for 10 trials of the experiment. The number of criteria is denoted by p, n stands for the number of variables, |S | denotes the number of subproblems, and av. s the average number, over trials, of variables per subproblem. |E| denotes the number of efficient solutions of the problem. Then, av. |E s | : p °sS |E s |q{|S | denotes the average number of non-dominated points over subproblems, and av. T. E s denotes the average time, over subproblems, for solving a subproblem. The T. No Dec column records the time spent computing the non-dominated set for the whole problem with the e-constraint method. Finally, T.Dec denotes the time spent solving the same problem using decomposition as previously described, and T. Pool is the time spent combining efficient solutions of the subproblems, and filtering it for dominance over the p original criteria.

	p	n	|S| av |s|	|E|		av |E s |	T. NoDec	T. Dec	av T.E s	T.Pool
	2	100	2	50	134.8	[18.4] 47.25	[5.63] 10.19	[2.25]	4.71	[0.77]	2.353 [0.385] 0.005 [0.001]
	2	100	3	33.3	124.8	[31.5] 23.91	[4.96]	9.92	[3.55]	3.06	[0.92]	1.011 [0.303] 0.026 [0.014]
	2	100	4	25	106.1	[29.0] 14.13	[3.50]	8.81	[3.03]	2.15	[0.72]	0.505 [0.154] 0.125 [0.124]
	2	120	2	60	199.4	[58.4] 67.45 [18.59] 22.74 [18.49]	8.37	[3.34]	4.181 [1.664] 0.011 [0.008]
	2	120	2	60	191.3	[54.1] 65.15 [17.58] 28.64 [15.39] 11.89	[5.84]	5.940 [2.920] 0.011 [0.006]
	2	120	4	30	160.3	[30.0] 19.33	[2.84] 15.69	[6.01]	3.12	[0.79]	0.693 [0.146] 0.351 [0.241]
	3	30	2	15	65.9	[24.9] 18.22	[6.34]	4.52	[1.94]	1.95	[0.63]	0.971 [0.315] 0.001 [0.001]
	3	30	3	10	59.9	[27.3]	8.67	[3.22]	3.96	[2.37]	1.63	[0.39]	0.544 [0.131] 0.001 [0.001]
	3	30	4	7.5	66.2	[37.3]	5.52	[1.28]	3.82	[2.12]	1.80	[0.07]	0.449 [0.016] 0.002 [0.001]
	3	40	2	20	218.2	[80.4] 44.75 [14.85] 22.55 [11.18]	5.03	[2.22]	2.509 [1.105] 0.008 [0.005]
	3	40	2	20	202.6 [119.7] 40.50 [20.77] 19.63 [14.38]	4.37	[2.57]	2.178 [1.281] 0.008 [0.009]
	3	40	4	10	180.1 [113.8] 10.05	[4.42] 15.81 [11.49]	2.49	[0.80]	0.619 [0.199] 0.013 [0.010]
	4	20	2	10	70.5	[37.5] 18.65	[8.48]	6.39	[4.14]	2.55	[1.25]	1.273 [0.626] 0.001 [0.001]
	4	20	3	6.7	73.1	[60.6]	6.24	[2.27]	6.13	[5.72]	1.51	[0.19]	0.503 [0.064] 0.002 [0.002]
	4	20	4	5	37.9	[19.5]	3.71	[0.76]	2.50	[1.36]	1.73	[0.08]	0.432 [0.020] 0.001 [0.001]
	4	25	2	12.5	230.1 [196.2] 35.80 [25.08] 36.07 [41.15]	5.50	[4.79]	2.744 [2.391] 0.007 [0.008]
	4	25	2	12.5	217.9	[98.3] 29.95	[9.65] 30.92 [18.84]	3.86	[1.86]	1.927 [0.927] 0.006 [0.004]
	4	25	4	6.25	84.6	[44.7]	5.38	[1.40]	8.05	[4.59]	1.98	[0.23]	0.494 [0.059] 0.002 [0.001]
	Table 1.1: Running time (in seconds) for solving an uncoupled problem without, or with decom-
	position (average values for 10 trials).
	4 5	Pooled : r for e E ÝÝÑ s ¡ 1 do for e I E s do
	6 7 8 9			/*ē is a new solution object solpēq Ð psolpeq, solpe I qq valpēq Ð valpeq valpe I q Pooled Ð Pooled tēu					*/
	10 11 return EpE Ý Ñ m q E Ý Ñ s Ð Pooled

Table 2 .

 2 2 reports step-wise measurements for the IF version, on the same instances. For s S ,

  CHAPTER 2. COMPUTING EFFICIENTLY THE NON-DOMINATED SUBSET OF THE MINKOWSKI SET SUM

	p |S| |Y s |	|ND|	S R p%q		NA T ime	T ime	IF	TG p%q
	3	100	552.60 [125.96]	99.95	[0.01]	3.17	[0.08]	0.06	[0.01]	98.22	[0.31]
	3	100 1523.82 [295.30]	99.85	[0.03]	4.28	[0.62]	0.41	[0.11]	90.54	[2.16]
	3	100 5265.10 [1038.13]	99.47	[0.10]	7.27	[0.93]	1.17	[0.28]	83.76	[3.96]
	3	100 14798.60 [3610.59]	98.52	[0.36]	12.55	[1.70]	3.48	[0.86]	72.42	[4.87]
	3	125	654.18 [153.81]	99.97	[0.01]	10.11	[0.24]	0.16	[0.02]	98.42	[0.23]
	3	125 1821.82 [367.57]	99.91	[0.02]	14.53	[1.36]	0.74	[0.13]	94.92	[0.92]
	3	125 7752.73 [1595.51]	99.60	[0.08]	21.40	[2.16]	2.40	[0.59]	88.81	[2.42]
	3	125 21380.10 [4460.15]	98.91	[0.23]	32.81	[2.92]	6.73	[1.45]	79.55	[3.45]
	4	60	544.08 [150.59] 100.00	[0.00] 165.89	[26.25]	0.18	[0.06]	99.89	[0.03]
	4	60 2486.45 [584.10]	99.98	[0.01] 277.19	[11.55]	1.14	[0.25]	99.59	[0.09]
	4	60 9846.20 [2872.89]	99.92	[0.02] 309.14	[13.56]	3.82	[1.13]	98.77	[0.33]
	4	60 35223.10 [9897.77]	99.73	[0.08] 435.65	[63.42]	15.65	[4.92]	96.39	[1.06]

Table 2

 2

				IF T ime			IF S ize Reduction
			Step 1	Step 2	Step 3	Step 0	Step 1	Step 2
	3	100	0.012 [0.002] 0.037 [0.009]		97.178 [0.540] 98.020 [0.414]
	3	100	0.045 [0.009] 0.351 [0.101]		94.604 [0.914] 97.164 [0.349]
	3	100	0.078 [0.013] 1.065 [0.271]		89.014 [2.025] 95.177 [0.674]
	3	100	0.136 [0.026] 3.219 [0.817]		80.228 [3.030] 92.563 [1.166]
	3	125	0.028 [0.002] 0.088 [0.018]		97.814 [0.358] 98.469 [0.247]
	3	125	0.083 [0.013] 0.573 [0.124]		95.902 [0.619] 97.724 [0.284]
	3	125	0.143 [0.018] 2.053 [0.546]		90.302 [1.726] 95.898 [0.495]
	3	125	0.233 [0.027] 6.020 [1.370]		83.333 [2.460] 93.432 [0.915]
	4	60	0.007 [0.003] 0.034 [0.012] 0.071 [0.026] 95.215 [1.251] 96.646 [0.574] 97.341 [0.333]
	4	60	0.025 [0.005] 0.194 [0.048] 0.704 [0.164] 90.114 [1.812] 94.843 [0.759] 96.210 [0.418]
	4	60	0.032 [0.008] 0.500 [0.162] 2.673 [0.849] 82.689 [3.257] 91.881 [1.127] 94.589 [0.871]
	4	60	0.045 [0.009] 1.312 [0.297] 11.864 [3.998] 72.479 [4.155] 87.243 [1.571] 92.334 [1.193]

.1: Global performance improvement due to filtering at each step, as a function of the number of objectives, the number of sets and the number of points in each set. Average values for trials. p |S| |Y s |

Table 2

 2

.2:

  CHAPTER 2. COMPUTING EFFICIENTLY THE NON-DOMINATED SUBSET OF THE MINKOWSKI SET SUM

					S qIF	T ime	PW Pool T ime Gain (%)
	2	2	400	1133.3 [275.1]	1.376 [0.172]	1.373 [0.188]	0.292 [4.254]
	3	2	400	2055.2 [361.7]	2.86 [0.645]	2.794 [0.574]	1.874 [3.904]
	4	2	400	4712.4 [536.8]	10.879 [2.85]	10.809 [2.691] -0.785 [15.921]
	2	3	100	550.7 [120.5]	0.248 [0.058]	0.244 [0.053]	1.519 [3.647]
	3	3	100	1540.5 [264.4]	1.371 [0.359]	1.357 [0.354]	0.91 [2.132]
	4	3	100 5859.8 [1005.5] 14.218 [5.493]	14.02 [5.24]	0.824 [3.511]
	2	4	40	329.8 [77.6]	0.078 [0.014]	0.082 [0.018]	-5.065 [14.5]
	3	4	40	1589 [446.6]	1.431 [0.681]	1.381 [0.634]	1.897 [10.61]
	4	4	40	6000.9 [1360.4] 15.401 [6.486] 15.157 [6.277]	1.308 [4.646]
	2	5	20	280.4 [85.9]	0.042 [0.011]	0.043 [0.01]	-2.726 [11.01]
	3	5	20	1384.1 [216.1]	0.971 [0.302]	0.965 [0.287]	0.247 [2.989]
	4	5	20	5295.7 [1155.9] 12.898 [5.379] 12.776 [5.255]	0.312 [7.079]

p |S| |Y s | |ND|

Table 2 .

 2 3: Time Gain for Pairwise Pooling relative to Sequential Pooling with filtering at each step. Average values for 20 trials.

  |Y|u do In the following, using the same sets Y and Z and notations as in Example C.1, we develop an execution of the algorithm for illustrative purposes. Let H represent the heap, and P represent the set NpY Zq being constructed.

	22 return Y	Ý Ñ |S |
	Example 2.2 shows an execution of Algorithm 9.
	Example 2.2. 0. We begin with P t p13, 7qu and H t p1,1q p12, 8q, p1,2q	p2,1q

14 for j t2, ..., |Z| do 15 T pi, jq Ð 2 16 while H r do 17 h Ð rootpHq 18 if notDominatedIn(h, P) then 19 P Ð P thu /*For p 2, insertion is at the back. For p ¥ 3, it may be at the front or at the back depending on the chosen heuristic */ 20 updateph, H, T, Y, Zq 21 Y Ý Ñ s Ð P

  11qu. p8, 12q is dominated by p8, 14q, so that P remains unchanged. p2, 4q is the only child of p1, 4q, thus T p2, 4q Ð 2 ¡ 1 1 11qu. p8, 12q is dominated by p8, 14q, so that P remains unchanged. T p4, 2q Ð 1 ¡ 1 0, T p3, 3q Ð 2 ¡ 1 1 15qu. p8, 11q is dominated by p8, 14q, so that P remains unchanged. T p3, 3q Ð 1 ¡ 1 0, T p2, 4q Ð 1 ¡ 1 0. 14qu. p7, 15q is not dominated in P, so that P Ð P tp7, 15qu. T p4, 3q Ð 2 ¡ 1 1. 14qu. p7, 13q is dominated by p7, 15q, so that P remains unchanged. T p4, 3q Ð 1 ¡ 1 0, T p3, 4q Ð 2 ¡ 1 1 14qu. p6, 16q is not dominated in P, so that P Ð P tp6, 16qu. T p4, 4q Ð 2 ¡ 1 1. 14qu. p5, 14q is dominated by p6, 16q, so that P is unchanged. T p3, 4q Ð 1¡1 0. 16qu. p4, 16q is dominated by p6, 16q, so that P is unchanged. T p4, 4q Ð 1¡1 0. 19qu. p3, 19q is not dominated in P so that P Ð tp13, 7q, p12, 8q, p9, 11q, p8, 14q, p7, 15q, p6, 16q, p3, 19qu 16. H r, end.

	7. H t p8, 8. H t p1,4q p8, 12q, p3,2q p8, 12q, p2,3q p3,2q p8, 12q, p2,3q p8, 9. H t p2,3q p8, 11q, p4,2q p7, CHAPTER 2. COMPUTING EFFICIENTLY THE NON-DOMINATED SUBSET OF THE MINKOWSKI SET SUM 10. H t p4,2q p7, 15q, p3,3q p7, 13q, p2,4q p5, 11. H t p3,3q p7, 13q, p2,4q p5, 12. H t p4,3q p6, 16q, p2,4q p5, 13. H t p2,4q p5, 14. H t p3,4q p4, 15. H t p4,4q p3,

  Table 2.4.

			IF UF	UPool UF	IF CH	UPool CH
			T ime	T ime	TG{UF(%)	T ime	T ime	TG{CH(%)
	2	200	0.197 [0.05]	0.072 [0.02]	63.322 [8.87]	0.383 [0.11]	0.086 [0.02]	77.668 [5.67]
	2	400	0.466 [0.12]	0.132 [0.03]	71.675 [6.70]	0.717 [0.14]	0.171 [0.04]	76.100 [6.61]
	2	600	0.923 [0.06]	0.230 [0.02]	75.085 [1.27]	1.305 [0.25]	0.334 [0.02]	74.372 [4.97]
	3	200	0.314 [0.06]	0.099 [0.02]	68.635 [4.21]	0.805 [0.27]	0.121 [0.02]	85.014 [4.42]
	3	400	1.294 [0.30]	0.346 [0.07]	73.268 [2.38]	3.171 [0.83]	0.485 [0.11]	84.712 [3.23]
	3	600	4.357 [1.59]	0.861 [0.17]	80.247 [6.30]	7.467 [1.64]	1.342 [0.22]	82.030 [3.68]
	4	200	0.952 [0.38]	0.241 [0.05]	74.690 [4.99]	2.764 [1.00]	0.316 [0.08]	88.573 [2.66]
	4	400	5.178 [2.16]	1.056 [0.39]	79.604 [4.17]	13.173 [6.09]	1.580 [0.64]	88.008 [2.49]
	4 2 2 2 3 3 3 4 4 4	600 200 300 400 200 300 400 200 300 400	11.534 [3.19] 0.158 [0.03] 0.339 [0.04] 0.668 [0.09] 1.366 [0.51] 2.902 [0.85] 5.886 [1.77] 7.049 [3.11] 15.396 [4.78] 29.178 [12.22]	2.117 [0.62] 0.115 [0.02] 0.220 [0.05] 0.393 [0.08] 0.906 [0.39] 1.546 [0.54] 3.013 [1.23] 4.522 [2.43] 8.342 [3.68] 14.785 [9.27]	81.648 [2.27] 26.872 [7.92] 34.976 [7.04] 41.138 [6.90] 33.664 [5.96] 46.716 [4.69] 48.815 [6.19] 35.852 [4.64] 45.813 [5.15] 49.328 [6.68]	30.086 [14.25] 0.192 [0.04] 0.400 [0.07] 0.715 [0.11] 1.736 [0.47] 3.731 [1.04] 7.048 [1.48] 8.131 [2.51] 17.755 [4.01] 34.559 [10.53]	3.518 [1.21] 0.387 [0.06] ¡101.879 [39.44] 88.308 [2.91] 0.859 [0.13] ¡114.893 [39.70] 1.774 [0.42] ¡148.161 [53.44] 3.393 [1.03] ¡95.447 [41.16] 7.536 [2.18] ¡101.997 [35.78] 14.754 [4.72] ¡109.355 [43.69] 16.123 [6.80] ¡98.289 [39.88] 37.823 [10.29] ¡113.029 [27.41] 67.384 [25.11] ¡94.984 [30.92]

p |S| |Y s |

Table 2 .

 2

			IF UF	UPool UF	IF CH	UPool CH
	2 2 2 3 3 3 4 4 4 2 2 2	100 150 200 100 150 200 100 150 200 50 100 150	T ime 0.089 [0.02] 0.249 [0.06] 0.436 [0.09] 2.247 [0.94] 6.420 [2.20] 10.247 [4.41] 20.134 [9.85] 92.005 [36.06] 149.676 [128.74] 0.046 [0.01] 0.277 [0.07] 0.784 [0.19]	T ime 0.109 [0.02] 0.268 [0.07] 0.440 [0.09] 2.034 [0.88] 5.771 [2.04] 9.316 [4.11] 19.295 [9.65] 89.073 [37.87] 138.579 [129.90] 0.099 [0.02] ¡113.898 [20.57] TG{UF(%) ¡22.409 [11.36] ¡7.712 [8.90] ¡0.976 [7.68] 9.484 [5.21] 10.115 [5.08] 9.087 [2.46] 4.170 [2.30] 3.188 [5.69] 7.414 [3.94] 0.358 [0.08] ¡29.349 [11.60] 0.861 [0.18] ¡9.923 [9.57]	T ime 0.082 [0.02] 0.184 [0.03] 0.320 [0.05] 1.206 [0.33] 2.866 [0.61] 4.733 [1.12] 6.324 [1.93] 25.035 [5.68] 42.377 [15.59] 0.029 [0.00] 0.142 [0.02] 0.327 [0.04]	T ime 0.097 [0.01] 0.215 [0.02] 0.368 [0.05] 1.229 [0.30] 3.239 [0.65] 5.850 [1.29] 7.752 [2.11] 27.216 [6.02] 43.912 [15.47] 0.073 [0.01] ¡154.276 [38.12] TG{CH(%) ¡17.112 [18.30] ¡17.053 [22.79] ¡15.029 [27.05] ¡1.953 [12.71] ¡13.003 [10.57] ¡23.590 [13.96] ¡22.588 [16.89] ¡8.713 [14.92] ¡3.622 [10.94] 0.173 [0.01] ¡22.214 [17.77] 0.307 [0.03] 6.249 [16.69]
	3	50	2.363 [0.80]	2.287 [0.79]	3.229 [5.90]	0.786 [0.20]	0.572 [0.09]	27.146 [17.02]
	3 3 4 4 4	100 150 50 100 150	14.074 [6.50] 55.894 [25.22] 72.060 [59.94] ¡ 500 ¡ 500	13.593 [6.26] 57.550 [27.14] 77.700 [71.43] ¡ 500 ¡ 500	3.413 [3.15] ¡2.963 [2.08] ¡7.827 [14.81] n.a. n.a.	3.564 [1.01] 8.430 [1.68] 7.346 [4.79] 35.241 [8.89] 111.945 [34.36]	2.657 [0.59] 7.642 [1.39] 6.175 [3.16] 35.016 [9.33] 118.006 [27.95]	25.442 [11.07] 9.348 [8.24] 15.942 [16.57] 0.638 [13.80] ¡5.414 [17.37]

4: Computing time for the sequential method with intermediary filtering, the UPool method, and relative time gain as a function of |Y s | and |S |, for p 2, 3. Average values for trials. p |S| |Y s |

Table 2 .

 2 5: Computing time for the sequential method with intermediary filtering, the UPool method, and relative time gain as a function of |Y s | and |S |, for p 4, 5. Average values for

	trials.

  Λ Ð generateWeightSetpn λ q 2 P Ð getS upportedpY ÝÝÑ s ¡ 1 , Y s , Λq 3 for pi, jq BoxCombinations do

	CHAPTER 2. COMPUTING EFFICIENTLY THE NON-DOMINATED SUBSET OF THE
	MINKOWSKI SET SUM
	2.4.2.2 Box-to-box dominance
	7	upported then S upported Ð S upported ty ¦ u
	8 return S upported
	Algorithm 13: PointFiltBox
		input : BoxCombinations, n λ
		output: BoxCombinations
	4 5 6	for ȳ P do if ȳ © u Y i ÝÝÑ s ¡ 1 BoxCombinations Ð BoxCombinationszpi, jq u Y s j then
	7 return BoxCombinations
	50

1

  Table 2.6), box-based filtering with box-to-box dominance performed worse than the best version of UPool in this case (column TG{UPUF), by a large margin. However, for p 4 (Table 2.7), where UPool tended to perform rather worse than IF, CHAPTER 2. COMPUTING EFFICIENTLY THE NON-DOMINATED SUBSET OF THE MINKOWSKI SET SUM

			IF UF		Box ¡ to ¡ Box
	2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3	2 2 2 3 3 3 4 4 4 2 2 2 3 3 3 4 4 4	T ime 0.197 [0.05] 0.466 [0.12] 0.923 [0.06] 0.314 [0.06] 1.294 [0.30] 4.357 [1.59] 0.952 [0.38] 5.178 [2.16] 11.534 [3.19] 0.158 [0.03] 0.339 [0.04] 0.668 [0.09] 1.366 [0.51] 2.902 [0.85] 5.886 [1.77] 7.049 [3.11] 15.396 [4.78] 29.178 [12.22]	T ime 0.117 [0.03] 0.184 [0.06] 0.362 [0.06] 0.166 [0.04] 0.620 [0.34] 1.802 [0.56] 0.562 [0.22] 2.850 [1.94] 6.700 [3.02] 0.163 [0.03] 0.305 [0.06] 0.498 [0.13] 1.339 [0.51] 2.196 [0.77] 4.046 [1.82] 6.870 [3.42] 12.854 [5.87] 23.715 [14.31]	TG{UF 40.513 [10.79] 60.586 [14.60] 60.737 [7.26] 47.056 [8.93] 52.112 [14.42] 58.649 [11.43] 40.980 [13.07] 44.959 [23.14] 41.914 [12.05] ¡3.141 [9.95] 9.873 [10.32] 25.496 [11.08] 1.922 [12.48] 24.325 [9.91] 31.261 [11.71] 2.534 [9.91] 16.509 [14.76] 18.723 [14.29]	av.Av 1.686 [0.10] 3.040 [0.17] 4.572 [0.36] 2.348 [0.18] 4.322 [0.45] 6.770 [0.80] ¡109.336 [46.55] TG{UPUF (%) ¡62.188 [23.45] ¡39.152 [44.54] ¡57.590 [27.38] ¡68.798 [25.70] ¡79.141 [51.54] 3.183 [0.39] ¡133.188 [41.00] 5.972 [1.00] ¡169.858 [95.96] 8.716 [1.07] ¡216.514 [63.98] 2.244 [0.15] ¡41.041 [11.24] 3.054 [0.22] ¡38.606 [12.43] 3.874 [0.28] ¡26.574 [12.87] 5.123 [0.69] ¡47.850 [23.87] 7.363 [1.06] ¡42.022 [13.67] 9.381 [1.03] ¡34.296 [13.83] 9.464 [2.02] ¡51.938 [16.54] 13.467 [2.46] ¡54.079 [17.40] 15.829 [3.62] ¡60.398 [15.73]

p |S| |Y s |

Table 2

 2

			IFCH		Box ¡ to ¡ Box
			T ime	T ime	TG{CH	av.Av	TG{UPCH (%)
	4	2	0.082 [0.02]	0.081 [0.01]	0.081 [11.74]	2.346 [0.13]	16.477 [11.25]
	4	2	0.184 [0.03]	0.174 [0.03]	0.174 [14.12]	2.803 [0.25]	19.223 [13.30]
	4	2	0.320 [0.05]	0.293 [0.09]	0.293 [19.14]	3.404 [0.43]	20.426 [16.97]
	4	3	1.206 [0.33]	1.214 [0.37]	1.214 [17.04]	3.490 [0.67]	1.273 [19.42]
	4	3	2.866 [0.61]	2.746 [0.62]	2.746 [11.10]	3.881 [0.97]	15.222 [8.97]
	4	3	4.733 [1.12]	4.667 [1.17]	4.667 [12.20]	3.846 [1.27]	20.221 [10.24]
	4	4	6.324 [1.93]	6.384 [1.80]	6.384 [10.63]	2.328 [0.46]	17.646 [11.06]
	4	4	25.035 [5.68]	24.228 [5.26]	24.228 [11.74]	2.458 [0.65]	10.981 [11.05]
	4	4	42.377 [15.59]	41.476 [16.14]	41.476 [6.96]	2.979 [1.18]	5.548 [11.44]

.6: Computing time for the sequential method with intermediary filtering, for box-to-box dominance, and relative time gain as a function of |Y s |, for p 2, 3. Average values for 20 trials. p |S| |Y s |

Table 2 .

 2 7: Computing time for the sequential method with intermediary filtering, for box-to-box dominance, and relative time gain as a function of |Y s |, for p 4. Average values for 20 trials.

	|S| |Y s | 2 200 2 400 2 600 3 200 3 400 3 600 4 200 4 400 4 600	IF UA T ime 0.197 0.466 0.923 0.314 1.294 4.357 0.952 5.178 11.534	[0.05] [0.12] [0.06] [0.06] [0.30] [1.59] [0.38] [2.16] [3.19]	T ime 0.100 0.207 0.401 0.140 0.552 1.420 0.443 1.987 5.148	[0.05] [0.07] [0.14] [0.05] [0.21] [0.48] [0.14] [1.08] [2.15]	TG{UF 49.098 [19.05] 55.545 [18.76] 56.581 [15.94] 55.379 [11.93] 57.367 [12.26] 67.419 [14.60] 53.467 [13.17] 61.633 [12.08] 55.368 [9.49]	Point ¡ to ¡ Box av. Av 30.667 [6.82] 102.907 [32.94] 159.007 [48.80] 64.344 [12.70] 172.334 [53.46] 381.415 [103.48] 109.464 [28.79] 344.285 [118.76] 655.786 [188.35]	n λ 9.150 10.350 10.550 13.350 14.400 14.900 15.800 17.700 18.150	TG{UPUF (%) ¡38.781 [51.13] ¡56.948 [45.88] ¡74.271 [63.22] ¡42.263 [33.19] ¡59.483 [42.72] ¡64.939 [48.52] ¡83.852 [33.32] ¡88.108 [44.97] [3.73] ¡143.201 [43.58] [2.61] [2.51] [2.72] [2.59] [3.08] [2.35] [3.57] [2.93]

Table 2 .

 2 8: Computing time for the sequential method with intermediary filtering, for point-to-box dominance, and relative time gain as a function of |Y s |, for p 2. S up E f 100 in each case.

	Average values for 20 trials.
	|S| |Y s | 2 200 2 300 2 400 3 200 3 300 3 400 4 200 4 300 4 400	IF UA T ime 0.158 0.339 0.668 1.366 2.902 5.886 7.049 15.396 29.178 [12.22] [0.03] [0.04] [0.09] [0.51] [0.85] [1.77] [3.11] [4.78]	T ime 0.168 0.310 0.600 1.313 2.592 5.003 6.840 14.240 26.724 [12.64] [0.03] [0.07] [0.13] [0.52] [0.86] [1.92] [3.18] [4.96]	TG{UF ¡6.345 [12.77] 8.604 [12.96] 10.146 [11.92] 3.892 [4.95] 10.675 [9.14] 15.003 [11.37] 2.956 [5.01] 7.504 [5.72] 8.409 [6.78]	Point ¡ to ¡ Box av. Av 3.172 [0.61] 5.310 [1.20] 6.463 [1.46] 8.697 [2.08] 11.960 [3.49] 18.678 [4.96] 14.164 [3.93] 21.873 [9.11] 28.136 [9.50]	S up E f 21.050 20.450 20.350 29.400 29.200 28.650 41.050 40.100 37.850	[5.06] [4.03] [4.03] [6.03] [5.22] [5.07] [5.47] [4.78] [4.70]	TG{UPUF (%) ¡45.423 [15.16] ¡40.558 [19.40] ¡52.652 [19.84] ¡44.881 [9.77] ¡67.640 [16.74] ¡66.059 [19.23] ¡51.281 [10.41] ¡70.698 [11.06] ¡80.751 [18.89]

Table 2

 2

.9: Computing time for the sequential method with intermediary filtering, for point-to-box dominance, and relative time gain as a function of |Y s |, for p 3. n λ 100 in all instances.

  Ef but smaller TG, probably by lack of diversity or coverage by CHAPTER 2. COMPUTING EFFICIENTLY THE NON-DOMINATED SUBSET OF THE MINKOWSKI SET SUM supported points necessary for point-to-box filtering. Higher n λ will put additional burden on the filtering operation. Whether we should expect optimal n λ to increase with |Y ÝÝÑ s ¡ 1 ¢ Y s | is in our opinion not obvious -as the proportion of supported points depends on other factors than the sheer number of points, and we could not confidently confirm it.

  Point y R p is an upper bound point on Y R p if and only if y I Y, y I © y or equivalently if Y yR p © r. A set of upper bound points will be called a weak upper bound set. The union of the unreachable zones defined by the elements of a weak upper bound set defines an unreachable region, which is especially useful in the context of Branch & Bound dominance relations. Definition 3.2. Point y R p is a lower bound point on Y R p if an only if y I Y, y © y I or equivalently if Y yR p © r. A set of lower bound points will be called a weak lower bound set.

Assuming we are maximizing we give the following definitions, for R p © ty R p | y © 0u: 60 CHAPTER 3. DECOUPLING A COUPLED PROBLEM TO OBTAIN BOUND SETS Definition 3.1. Definition 3.3. (Ehrgott and Gandibleux (2001)) Given Y R p , U B is a strong upper bound set to Y (in the maximization case) if and only if 5 dy Y hy I U B, y I © y dy U B, y is an upper bound point to Y (3.1) or equivalently if Y U B R p

Table 3

 3

	.1: Effect of decomposition on the resolution of an uncoupled problem. Average values for
	10 trials.
	69

  3. DECOUPLING A COUPLED PROBLEM TO OBTAIN BOUND SETS

	p	n	|S| |K| DV	T.ND	|ND|
		200	2	1	10	49.99	[13.16] 426.80	[54.63]
		200	3	1	10	59.27	[14.79] 382.60	[39.54]
		200	4	1	10	70.87	[23.50] 419.00	[65.95]
		200	2	2	10	90.80	[41.41] 420.20	[69.15]
		200	2	3	10	89.24	[20.39] 427.04	[59.82]
		200	2	1	20	82.04	[24.73] 450.50	[76.25]
		100	2	1	10 1663.06 [502.95] 3452.58 [577.36]
		100	3	1	10	933.32 [461.91] 2601.62 [886.65]
		100	4	1	10	859.59 [262.40] 2751.31 [556.69]
		100	2	2	10 1647.33 [1042.60] 3179.33 [1281.37]
		100	2	3	10 2727.04 [1192.15] 4372.53 [893.25]
		100	2	1	20 1971.90 [866.14] 3363.50 [997.37]
		60	2	1	10 2304.04 [1740.95] 3836.13 [1421.05]
		60	3	1	10 1908.62 [1735.73] 3858.21 [2565.50]
		60	2	2	10 2334.97 [1276.06] 3929.67 [1650.62]

Table 3 .

 3 2: Computation of the set of non-dominated points for the original reference instance. Average values for trials.

	p	n	|S| |K| DV	T.RS		T.RS T.ND (%)	|RS|	ApQ.(ε)	|NDRS| |ND|
	2	200	2	10	12.62	[2.33]	25.25 [17.71] 393.82 [62.68] 0.02 [0.01] 0.02 [0.06]
	2	200	3	10	6.69	[0.66]	11.29	[4.46] 361.72 [35.38] 0.03 [0.01] 0	[0]
	2	200	4	10	6.64	[1.69]	9.37	[7.19] 394.68 [70.39] 0.03 [0.01] 0	[0]
	2	200	2	10	13.40	[2.21]	14.76	[5.34] 364.44 [73.81] 0.04 [0.01] 0	[0]
	2	200	2	10	15.64	[3.38]	17.53 [16.56] 360.03 [57.99] 0.06 [0.02] 0	[0]
	2	200	2	20	13.91	[4.39]	16.96 [17.75] 394.15 [66.44] 0.05 [0.01] 0	[0]
	3	100	2	10	92.68 [23.08]	5.57	[4.59] 3148.90 [645.92] 0.02 [0.01] 0.26 [0.26]
	3	100	3	10	29.49 [15.74]	3.16	[3.41] 2258.62 [842.79] 0.03 [0.01] 0.24 [0.31]
	3	100	4	10	16.77	[5.52]	1.95	[2.10] 2559.37 [627.78] 0.03 [0.01] 0.17 [0.30]
	3	100	2	10	73.80 [28.27]	4.48	[2.71] 2703.27 [1086.52	0.05 [0.02] 0	[0]
	3	100	2	10	112.53 [51.92]	4.13	[4.36] 3350.13 [965.83] 0.06 [0.02] 0	[0]
	3	100	2	20	68.26 [20.07]	3.46	[2.32] 2492.41 [770.68] 0.05 [0.02] 0	[0]
	4	60	2	10	128.87 [58.60]	5.59	[3.37] 3910.98 [1697.70] 0.03 [0.02] 0.31 [0.36]
	4	60	3	10	30.29 [17.59]	1.59	[1.01] 3222.32 [1938.03] 0.03 [0.02] 0.25 [0.28]
	4	60	2	10	102.67 [61.53]	4.40	[4.82] 2456.53 [648.04] 0.04 [0.02] 0.18 [0.25]

Table 3 .

 3 3: Computation of the lower bound set associated with the restrict split restriction (RS lower bound). Average values for 10 trials.

	CHAPTER 3. DECOUPLING A COUPLED PROBLEM TO OBTAIN BOUND SETS
	p	n	|S| |K| DV	T.CS		T.CS T.ND (%)	|CS|	ApQ. (ε)	|NDCS| |ND|
	2 200	2	1	10	14.13	[1.74]	28.27 [13.22] 439.90 [62.97] 0.01 [0]	0.28 [0.34]
	2 200	3	1	10	7.52	[0.73]	12.69	[4.94] 396.63 [35]	0.01 [0]	0.21 [0.35]
	2 200	4	1	10	6.81	[1.58]	9.61	[6.72] 400.49 [78.36] 0.01 [0]	0.17 [0.22]
	2 200	2	2	10	17.31	[2.63]	19.06	[6.35] 444.11 [57.59] 0.01 [0]	0	[0.01]
	2 200	2	3	10	23.31	[4.99]	26.12 [24.45] 456.12 [60.25] 0.01 [0.01] 0	[0.00]
	2 200	2	1	20	16.85	[4.46]	20.54 [18.03] 453.32 [83.70] 0.01 [0.01] 0.03 [0.07]
	3 100	2	1	10	104.62 [21.90]	6.29	[4.35] 3741.15 [694.08] 0.01 [0.01] 0.35 [0.34]
	3 100	3	1	10	34.49 [14.35]	3.69	[3.11] 2930.96 [924.46] 0.01 [0.01] 0.29 [0.23]
	3 100	4	1	10	20.03	[5.14]	2.33	[1.96] 2888.02 [538.25] 0.01 [0.01] 0.47 [0.37]
	3 100	2	2	10	107.08 [38.20]	6.50	[3.66] 3639.71 [1304.65] 0.02 [0.01] 0.06 [0.07]
	3 100	2	3	10	146.64 [43.17]	5.38	[3.62] 4588.33 [968.82] 0.02 [0.01] 0.09 [0.16]
	3 100	2	1	20	113.28 [46.60]	5.74	[5.38] 3842.65 [1237.79] 0.01 [0.01] 0.26 [0.20]
	4	60	2	1	10	126.15 [45.67]	5.48	[2.62] 3913.12 [1385.05] 0.02 [0.01] 0.64 [0.38]
	4	60	3	1	10	33.85 [19.53]	1.77	[1.13] 4126.19 [2848.88] 0.01 [0.01] 0.75 [0.28]
	4	60	2	2	10	154.95 [66.09]	6.64	[5.18] 3924.00 [1480.57] 0.03 [0.01] 0.28 [0.36]

Table 3 .

 3 4: Computation of the upper bound set associated with the copy split relaxation (CS upper bound). Average values for 10 trials.

  CHAPTER 3. DECOUPLING A COUPLED PROBLEM TO OBTAIN BOUND SETS

	p	n	|S| |K| DV	T.LR		T.LR T.ND (%)	|LR|	ApQ.(ε)
	2 200	2	1	10	11.08	[1.73]	22.16 [13.15] 426.44 [61.71] 0.05 [0.01]
	2 200	3	1	10	6.54	[0.89]	11.03	[6.02] 416.20 [51.26] 0.04 [0.01]
	2 200	4	1	10	6.56	[2.23]	9.26	[9.49] 452.54 [88.52] 0.04 [0.01]
	2 200	2	2	10	11.01	[2.14]	12.13	[5.17] 421.10 [85.76] 0.09 [0.01]
	2 200	2	3	10	10.82	[2.57]	12.12 [12.61] 351.91 [53.51] 0.12 [0.02]
	2 200	2	1	20	9.62	[2.77]	11.73 [11.20] 419.91 [94.16] 0.09 [0.02]
	3 100	2	1	10	1.21	[0.40]	0.07	[0.08] 3354.34 [787.92] 0.06 [0.01]
	3 100	3	1	10	0.97	[0.49]	0.10	[0.10] 2690.81 [1067.95] 0.04 [0.01]
	3 100	4	1	10	0.94	[0.35]	0.11	[0.13] 2749.46 [776.61] 0.05 [0.01]
	3 100	2	2	10	1.84	[0.63]	0.11	[0.06] 3251.36 [1013.56] 0.09 [0.02]
	3 100	2	3	10	2.37	[1.13]	0.09	[0.09] 3355.25 [1065.60] 0.14 [0.01]
	3 100	2	1	20	2.33	[1.08]	0.12	[0.12] 2932.84 [1195.43] 0.11 [0.02]
	4	60	2	1	10	1.88	[1.37]	0.08	[0.08] 4128.46 [1761.28] 0.05 [0.01]
	4	60	3	1	10	1.80	[1.52]	0.09	[0.09] 3720.42 [2406.99] 0.06 [0.02]
	4	60	2	2	10	1.19	[0.57]	0.05	[0.04] 2920.16 [1112.89] 0.09 [0.02]

Table 3 .

 3 5: Computation of the upper bound set associated with the locally relaxed relaxation (LR upper bound). Average values for 10 trials.

  Graph p1q provides a representation of all possible assignments, as well as the grouping of machines m 1 and m 2 on site s 1 , and machines m 3 and m 4 on site s 2 . Task t 1 can be performed either on s 1 or on s 2 , task t 2 is local to site s 1 , and task t 3 is local to site s 2 . Graph p2q shows a feasible solution, and Graph p3q and p4q show two unfeasible solutions. Solution p3q is unfeasible because task t 1 is assigned twice. Solution p4q in unfeasible because the sum of weights of assignments to machines of site s 2 is 1.1, which is larger than 1.

						m 3
						m 4
	Figure 4.1: The REF problem can be written as the following 0-1 linear program.
	max	ţT mMpSptqq	g	j tm x tm	d j t1, ..., pu
	s.t.	mMpSptqq	x tm ¤ 1	d t T
		ţTpsq mMpSptqq	w tm x tm ¤ b s	d s S
		x tm t0, 1u			d pt, mq N

  REF is a coupled system MOCO, for which the set of decisions is N tpt, mq | t T, m MpS ptqqu which corresponds to the set of indices of decision variables. The set p1q of inequalities describe the coupling constraints, indexed by T C , which corresponds to the set of couplings K in problems formulated in Section C.1.2.1. S naturally corresponds to the set of subsystems.

						CHAPTER 4. APPLICATION PROBLEM
	max	şS ţTpsq mMpSptqq	g	j tm x tm	d j t1, ..., pu
	s.t.	mMpSptqq	x tm ¤ 1		d t T C p1q
		mMpSptqq	x tm ¤ 1		d s S , d t T L psq p2q
		ţTpsq mMpSptqq	w tm x tm ¤ b s	d s S p3q
		x tm t0, 1u			d pt, mq N
						These observations
	allow the following reformulation:

  Which defines an instance of the multiple knapsack problem (MKP,[START_REF] Kellerer | Knapsack Problems[END_REF]).For any s S , the problem of assigning tasks t T psq to machines m Mpsq is a Multi-Choice Thus when T C r, since for all t T , and tsu S ptq, we have MpS ptqq Mpsq, REF can be decomposed into |S | instances of MCKP. This variant is of interest because in practice, when we produce an uncoupled variant of an instance of REF, we obtain a collection of such MCKP instances.

		g t x ts j x ts ¤ 1 Sptq Sptq ş ţTpsq w t x ts ¤ b s	dj t1, ..., pu d t T d s S
		x ts t0, 1u	d t T, s S
	4.1.2.3 Multi-Choice Knapsack variant
	Knapsack Problem (MCKP), formulated as follows :
	max	ţTpsq mMpsq	g tm x tm j	dj t1, ..., pu
		mMpSptqq	x tm ¤ 1	d t T psq
		ţTpsq mMpSptqq	w tm x tm ¤ b s	d s S
		x tm t0, 1u		d pt, mq N

  Thus, for s S , x s K is not constrained in subproblem associated with site s. In REF, K tpt, mq | t T C , m MpS ptqqu, and K tpt, mq | t T L , m MpS ptqqu, Applied to problem REF, this relaxation scheme yields the following problem:

	max	ţT mMpSptqq	g	j tm x tm	d j t1, ..., pu
	s.t.

  DPpte N u, r, T L psq, Mpsq, tsuq 6 Local PoolpS ubResq

		Compute non-dominated complex tasks
		assignments	*/
	2 Local Ð r
	3 S ubRes Ð pr 1 , ..., r |S| q
	4 for s S do
	5	S ubRes s Ð

4

  S ubRes s Ð DPpte N u, r, T L psq T ¦ s , Mpsq, tsuq Weak lower bound sets obtained from greedily completing partial solutions As a weak LB, the image of weak LocRes, i.e. f p ẼpX L , f qq (or alternatively Np f pX L qq) would

	CHAPTER 4. APPLICATION PROBLEM
	4.2.2.3
	5 return PoolpS ubResq
	94

Table 4 .

 4 2: Computing time (second) and approximation quality of the weak lower bound using greedy completions of elements of the strong LocRes lower bound. Average values for 10 trials.

Table 4 .

 4 3: Computing time (second) and approximation quality of the weak lower bound obtained from solving the restrict-split restriction using DP. Average values for 10 trials.

  , wqq 4 /*Assume T C tt C , ..., tu

		*/
	5 for t tt C , ..., tu do
	6	for x Q t¡1 do
	7	for m MpS ptqq do
	8	x Ð xrx tm Ð 1s
	9	if wp xq 0 then
	10

Table 4 .

 4 8: Number of partial solutions eliminated at each step of the basic DP process, using greedy UB and RS LB incumbent set. Average values for 5 trials.enough assignments for the upper bound dominance test to succeed.

  Observation 4.2. Let pX i | Iq be a family of sets with X i X for each i I, and f : X Ñ R p . If iI X i X, then For some instance of REF with feasible set X and σ S T C defining a restrict split variant of REF with feasible set X σ , σS K X σ X Proof. : This follows straightforwardly from the fact that 0 is an admissible value for x tm . : Let x X. For each t, either for all m MpS ptqq, x tm 0, or not. If not, then consider t and m Mp tq such that xt m 1. Then, x X implies that for all m I m, x tm I 0. Thus for each t T C , there is a unique s S such that for all m Mpsq, x tm 0. This means that for each x, there is some σ S T C such that x X σ .Thus we can solve an instance of REF with the simple procedure described by Algorithm 31:

	max	ţT mMpSptqq	g	j tm x tm	d j t1, ..., pu
	s.t.	mMpSptqq	x tm ¤ 1	d t T
		ţTpsq mMpSptqq	w tm x tm ¤ b s	d s S
		x tm t0, 1u			d t T, m MpS ptqq
		x tm 0				d t T C , m Mpσptqq
	Proposition 4.7. Algorithm 25: Resolution algorithm enumerating all decouplings of REF
	1 Out Ð r
	2 for σ S T C do
			Ep	¤	EpX, f qq EpX, f q
					iI
	or equivalently for pY i | i Iq with Y i Y R p for each i I, if	iI Y i Y, then
			Np	¤

iI

NpY i qq NpYq

We consider the case of REF, with σ S T C , the restrict split restriction of which is formulated as follow. Further on, X σ will denote the feasible set of a restrict split variant of REF for σ.

3

Out Ð Out ẼpX σ , f q 4 return ẼpOutq

  |Mpsq| |T L psq| |T C | |S| |N| |S T C | |ND| T.e ¡ const

									NoDec	Enum&Dec
	2	5	5	2	2	70	4	25.725 [10.08]	0.826 [0.405] 0.713 [0.233]
	2	6	5	2	2	84	4	32.25 [10.539] 1.031 [0.441] 0.866 [0.302]
	2	6	6	2	2	96	4	37.775 [13.339] 1.688 [0.696] 1.425 [0.484]
	2	7	6	2	2	112	4	41.05 [15.717] 1.609 [0.707] 1.535 [0.625]
	2	7	7	2	2	126	4	53.775 [15.94]	2.415 [1.138] 2.248 [0.816]
	2	5	5	3	2	80	8	27.925 [10.088] 0.946 [0.372] 1.639 [0.536]
	2	5	5	2	3	105	9	44.525 [13.478] 2.123 [0.934] 1.904 [0.488]
	3	5	5	2	2	70	4	203.325 [106.873] 23.567 [12.061] 7.699 [3.419]
	4	4	4	2	2	48	4	234.9 [122.423] 45.449 [30.068] 14.632 [4.992]

p

Table 4

 4

	p	.9: Instances with few "big" subsystems (average time for 10 trials) |Mpsq| |T L psq| |T C | |S| |N| |S T C | T.e ¡ const |ND| NoDec Enum&Dec
	2	2	2	2	4	32	16	9.075 [5.47]	0.096 [0.084] 0.47	[0.045]
	2	2	2	1	6	36	6	8.425 [5.254] 0.06	[0.047] 0.216 [0.021]
	2	3	3	2	4	60	16	21.675 [10.047] 0.565 [0.253] 0.77	[0.11]
	2	3	3	1	5	60	5	18.625 [8.72]	0.398 [0.249] 0.236 [0.031]
	3	3	3	2	4	60	16	131.475 [72.204] 10.293 [6.849] 1.929 [0.647]
	4	3	3	2	4	60	16	483.6 [264.622] 117.289 [83.511] 37.164 [5.005]

Table 4 .

 4 10: Greater number of 'small' subsystems (average time for 10 trials)

Table 4 .

 4 Out s Ð DPpQ s L , r, T σ psq, Mpsqq ẼpOut σ ¢ Out s q 11 reports |S T C | the total number of decouplings, |ND| the size of the non-dominated set, T.Enum.DP the computing time of the enumerative method where the subproblems of each restrict-split decoupling are solved using the DP, and T.Enum.DP.PreC the enumerative method where the local parts are precomputed and passed to subproblems, at which point the DP procedure is merely completed for remaining tasks.

	Algorithm 26: Resolution algorithm enumerating all decouplings of REF with pre-
	computation of local steps)
	1 Q s L Ð ẼpX s L , pf s , ws qq for each s S 2 Out Ð r
	3 for σ S T C do
	4	Out σ Ð r
	5	for s S do
	6
	9	Out Ð Out Out σ
	10 return EpOutq
	4.4.2.1 Experimental Results

7

/*DP is initialized with precomputed state

Q s 0 */ 8 Out σ Ð

2

  X ν I rx t Ðτpx t ,δqs | δ ∆ t @ where X ν I rx t Ðas 2 x X ν | x t a @ Example 4.4. Consider a coupled problem with feasible set X , S a set of subsystems, T C a set of complex tasks. The partition T C is defined so that each of its elements is the set of variables associated with a coupling. Thus T C ttpt, mq | m MpS ptqqu | t T C u, i.e. each depth of the tree is associated with a complex task, and more precisely with the subset of variables coupled by the coupling. Let us define ∆ t S ptq the set of branching decisions associated with each t T C . Taking decision s S ptq will be represented by τpx t , sq xrx tm Ð 0 | m MpS ptqqzMpsqs, i.e. all assignments of task t to machines which are not in the chosen s are set to 0. Then given node ν at depth t ¡ 1, children nodes ν are defined by tX ν s | s δ t u, with X ν s : tx X ν | x t τpx t , squ. Figure C.13 contrasts this sort of decoupling branching with the binary branching that is classical in branch and bound.

  x tm ¤ 1 d t T C p1q mMpSptqq x tm ¤ 1 d s S , d t T L psq p2q ţTpsq mMpSptqq w tm x tm ¤ b s d s S p3q Figure C.11: Représentation graphique de copy split appliquée à REF C.4.2 Ensembles bornants pour REF C.4.2.1 Relaxation Copy-split Dans le cas de REF, cette relaxation revient à remplacer chaque affectation de tâche complexe par |S | copies de cette affectation, une chaque sous-système, comme illustré par la Figure C.11. On obtient alors un problèmes décomposable en |S | problèmes de la forme suivante, pouvant être résolu par programmation dynamique.

	t 1		m 1	t 1	m 1
			s 1		s 3
			m 2	t 2	m 2
	t 2
			m 3	t I 1	m 3
			s 2		s 4
	t 3		m 4	t 3	m 4
	max	ţTpsq mMpSptqq	g tm x tm j	d j t1, ..., pu
	s.t.	mMpsq	x tm ¤ 1		d t T psq
		ţTpsq mMpSptqq	w tm x tm ¤ b s
		x tm t0, 1u	d t T psq, m MpS ptqq
		x tm t0, 1u		d pt, mq N

C.4.2.2 Singleton borne supérieure forte.

This chapter is based on submitted article[START_REF] Kerbérénès | Computing efficiently the non-dominated subset of the minkowski set sum[END_REF]

An immediate development of our work on box-based dominance relations should provide a boxing algorithm with parameters that are independent of the value of p. One line of research would consider a trade-off between two criteria. The first one would be a desired minimum number of points per box, so as to control the number of points eliminated per box-based dominance test. The second one would be the maximum acceptable distance, over criteria, between a point and the box it belongs to, so as to keep the upper bound of a box tight enough.

This chapter is based on submitted article[START_REF] Kerbérénès | Bound sets for multiobjective optimization of complex systems[END_REF]

Earlier versions of this chapter were presented orally at the RAMOO

and ROADEF

conferences.

return PoolpLocal, Complexq

Remerciements Acknowledgement

This work was partly funded by the French National Agency of Research and Technology (ANRT), CIFRE grant number 2017/1000, and supported by NAVAL GROUP RESEARCH.

Chapter 4

Application Problem 1

Chapter Abstract

In this chapter, we will study an application problem called REF, a particular case of the generic coupled problem. REF is a multi-site assignment problem with coupling constraints, and admits a dynamic programming resolution algorithm. We apply the various notions of uncoupled lower and upper bounds to REF, and provide dynamic programming algorithms to compute them. We show that decomposition as well as other computational tricks can improve the initial dynamic programming approach. Facing some intrinsic limitations of the initial "bottom-up" approach, we switch the perspective and use dynamic programming as a subprocedure in a "top-down" approach which tries to enumerate, as efficiently as possible, all ways of decoupling the problem. Table 4.5: Using partial filtering over extensions built from tasks assignments in the same site, at each step of the DP process. Average values for 10 trials.

Results reported in table 4.4 show that our DP algorithm benefits most from using filtering using the RS lower bound as incumbent. Although it was in several cases more costly to compute, it led in all cases to a quicker resolution of the problem, which, in the biobjective case, led to solving the problem up to 15.6 times faster than the basic DP algorithm, and about 4 times faster for p 3 on the tested instances. Thus, filtering using an incumbent set obtained from a weak lower bound seems to be a very reliable way to improve the performance of the DP based approach.

Experiments on partial filtering presented in Table 4.5 showed that, for p 2, this trick provides a slight but systematic improvement of average performance. For p 3 however we did not see an advantage to partial filtering. In any case, the average improvement yielded by partial filtering was smaller than the standard deviation and thus probably not reliable in practice. We believe there may be ways of further improving this idea by keeping in memory, from the first partial dominance test, which dominance tests will lead to an incomparability when all partitions of the set of states are pooled together for the main dominance filtering. bound being of course counted in the computing time for KS DP) does allow some performance improvement. We achieved additional performance improvements by combining kickstarting with the use of filtering using the RS LB as set of incumbent solutions, and applying filtering both in the local subproblems and when assigning complex tasks. Kickstarted DP works especially well when the number of complex tasks is small in proportion of the number of local tasks. However, in the opposite case, it appears to perform very similarly to the basic DP algorithm, especially so when the number of criteria increases.

Finally, experiments reported in

Table 4.7 summarizes the results obtained in this subsection, to determine the best "bottomup" method available. The clear winner is the approach filtering with the RS LB incumbent set over the basic DP sequence, rather than the kickstarted one.

This may be surprising at first, but to understand it, one may note that occurrences of filtering of partial solutions using the incumbent set in the decoupled subproblems of the KS approach were virtually inexistent. Indeed, as partial solutions receive assignments, the upper bound over their extensions becomes tighter. In a short experiment reported in Let σ S T C identify one of the possible restrict-split instances, with feasible set denoted X σ , a restriction of the original feasible set X. If it was the case that for all σ S T C , Np f pX σ qq dominates Np f pX τ qq, then σ would obviously be the best decoupling. However, this is extremely unlikely if the profits associated for each criterion with task assignments are uncorrelated.

Conversely, we can determine a decoupling σ S T C to be undesirable a posteriori, either if it is dominated by Np f pXqq, or if there exists τ S T C such that Np f pX τ qq dominates Np f pX σ qq.

The second case of course implies the first case by transitivity. In both case, we may be able to use multiobjective branch and bound reasoning to avoid enumerating all decouplings. However, if only the first case is met, because it may be that only the union of non-dominated sets associated with a subset of leaf-instances can dominate Np f pX σ qq, a heavier computational effort will be required for a successful bound test.

In the following experiment, we will solve a number of instances and measure a posteriori the proportion of leaf instances such that their non-dominated set is dominated by the nondominated set of the whole problem (column Lvs dtd ND), and proportions of leaf instances such that their non-dominated set is dominated by the non-dominated set of another leaf-instance (column Lvs dtd Lvs).

Results are reported in Table 4.14. We found, as is necessary, that Lvs dtd Lvs is always smaller than Lvs dtd ND, but Lvs dtd ND was surprisingly large, although profit coefficients were chosen randomly and uncorrelated, it appears that an overwhelming majority (up to about 96% in our experiments) of decouplings can be dominated by a single other decoupling. This makes possible to consider a branch and bound (BB) approach to REF in which we would branch over the decisions to assign complex tasks to sites, rather than to machines. Thus, we now present branch and bound (BB) notions directly in the multiobjective case, and in a way that will allow us to easily introduce the unusual 'decoupling" branching scheme particularily suited for decomposition.

Application to the resolution of biobjective REF

The enumeration of all decouplings described in Section C.4.4.1 will be performed as an arborescent search, in which leaf-nodes are associated with feasible sets X σ for σ S T C . First, we compute, as an incumbent set, the non-dominated set of one restrict-split variant, that which maximizes the heuristic described in Section C.3.4. Of course, we will not be solving this variant again during the arborescent search. We also pre-compute, as in Section C.4.4.2, the local restriction which we will use to initialize each leaf-node problem. To use our BB approach, the following assertion is needed: Observation 4.3. Proposition C.5 implies that the decoupling branching scheme is exhaustive.

We explore this arborescence with a depth first strategy. At each intermediary node, we try to eliminate active search zones at this node, separating them with the supporting hyperplanes of solutions to weighted sum scalarizations of the linear programming relaxation of the node problem. The weight vectors orienting these hyperplanes are generated according to the dichotomic method for the generation of supported solutions. Because we need to limit the number of supported solution we compute at each node, we defined an ε-dominance threshold. When two consecutively generated supported solutions ε-dominate one another for a value of ε smaller than this threshold, we stop. We tuned the value of the threshold manually so as to maximize the ratio of the number of zone eliminations over the number of generated supported points.

Note that we do not update the search regions when new feasible solutions are found in leaf nodes, because this process proved too costly in preliminary experiments. The overall procedure is described by Algorithm 28, where σ ν S T C denotes the assignment of complex tasks to sites which is associated with leaf node ν, and DPpreComp denotes the sub-procedure associated with lines 6 to 9 of Algorithm 26, where a restrict-split variant is solved using DP initialized with the precomputed decision subsequences associated with local tasks.

In the following experiment, we compare the time needed to solve REF using the e-constraint method, enumerating all restrict-split variants with precomputation of the local subproblems as in C.4.4.2, with the running time of the decoupling branching approach. This experiment is limited to the biobjective case, so as to take advantage of the dichotomic scheme for generating pieces of the separating hypersurface. Instances are generated in the same way as for previous experiments with REF. In Table 4.15, T. Enum. BB denotes the decoupling branch and bound approach, with #S aved denoting the sum of the number of leaf nodes in the subtrees which were avoided thanks to bounding.

We can observe that for all tested series of parameters except the last one, we found the approach with bounding to perform slightly better than the enumeration of restrict-split variants with partially precomputed dynamic programming. In cases where the number of restrict-split variants to be solved was too high to allow the enumerative approach to beat the e-constraint method, using BB tricks could not reverse the tide, and e-constraint remains superior. Because we did not attempt to enrich the search region with feasible solutions obtained at leaf-nodes, we were limited Algorithm 28: Decoupling branch and bound input : L an incumbent set, Q L pQ s L | s S q the precomputed last steps of each local tasks subproblem

Appendices

Appendix A

The double bound as a solution concept in MOCO

A.1 Definitions

In some operational situations, knowledge of all non-dominated points may not be the most useful information to a decision maker, especially in contexts that require quick decision making. In such contexts, the decision maker may be interested in being presented with bounds to the nondominated set, so that she can use this information to trigger another multicriteria decision method.

For example, the information provided by either a strong or weak upper bound, and by either a weak or strong lower bound, may help the decision maker select an aspiration point to be approached by means of a scalarized method, knowing which level of performance are atteinable. In other contexts, when both a feasible, lower bound solution and an upper bound point are provided, the latter may be used as a quality estimate of the former. The decision maker may decide to settle for the feasible point if the closest upper bound point that dominates is sufficiently close to it. Such a double bound materializes the uncertainty of the decision maker regarding the location of the non-dominated set, in the form of an area or hypervolume. One way to quantify this uncertainty, given X the feasible set of the problem and f the vector valued objective function, is to produce an estimate of the probability that a point drawn randomly and bounded by an upper bound and a lower bound on f pXq will fall "between" the LB set and the U B set we computed.

What we mean by "between" needs to be made precise and depends on whether the U B and LB sets are weak or strong bounds to Np f pXqq. The correctness of our measurement requires that the lower bound of the hyperrectangle in which we sample points be a lower bound point to f pXq, and accordingly that the upper bound of the same hyperrectangle be an upper bound on f pXq. Assuming that for each j t1, ..., pu, f j is linear with positive coefficients, we can take 0 as a lower bound. As for the upper bound on the hyperrectangle, we propose two alternatives. They differ with regards to the sort of a priori information about f pXq that they require from the decision maker. In the first case, under the same hypotheses about f pXq, the only required information are that coefficients p i j be known. Then we set u max such that for each j t1, ..., pu, u max j °n i1 p i j . It is easily seen that this is the maximum of f j pxq for x t0, 1u n , i.e. if the optimization of f j was unconstrained, and thus y ¨umax for any y f pXq. In the second case, we use the ideal point of f pXq as upper bound.

The use of the ideal point as an upper bound on random points is well suited to the case where our measurement is meant to quantify the decision maker's uncertainty, but it is not well suited when it is meant to quantify the quality of bound concepts that do not involve the ideal point itself. Indeed, in some cases, the ideal point may be dominated by some point in the upper bound set, and thus constitute by itself a tighter upper bound, the knowledge of which reduces uncertainty. In much the same way, the lexicographic optima obtained by computing the ideal point constitute lower bound points which may dominate points in the (weak) lower bound under scrutiny. In a context when uncertainty is to be quantified, one should add these points to the bound sets, and add the time spent computing the ideal point to the time spent computing the double bound.

A.2 Experimental Results

In practice, we draw 10000 points uniformly randomly in a hypersquare box of lower bound 0 and upper bound u max or ideal point. In this experiment, LB is always a weak lower bound set, so for each random point r drawn, we call the draw a success if there is no y LB such that r ¨y, and if there is a y I U B such that r ¨yI . Our metric is then simply the proportion of successful draws For any value of p, and any parameter set, we find that a double bound obtained using the RS and lexicographic optima weak lower bound, and the CS and ideal point upper bound yielded the smallest area of uncertainty of all four combinations of bound notions. Expectedly, the size of the area of uncertainty for any notion of double bound clearly increases with an increase in |K| and in DV, while it appears to slightly increase with an increase of n, although for RS {CS , this trend it tainted with great variability.

Appendix B

Weak lower bounds obtained from correcting weak upper bound solutions

B.1 Correction of solutions to the copy-split relaxation

Recall that the copy-split relaxation of REF, which is equivalent to a variant in which for each t T C , the collection of assignments variables px tm | m MpS ptqq is replaced by ppx tm | m MpS ptqqq | s S q i.e. by a collection of copies of the complex task assignments, each restricted to a subsystem. This relaxation is formulated as max ţT mMpSptqq This relaxation, in practice, gives an upperbound which is a lot looser than the relaxation of coupled variables from local knapsack constraints, and contrary to the latter, its quality decreases as the number of subsystems increases. However we can use it to another end. Solutions to this problem may be infeasible as they violate the assignment constraints of complex tasks assignment variables, but they can be made feasible easily by just choosing one among the several assignments of a complex task. Given a solution for the relaxation, we make the choice of which task assignment to undo in the following way. Contrarily to greedy approaches where we want to add assignments, here we sort assignment, for each j t1, ..., pu, by increasing order of g j tm w tm , and then by increasing sum of their objective-wise ranks. Thus, we undo the "worst" assignments first, freeing more budget by sacrificing less profit. Until all complex tasks are assigned at most once, we consider in order each assignment pt, mq. If it has been made in the solution, while the task has also be assigned to another machine m I , we undo assignment pt, mq. This is described by UPPER BOUND SOLUTIONS Algorithm 29.

Algorithm 29: Greedy correction of solutions of the copy-split relaxation input : Asg, EpX I , f q output: Corrected X 1 Corrected Ð r 2 for x EpX I , f q do 3 while x X do

Solutions to this relaxation may be such that the original knapsack constraint, taking into account the coupled variables, are violated. To obtain a feasible solution, we therefore need to undo some of the assignments done in these solutions. Since we want to obtain a feasible solution that is as close as possible to the Pareto set, we wish to undo assignments so as to decrease objective value as little as possible. Here consider assignments pt, mq for t T, m MpS ptqq in the same order as for the previous correction concept, yielding sequence Asg. Until the knapsack constraint for each of the subsystems is satisfied, we consider in order each assignment, and if it is done in a subsystem for which the knapsack constraint is violated, we remove it. This procedure is described by Algorithm 30. UPPER BOUND SOLUTIONS Algorithm 30: Greedy correction of solutions of the LR relaxation input : Asg, EpX I , f q the set of efficient solutions of the LR Relaxation. 1 Corrected Ð r 2 for x EpX I , f q do 3 while x X do

B.3 Experimental results

We perform experiment regarding the computation of an incumbent set from corrections of solutions to the CS and LR upper bounds, a process which we have described in Sections B.1 and B.2 respectively. We compare the time required to obtain these incumbent sets from the upper bound computed either with the e-constraint method or the DP method, and we measure how well they approximate the efficient set of the original problem. From Table B.1 we can observe that, because the correction process is quick and linear, these incumbents are obtained in virtually the same time as the bound sets they originate from. However, this time is itself rather long, and generally less than 10 times quicker than the resolution of the whole problem, except of course when p ¥ 3. Finally, we can observe that correcting solutions from the LRU B usually provides a tighter bound, which makes it the better of the two options in cases the LRU B is computed quicker than CS U B (i.e. for low values of |S |). Definition C.3. Une solution x X est efficace par rapport à f si et seulement si il n'existe acun x I X, x I tel que f px I q © f pxq, et est faiblement efficace si et seulement si il n'existe aucun x I X tel que f px I q ¡ f pxq.

Ainsi, il est fréquent de considérer comme concept de solution l'énumération de l'ensemble des solutions effiaces ou l'ensemble de ponts non dominés:

EpX, f q : tx X | x est efficace relati f uNpYq : ty Y | y is non-dominatedu Definition C.4. Une solution efficace x I EpX, f q est supportée s'il existe λ pλ j | j t1, ..., puq R p ¡ tel que x I arg max xX °p j1 λ j f j pxq. Quand p 2, la méthode dichotomique de [START_REF] Aneja | Bicriteria transportation problem[END_REF] énumère tous les points supportés.

Definition C.5. Pour y py 1 , y 2 , ...y p q et y I py I 1 , y I 2 ..., y I p q R p , y © Lex y I si et seulement si y 1 ¡ y I 1 , ou y 1 y I 1 et py 2 , ...y p q © Lex py I 2 , ..., y I p q. Alternativement, le filtrage peut être fait par insertion dans un KDTree au sens de Chen et al. (2012), au sein d'une méthode en deux phases décrite par l'Algorithme 2. Un point dénommé sieve permet de commencer la série de tests par un test ayant de bonnes chance de réussir. Une routine de pruning permet un filtrage intermédiaire du KDTree et une amélioration des performances.

C.1.1.2 Choix d'un algorithme de filtrage par dominance

Quand p 2, on utilise l'algorithme unidirectionnel décrit par [START_REF] Kung | On finding the maxima of a set of vectors[END_REF], avec filtrage en temps constant. Quand p ¥ 4, on utilise la méthode en deux phases avec insertion dans un KDTree de [START_REF] Chen | Maxima-finding algorithms for multidimensional samples: A two-phase approach[END_REF]. Qand p 3, il faudra décider au cas par cas entre unidirectionnel et insertion dans un KDTree.

C.1.2 Problèmes coupés et décomposition C.1.2.1 Formulations décomposables

Soit N t1, ..., nu l'ensemble des variables décisions partionné par S . Un sous-système s S est identifié par un sous-sytèmes de variables. Thus, x px 1 , ..., x n q, ou px s |s S q, où x s px i | i sq px s 1 , ..., x s n s q. Les interactions entre sous-systèmes sous représentées par des contraintes associant des variables de différents sous-systèmesL Soit K PpNq un ensemble de couplages représentés par les ensembles d'indices de variables couplées.

Definition C.6. s, s I S , s s I sont couplés par la contrainte k s'il existe i s k et j s I k. Un problème d'optimisation est non couplés si et seulement si K r.

C.3.2 Restrictions décomposables

Pour obtenir une restriction non-couplée du problème original, on fixe toutes, ou certaines variables couplées, à des valeurs admissibles ou neutres.

Definition C.14. e i est une valeur admissible pour x i si il existe x X tel que x i e i , i.e. t.q. px ¡i , e i q X. Definition C.15. e i est une valeur neutre pour x i si pour tout x X, px ¡i , e i q X.

x s K px s k | k Kq dénote les variables du système s couplées par la contrainte couplante k. x s px s K , x s K q, avec x s K le vecteur des variables non couplées de s. Soit x K le vecteur des variables non couplées. On obtient la restriction LocRes et la borne inférieure associée en fixant toutes les On fixe les variables apparaissant dans la contrainte k, sauf celles de s.

Example C.2. L'exemple suivant illustre la variante restrict split découplée du problème original, obtenue en fixant les variables des sous-systèmes autres que s 1 qui apparaissent aussi dans la contrainte couplate k.

C.3.3 Relaxations décomposables

Une relaxation est obtenue en remplaçant les fonctions objectif f j par des fonctions relâchées f I j ou en considérant un sur-ensemble X I de X, obtenu en relâchant les contraintes de X. Le plus simple est d'ignorer ces contraintes. Pour obtenir la décomposabilité, l'on peut se limiter à ignorer les termes de fonctions de contraintes additives associées avec certains sous-systèmes. Ces contraintes sont "relâchabels" en ce sens. Par exemple les contraintes de budget linéaires sont relachable en ignorant un sous ensemble quelconque de variables.

C.3.3.1 Copy-splitting des contraintes couplantes

Si la contrainte couplante k est relâchable pour chaque s γpkq, soit g Is k px s k q ¤ b k la relaxation de g k px k q ¤ b k associée à s |γpkq|. Alors le problème suivant est une relaxation de pPq, décompos-

Example C.3. Le schéma suivant illustre le remplacement de la contrainte k par deux contraintes k I and k P , copies de k restreinte à chaque sous-système.

C.3.3.2 Relaxations locales des variables couplées

On suppose que chaque x s X s est du type g s px s q ¤ b s . Pour K l'ensemble de variables non couplées g s px s q ¤ b s est relâchée en g Is px s K q ¤ b s . On obtient: max f pxq

Le problème est alors décomposé d'une part en la partie ne comprenant que les variables couplées, et d'autre part en chaque sous-système du problème original.

C.3.4 Résultats expérimentaux

On résout un problème générique couplé, puis on calcule à la fois des ensembles bornant supérieurs et inférieurs par décomposition. Le but est de comparer le temps de calcul nécessaire à l'obtention de ces bornes avec le temps de calcul nécessaire à la résolution du problème et d'estimer la qualité d'approximation de la solution du problème original permise par les ensemble bornants issues des relaxations et restrictions. La Figure C.9 montre l'encadrement de l'ensemble de solutions non-dominées original par des ensembles bornant inférieurs et supérieurs faibles.

Résoudre une variante restrict-split demande de produire une affectation des contraintes couplantes aux sous-systèmes. Ce choix est réalisé de manière heuristique en calculant:

Concernant les contraintes relâchées de la variante copy-split, pour tout k K, si °sS °isk a i x i ¤ b k , alors, comme les a i sont positifs, on a pour tout s S , °isk a i x i ¤ b k . D'où l'expression: max f pxq

Les contraintes de la variante locallement relachée on considérant que du fait de la positivité des c i on a °is s K c i x i ¤ °is c i x i , où s s K signifie z kK k. Ainsi °is s K c i x i ¤ d s est une relaxation de °is c i x i ¤ d s , qu'on lui substitue dans la variante: max f pxq

Pour mesurer la qualité d'approximation permise par un ensemble bornant, on utilise la notion de ε-dominance a posteriori: In the following tables, RS denotes the restrict split restriction (Section 3.2.2) of the original problem, which provides a lower bound set. CS denotes the copy split relaxation (Section C.3.3.1) of the problem, which provides an upper bound set. LR denotes the relaxation obtained by omitting coupled variables in local subproblem constraints (Section C.3.3.2), which also provides an upper bound set. T. denotes computing time, |.| denotes the size of a set of points, and ApQ.pεq denotes the quality of approximation in terms of a posteriori ε dominance.

Dans le cas bi-objectif, la solution est calculée rapidement, donc le coût du calcul des ensem-bles bornant peut être considéré comme trop lourd. Au delà cependant, il devient négligeable.

Empiriquement, on observe que si le nombre de sous-systèmes augmente à nombre total de variables constant, le temps de calcul de toute variante décomposable diminue, puisque ces variantes se divisent en plus petit problèmes. Ainsi la borne inférieure associée à restrict-split est obtenue en moins de 2% du temps nécessaire à la résolution du problèmes quand |S | 4. Example C.5. Admettons T tt 1 , t 2 , t 3 u, M tm 1 , m 2 , m 3 , m 4 u et S ts 1 , s 2 u, tels que S pt 1 q ts 1 , s 2 u, S pt 2 q ts 1 u, S pt 3 q ts 2 u, Mps 1 q tm 1 , m 2 u, Mps 2 q tm 3 , m 4 u. Supposons que le budget de ressource de chaque site est de 1.

C.4.1.2 Algorithme de programmation dynamique

Soit 0 N la solution telle que pour tout i N, 0 i 0, faisable par définition. T t1, ..., tu, l'ensemble des tâches dans REF, correspond à l'nesemble des étapes de décision. Soit Q 0 : t0 N u l'état initial du processus de décision. Les états des étapes suivantes associées à chaque t T sont obtenus en construisant les solutions réalisables par affectation de chaque t aux machines m MpS ptqq. τpx, δq désigne la solution obtenue en prenant la séquence de décisions δ.

Comppxq : f j pxq ¥ f j px I q ñ f j pτpx, δqq ¥ f j pτpx I , δqq dj t1, ..., pu and ws pxq ¥ ws px I q ñ ws pτpx, δqq ¥ ws pτpx I , δqq ds S La correction de l'algorithme de programmation dynamique est garantie par Proposition C.3. L'algorithme de programmation dynamique produit l'ensemble Q t tel que

Une instance d'un problème de programmatino dynamique est décrite par pQ 0 , T C , T L , M, S q, et on écrit B Ð DPpQ 0 , T C , T L , M, S q quand on obtient la préimage ẼpX, f q de Np f pXqq par programmation dynamique.

C.4.1.3 Experimental results

On teste l'algorithme basique de programmation dynamique sur des instances où les taches complexes ont accès à tous les sites. Les budget de ressources des sites sont définis de manière à correspondre à des instances de sac à dos difficiles.

On trouve que l'algorithme de programmation dynamique est préférable à la méthode générique e-contrainte quand p augmente, mais pour p 2, quand |T C | ou |S | augmente, il devient beaucoup plus lent, du fait de la p |S | dimensionalité des états du processus de DP, qui induit de nombreuses incomparabilités et l'accumulation du nombre de solutions partielles. Algorithm 31: Algorithme de résolution énumérant les découplages de REF

Empiriquement, on observe que même l'énumération de S T C peut être compétitive et plus rapide que l'exécution d'une méthode générique, quand le nombre de sous-système est relativement faible relativement à la taille des sous-systèmes.

C.4.4.2 Pre-computation of independent dynamic programming states

La variante restrict-split est résolue par décomposition, et chacun de ses sous-problèmes est résolu par programmation dynamique avec p 1 critèrs seulement, et le filtrage des combinaison se limite à p critères. Ces variantes comprennent certaines tâches complexes devenues simples, mais les tâches déjà simple dans le problème original devraient être traitées de la même manière dans toutes les variantes. La sous-séquence de programmation dynamique associée avec ces tâches peut donc être précalculée et utilisée pour initialiser chaque variante découplée considérée. On considère l'ensemble réalisable X s L du sous-problème s réstreint à ses tâches locales, défini comme :

On calcule d'abord Q s L ẼpX s L , pf s , ws qq pour tout s S . Puis tout tous σ S T C , on complète l'ensemble des |S | sous-problèmes indépendants avec les tâches rendues simples pour chaque sous-problème par σ. On combine et filtre les solutions des sous-problèmes pour obtenir la solution de la variant restrict split. On calcule ensuite l'union des ensembles non-dominés pour tous les σ S T C , que l'on filtre.

Empiriquement, on observe que l'utilisatino de la programmation dynamique et du précalcul amène à une amélioration significative de la performance: pour des instances faiblement couplées, même dans le cas bi-objectif, le problème est résolu 4 à 8 fois plus vite qu'avec le l'algorithme générique. Pour p 3, de 18 à 52 fois plus vite. Cette approche permet également une résolution beaucoup plus rapide qu'avec l'algorithme de programmation dynamique classique, même amélioré par le filtrage par une borne inférieure. Toutefois, si le degré de couplage de problème augmente (nombre de contraintes couplantes ou de sous-problèmes, donc valeur de S T C ), la méthode générique redevient préférable dans le cas biobjectif.

C.4.4.3 Recherche arborescente et branching

Pour améliorer la méthode d'énumération des découplages, on la représente sous forme d'une arborescence, telle que chaque niveau correspond à une tâche complexe, et chaque décision de branchement à un sous-système auquel associer la tâche complexe. Les feuilles de cette arborescence correspondent à l'ensemble réalisable d'une variante restrict split du problème original, que l'on peut résoudre par décomposition et programmation dynamique. Un algorithme de branch & bound est correct si le schéma de branchement est exhaustif au sens suivant (ce qui est impliqué par la proposition C.5).

Definition C.17. E is exhaustive if and only if for all ν V ¤ ν I s.t.pν,ν I qE

L'énumération est faite en profondeur d'abord. On calcul d'abord un ensemble de solutions candidates, qui correspondent à la solution de la variante restrict-split choisie par l'heuristique de la Section C.3.4, qui sera exclue de la recherche arborescente. On tire aussi partie du précalcul des sous-séquences de programmation dynamique décrit précédemment. La question est de savoir si des élimination de sous-arbre de type branch & bound peuvent être effectuées dans de telles arborescences.

Dans le branch & bound monocritère, il s'agit de montrer qu'une solution réalisable déjà connue est meilleure qu'une borne supérieure sur les solutions réalisables d'un sous arbre. La généralisation de cette notion au cas multiobjectif nécessite la définition de la "région de recherche", qui Pour toutes les instances testées à part une, l'approche branch & bound est légèrement plus rapide que l'approche énumérant les découplages avec décomposition, programmation dynamique et précalcul. Dans tous les cas, on voit qu'une proportion significative des noeuds feuilles de l'arborescence ont été évités. Des travaux supplémentaires sont nécessaires pour rendre cette exploration plus efficace et généraliser la séparation au delà de deux critères.

MOTS CLÉS

Optimisation Multiobjectifs, optimisation combinatoire, systèmes complexes, programmation dynamique, décomposition RÉSUMÉ Un système complexe peut être vu comme une collection de sous-systèmes irréductiblement liés mais assez indépendants pour être distingués. Cette thèse considère l'optimisation de systèmes complexes dans le cadre multi-objectifs. L'interaction entre sous-systèmes d'un système complexe y prend la forme de contraintes couplantes, c'est-à-dire faisant intervenir des variables issues de différents sous-systèmes. Après avoir rappelé les fondements de l'optimisation multi-objectifs et de l'algorithmique de filtrage par dominance, nous présentons la notion de système couplé, et définissons dans le cas multi-objectifs celle, centrale, de décomposition. Une implémentation simple de la décomposition suffit à améliorer le temps de résolution de problèmes non-couplés. Nous proposons de surcroit des méthodes algorithmiques avancées pour la combinaison de solutions de sous-problèmes et l'élimination des combinaison dominées, utilisant les notions de boîte bornante et d'algorithme unidirectionnel de filtrage par dominance. Le défi principal de l'optimisation de systèmes complexes reste de prendre en compte les contraintes couplantes, tout en évitant de considérer l'entièreté du problème original en même temps. Nous proposons des restrictions et relaxations génériques des contraintes couplantes de problèmes couplés, permettant d'obtenir des ensembles bornant supérieurement et inférieurement l'ensemble de solutions non-dominées. Nous montrons que ceux-ci peuvent être calculés en tirant parti de la décomposition. Enfin, nous présentons un problème d'application : une affectation multi-site multi-objectifs sous contraintes de ressources. Nous montrons que ce problème admet un algorithme de résolution par la programmation dynamique, et comment la décomposition peut être utilisée pour améliorer cette méthode initiale. D'une part, le processus séquentiel de décision peut lui-même être décomposé en sous-séquences indépendantes. D'autre part, des bornes ou des ensembles bornants obtenus par décomposition peuvent être utilisés pour accélérer le processus séquentiel de décision par l'élimination précoce de solutions partielles. ABSTRACT A complex system is a collection of subsystems which are independent enough to be distinguished but linked together in significant ways. This thesis considers the optimization of complex systems within the framework of multiobjective optimization and focuses on a representation of complex systems as coupled systems, meaning that the interaction between subsystems is modeled by coupling constraints, i.e. constraints involving variables from different subsystems. After having recalled the basic notions of multiobjective optimization and of dominance filtering algorithms, we introduce coupled problems, and define the key concept of decomposition in the multiobjective case. A simple implementation of decomposition already yields performance improvement in the resolution of uncoupled problems, but we provide further algorithmic improvements to the combination of solutions from subproblems and the elimination of dominated combinations, using notions of bounding boxes and of unidirectional dominance filtering algorithms. Beyond the uncoupled case, the main challenge of complex systems optimization remains to take coupling constraints into account, while never having to consider whole complex system optimization problem at once. We propose generic restrictions and relaxations of coupling constraints, which yield upper and lower bounds set on the set non-dominated set of a coupled problem. We show that bound sets can be obtained using decomposition. Finally, we present an application problem: a multiobjective multilocation assignment problem. We show that it admits a dynamic programming resolution method, and we show how decomposition can be used to improve on this initial resolution method. On the one hand, the sequential decision process can itself be broken down into independent subsequences. On the other hand, reasoning using bounds or bound sets obtained by decomposition can be used to speed up the sequential decision process by eliminating partial solutions early.

KEYWORDS

Multiobjecitve Opitmization, combinatorial optimization, complex systems, dynamic programming, decomposition