Camille Raymond
email: camille.raymond@altran.com

Daniel Prun
email: daniel.prun@enac.fr

Extending MBSE methodology and SysML formalism to integrate human considerations

Keywords: System Engineering, Human-System Integration, SysML, MBSE

In response to the challenges related to integrate human considerations within complex interactive systems, the Human-System Integration (HSI) initiative has been proposed.

Nevertheless, the HSI approach is difficult to fit efficiently into the system engineering processes because of the lack of shared methodologies, modeling semantics and tools. This paper focuses on the needs for a new methodology and an adapted formalism to improve the collaboration between the system engineers and HSI practitioners. It concludes with the presentation of our methodology based on the Model-Based System Engineering (MBSE) and the SysML modelling language.

INTRODUCTION

Integrating human aspects within complex systems raised many issues in system design. Usually, industry tends to refer to work procedures to supervise operators during their activities. These measures showed some limits as they considered human as a limitation factor in a system. Indeed, in seeking to restrict human influence, they consequently limit human capacities of reflection and adaptation.

Thus, the Human-System Integration (HSI) approach has been proposed in favor of improving the considerations of the human characteristics in interactive systems design. Nevertheless, this alternative is difficult to fit efficiently into system engineering processes due to the inadequate methodologies, formalisms and tools allowing an efficient inter-integration of both disciplines [START_REF] Booher | Handbook of human systems integration[END_REF]. This paper focuses on the need for a methodology to enhance the collaboration between system engineers and HSI practitioners through MBSE methodologies and SysML formalism. It presents our contribution lines through an aeronautical maintenance study case.

SYSTEM ENGINEERING AND MBSE

The system engineering is defined as "an engineering discipline whose responsibility is creating and executing an interdisciplinary process to ensure that the customer and stakeholder's needs are satisfied in a high quality, trustworthy, cost efficient and schedule compliant manner throughout a system's entire life cycle" [START_REF]INCOSE. Systems Engineering Handbook: A Guide for System Life Cycle Processes and Activities[END_REF]. The ISO/IEC 15288 [START_REF]Systems and software engineering -System life cycle processes[END_REF] defines the processes involved in a project life cycle including: requirements analysis, architectural design, implementation, integration, verification and validation.

A system is expressed by a set of concerns (technical, organizational, financial, safety, human considerations etc.) and implies various stakeholders during the system life cycle (client, engineers, human factors practitioners, developers, testers, etc.). It is necessary to have standards, methodologies and tools to carry the system engineering activities. They allow to express and evaluate different system architectures tradeoffs, and elect afterthought the most appropriate solution satisfying all the concerns according to the resources available. Thus, the system engineering community proposes the Model-Based System Engineering (MBSE) methodologies to support the system engineering processes. MBSE is defined as the "formalized application of modelling to support system requirements, design, analysis, verification and validation activities beginning in the conceptual design phase and continuing throughout development and later life cycles phases" [4]. MBSE methodologies depict the common system architecture breakdowns: operational layers, functional layers, logical layers and physical layers.

The MBSE methodologies substitute the traditional document-centric approach to a model-centric approach. The document-centric approach tends to produce a large amount of documents with various types of artifacts, potentially used in an inconsistent manner and ambiguous due to the use of natural language. The MBSE methodologies rely on centralized and up to date models, shared with all the stakeholders. They provide greater system design consistency and completeness, help to manage architecture complexity and improve the communication between stakeholders [START_REF] Estefan | Survey of model-based systems engineering (MBSE) methodologies[END_REF].

Contrary to the document-centric approach, models provide the advantage to be dynamically manipulated and transformed. Through the viewpoint mechanism, MBSE allows the representation of models from the perspective of a related set of concerns [START_REF]IEEE Recommended Practice for Architectural Description of Software-Intensive Systems. IEEE Std[END_REF]. Applying a viewpoint reduces the scope of an entire model to only present the interesting parts from a particular perspective. Besides, the model transformation mechanism produces new models from other models in input. For example, an analytical model for a simulation purpose can be produced from a descriptive model.

The MBSE traceability mechanism ensures the navigation through the layers. It allows system engineers to guarantee the consistency and the completeness of the models.

As presented, MBSE methodologies present many advantages to cover the system engineering processes. Thanks to the viewpoint, the model transformation and the traceability mechanisms, they allow various stakeholders to collaborate to the entire system architecture during the project life-cycle.

Nevertheless, studies [START_REF] Mavor | Human-System Integration in the System Development Process: A New Look[END_REF] [8] highlight the fact that the human element is not adequately taken into consideration in systems design. They point the recurrent issues of the system engineering to usually address human considerations afterthought the system architecture has been already specified and designed. The ISO/IEC 15288 remains vague regarding human involvements within a system. Mavor and Pew [START_REF] Mavor | Human-System Integration in the System Development Process: A New Look[END_REF] lists the shortcomings that lead to human integration problems:

• Human factors integration trends and standards are not captured • Human performance metrics, targets, and limitations are not specified • Human role design, job design, and organizational design are insufficiently captured • Team activity and team requirements are insufficiently captured

In response to the challenges, the Human-System Integration (HSI) initiative has been proposed.

HUMAN-SYSTEM INTEGRATION

Human-System Integration (HSI) is defined as "a multidisciplinary approach of integrating human, technical and organizational aspects at a system level, with full considerations given to the human requirements" [START_REF] Booher | Handbook of human systems integration[END_REF]. HSI frameworks give human factors practitioners the means to define human characteristics and constraints. Besides, they give systems engineers the tools to integrate these considerations into systems design. HSI extends tradeoff evaluations and agreements comparatively to the technical, organizational and human concerns.

HSI covers the following domains: human factors engineering, system safety, health, environment, habitability, survivability, training, personnel and manpower. These domains reflect how the human component impacts on the system performances and reciprocally, how a system influences the human component. HSI needs to address both human capabilities and system design within the system engineering processes.

Figure 1 illustrates the HSI model in terms of inputs and outputs [START_REF] Booher | Handbook of human systems integration[END_REF]. HSI process is implied during different stages of the system development (systems definition, systems development and systems deployment). It also focuses on the user needs during all the system design steps (highly concentrated user focus) and exploits disciplines related to human studies (human related technologies & disciplines).

The HSI process output is the system designed (technology), integrating human aspects (people and organization).

Figure 1 HSI Process [1]

Although relevant in many aspects, several studies highlighted major limitations to effectively integrate HSI in the system engineering processes.

Firstly, studies [START_REF] Mavor | Human-System Integration in the System Development Process: A New Look[END_REF] [9] identified a lack of shared methods, tools and formalisms precluding a meaningful communication and collaboration across system engineers, HSI practitioners and other personnel implied in system engineering. Current HSI tools do not take sufficiently into account the system development process. Without these means, system design efforts will continue to be inconsistent, incomplete and redundant since these domains are not able to collaborate.

Secondly, [START_REF] Orellana | Human System Integration Ontology: Enhancing Model Based System Engineering to Evaluate Human-System Performance[END_REF] asserts that the integration of HSI within system engineering is not just a technical issue, but also cultural and organizational. New semantics are needed to enlarge the modeling semantics scope to the HSI semantics. This new semantics will allow human elements to be considered in the whole system perspective.

THE NEED FOR AN EXTENDED FORMALISM

As previously explained, one major limitation of HSI and system engineering disciplines to collaborate is the lack of formalism adapted.

In one hand, SysML modeling language is standardized and widely adopted in the system engineering field. In the other hand, HSI discipline deals with a large variety of nonstandardized formalisms, each framework providing its own.

SysML provides the necessary means to encompass the HSI formalisms by extending its meta-model and including dedicated profiles. Several projects focus on the SysML extension to include HSI formalisms [10] [11] [12] [START_REF] Prun | Extending eFFBD formalism to task model[END_REF].

One of the most interesting features of SysML is its ability to support the simulation of the designed system. The existing solutions to simulate a system regarding both technical and human aspects are currently restricted.

Traditionally, system engineers use simulation tools to analytically test the system behavior and properties according to the requirements and specifications. Besides, human factors engineers developed models [START_REF] Paterno | Model-based design and evaluation of interactive applications[END_REF] [15] [START_REF] Lebiere | Multitasking and cognitive workload in an ACT-R model of a simplified air traffic control task[END_REF] which simulate human behaviors according to different points of view: fatigue, stress, cognitive workload, performance etc.

The sticking point is that there are very few links between system engineering and human factors models for simulation. They remain fragmented on their objectives, semantics and tools. The possibilities of interfacing them remain poorly studied. This constitutes obstacles to fulfill the requirements and the quality of the system.

A great effort needs to be done to implement the model transformation mechanism, to transform descriptive models to analytical models allowing the simulation objectives.

In response to these challenges, the U.S. Army Research Lab developed the Improved Performance Research Integration Tool (IMPRINT) to model human-system tasks in order to analyze functional allocation and workload estimation [START_REF] Allender | Modeling human performance: Impacting system design, performance, and cost[END_REF]. The IMPRINT tool provides descriptive and analytical models. The descriptive models are comparable to the Activity Diagram SysML diagram since it defines the tasks process fall to the system and the actors. The analytical models run simulation based on the descriptive models to analysis the task allocation in terms of workload levels.

The IMPRINT approach is interesting in a MBSE context since it uses models to evaluate a design architecture in terms of technical and human factors points of view.

Efforts must be done in this direction because it is not possible yet to fully integrate IMPRINT in a systemic models supported by MBSE.

THE NEED FOR AN EXTENDED METHODOLOGY

As Orellana [START_REF] Orellana | Extending Model-Based System Engineering for Human Machine Interaction Analysis and Fault Tolerant Design[END_REF] explained: "extending the modeling language is just a part of the extension needed. To extend MBSE to the specialty system engineering fields, not only do engineers need common semantics, but a way to integrate these semantics into current methods, processes, and tools being used by specialty system engineers".

This assertion highlights the need of a methodology to support system management and system design integrating technical and human aspects.

New methodologies are needed to guide the collaboration enhancement between system engineering and HSI disciplines.

First of all, an organizational plan needs to be established for monitoring both system engineering and HSI disciplines to ensure an effective collaboration. Muralidhar [START_REF] Muralidhar | How human systems integration and systems engineering can work together[END_REF] warns against the fact that one discipline prevails over the other one. Typically, system engineers tend to relegate aside human considerations aside because they do not know how to integrate these aspects into the system design. Developing a close relationship between system engineers and HSI practitioners is an important leverage to ensure the development of safe and effective systems. Dedicated supervisors should be allocated to manage the risk of disparities of concern in the design process.

Secondly, it is necessary to manage the systemic HSI process (Figure 1). The MBSE methodologies provide a favorable material to support it. They allow engineers to model the whole system, from the requirements analysis to the physical integration. Within this frame, HSI should focus on integrating their efforts on these models. The extended formalism and the traceability mechanism allow them to inspect the human considerations consistency and completeness through the entire system design.

Then, the methodology needs to be extended to guarantee that the Verification and Validation processes (V&V) consider human capabilities and constraints. The V&V activities are critical to ensure the overall quality of the system. Verification refers to the process of determining whether the system fulfill the specified requirements.

Validation refers to the evaluation of the design applied complies with the functional and the performance requirements. Simulation, as described early, can assist system designers to perform these process.

CONTRIBUTIONS

The objective of this study is to highlight the lack of formalism and methodology precluding a meaningful collaboration between system engineering and HSI. In this context, it is difficult to design a system guaranteeing the mandatory criteria of performance and safety within modern complex interactive systems.

Our research is guided by the need of a system engineering methodology adapted fora the integration of human factors issues, paying particularly attention to take these aspects into account across all the system layers. We consider that MBSE provides a relevant frame to support this methodology. As presented, MBSE guides the modelling activities through all the system engineering process and ensures the completeness and the consistency of the system design.

Then, we regard SysML as an adapted formalism to support our future methodology. Already widespread in the system engineering community, SysML provides the means to integrate HSI modelling semantics thanks to the model transformation, profile and viewpoint mechanisms. An extended SysML will facilitate the involvement of human factors practitioners within system design.

Our research will be proved on an aeronautical maintenance study case. Maintenance is central to navigation safety since it prevents aircraft failures. This activity represents an important financial cost for airlines, and it is implied in flights delays, cancellation and also accidents. Maintenance activities involve a high level of cognitive tasks within a stressful environment. Operators are potentially subject to high stress and high work load levels which deteriorate their performance and reliability.

Our methodology proposes to:

• Carry on the extension of SysML with HSI to better integrate human aspects. In our case, a maintainer can be defined by its experience, its physical and cognitive tasks involved in its activity, its stress, its fatigue, its relation with its organization and so on. • Trace these human aspects through all the MBSE layers. For example, the system used by the maintainers needs to be adapted to the current state of its user to ensure its performance. Then, the system could provide a dynamic functional allocation mechanism according to these aspects. Such an implementation impacts simultaneously all the system architecture layers (operational, functional, component and physical). • Extend simulation models to include human aspects.

As presented, the IMPRINT tool provides an interesting approach to evaluate human factors aspects within a system. Based on it, we will carry on the implementation of a model transformation mechanism to produce analytical models from descriptive models to simulate simultaneously technical and human aspects. This effort will allow to trade the system off to elicit the most adapted system design. • Develop a methodology to supervise the HSI efforts within the system engineering processes.


François a,b, , Guillot, Jérémy a,b, , Pesneau, Pierre Clautiaux

Exact approaches for solving a covering problem with capacitated subtrees

Keywords: Covering problems, Column generation, Dantzig-Wolfe decomposition, Branch-and-bound, Dynamic programming

In this document, we present a covering problem where vertices of a graph have to be covered by rooted subtrees. We present three mixed-integer linear programming models, two of which are compact while the other is based on Dantzig-Wolfe decomposition. In the latter case, we focus on the column generation subproblem, for which we propose several algorithms. Numerical results are obtained using instances from the literature and instances based on a real-life districting application. Experiments show that the branch-and-price algorithm is able to solve much bigger instances than the compact model, which is limited to very small instance sizes.

Introduction

In this document, we present a covering problem that arises in the application field of solid waste collection. Our problem belongs to the family of districting problems. The term "districting" refers to the operation of regrouping small geographical areas, called basic units, in bigger groups, called districts. In this work, basic units are presented as vertices of a graph, and districts are rooted subtrees.

For fifty years, researchers have been studying this family of problems for political districting [START_REF] Hess | Nonpartisan political redistricting by computer[END_REF][START_REF] Mehrotra | An optimization based heuristic for political districting[END_REF], sales territory design [START_REF] Zoltners | Sales territory alignment: A review and model[END_REF][START_REF] Ríos-Mercado | A reactive grasp for a commercial territory design problem with multiple balancing requirements[END_REF][START_REF] Elizondo-Amaya | A dual bounding scheme for a territory design problem[END_REF] and the optimization of public services [START_REF] Hanafi | Municipal solid waste collection: An effective data structure for solving the sectorization problem with local search methods[END_REF][START_REF] Muyldermans | District design for arc-routing applications[END_REF] (see, for example, Kalcsics' recent literature survey [START_REF] Kalcsics | Districting problems[END_REF]). In any districting problem, districts are required to be "connected", "compact" and "balanced". An important issue in this field of research is to formally define these three concepts, whose meaning strongly depends on the industry they apply to.

The connectivity constraint implies that any pair of vertices belonging to the same district can be connected by some path whose vertices also belong to that district. This constraint is usually meant to prevent intersecting districts, but it can also be used to represent geographical obstacles. In [START_REF] Ríos-Mercado | A reactive grasp for a commercial territory design problem with multiple balancing requirements[END_REF], the authors used the original connectivity constraint, which leads to linear programs with an exponential number of constraints. In [START_REF] Zoltners | Sales territory alignment: A review and model[END_REF] and [START_REF] Mehrotra | An optimization based heuristic for political districting[END_REF], the authors proposed a more restrictive version of connectivity, which uses the notion of center. For each district, one must set a center j; when the problem is defined on a graph, the center corresponds to a vertex. The connectivity constraint is verified if every district is a subtree of a spanning tree rooted in j. In [START_REF] Zoltners | Sales territory alignment: A review and model[END_REF] and [START_REF] Mehrotra | An optimization based heuristic for political districting[END_REF], the authors used a shortest path tree to define the spanning tree, and thus the connectivity constraint. The use of the shortest path tree is motivated in [START_REF] Mehrotra | An optimization based heuristic for political districting[END_REF] by the fact that a connected and compact sector is likely to contain the shortest path from the center to every node it contains. We use the same approach in this document.

Out of the three concepts (connectivity, compactness, balance), compactness is the one that is most open to interpretation. It is generally part of the objective function. Compactness can be modeled in a p-median problem fashion [START_REF] Zoltners | Sales territory alignment: A review and model[END_REF][START_REF] Mehrotra | An optimization based heuristic for political districting[END_REF], where one minimizes the sum of the distances between any vertex and the center of the district it belongs to, or in a p-center problem fashion [START_REF] Ríos-Mercado | A reactive grasp for a commercial territory design problem with multiple balancing requirements[END_REF][START_REF] Elizondo-Amaya | A dual bounding scheme for a territory design problem[END_REF], where one only takes into account the maximum distance, among all districts, between a vertex and its corresponding center. Our approach would be closer to the p-center model since we consider, for each district, the maximum shortest path distance between a vertex and its center. In this work, that distance is called the radius of the district. However, we define the total cost of a solution as the sum of the district radiuses instead of the maximum district radius. In this regard, our model can be compared to that of [START_REF] Muyldermans | District design for arc-routing applications[END_REF] and [START_REF] Silva De Assis | A redistricting problem applied to meter reading in power distribution networks[END_REF]. One of the criteria of [START_REF] Muyldermans | District design for arc-routing applications[END_REF] was the sum of the nor-malized radiuses for all districts, while [START_REF] Silva De Assis | A redistricting problem applied to meter reading in power distribution networks[END_REF] used the sum of the district diameters, i.e. the maximum distance between any two vertices belonging to the same district.

The balance constraint implies that the difference of volume between any two districts must be small, with regards to workload, service time or population size, for example. In most cases, this constraint requires, for each attribute being considered, that the volume of any district is close to the mean, by a given margin. This constraint can also be represented by an objective function minimizing the distance between each district and the mean [START_REF] Silva De Assis | A redistricting problem applied to meter reading in power distribution networks[END_REF].

In our problem, districts are not subject to a minimum workload. In the real-life instances we studied, balance was not an industrial requirement. However, workload is bounded from above, which prevents the districts from having an excessive workload. The number of districts to compute is given, and each district is subject to two capacity constraints. Another specificity of our problem is that it consists in finding a vertex cover in the graph instead of a partition. In cases where some vertices belong to many shortest paths, it helps ensuring that a feasible solution is computed.

Several papers on districting problems describe heuristics based on local search [START_REF] Hanafi | Municipal solid waste collection: An effective data structure for solving the sectorization problem with local search methods[END_REF][START_REF] Ríos-Mercado | A reactive grasp for a commercial territory design problem with multiple balancing requirements[END_REF][START_REF] Silva De Assis | A redistricting problem applied to meter reading in power distribution networks[END_REF]. There are works based on mathematical programming as well, and more specifically Dantzig-Wolfe decomposition [START_REF] Mehrotra | An optimization based heuristic for political districting[END_REF]. Column generation methods are well-suited for this family of problems, since they can easily be decomposed in one partitioning/covering master problem and a subproblem for each district.

The aim of our work is to solve exactly this problem using algorithms based on mathematical programming. In particular, we focus on a branch-and-price method. The main difficulty is to design an efficient subproblem oracle which, for a given center, computes the subtree rooted in this center that minimizes the reduced cost. This subproblem can also be seen as an extension of the knapsack problem with precedence constraints, two capacity constraints, and a profit that depends on the maximum distance from the center (radius). We propose several methods for solving this problem. The first method is based on executing, for all relevant radius values, a two-dimensional generalization of the branch-and-bound algorithm proposed in [START_REF] Cho | The critical-item, upper bounds, and a branch-and-bound algorithm for the tree knapsack problem[END_REF]. The second method solves a dynamic program, based on the algorithm proposed in [START_REF] Cho | A depth-first dynamic programming algorithm for the tree knapsack problem[END_REF], that includes the radius in its recurrence relation.

We validate our methods numerically using real-life instances and instances from the vehicle routing literature. Our experiments show that the model based on column generation exhibits much better performances than the compact models. Moreover, the enumerative branch-andbound oracle seems to be the best method to solve the column generation subproblems, although the dynamic programming algorithm can be faster in a significant number of cases.

In section 2, we present the problem formally, along with the different models studied in this document. We dedi-cate section 3 to the column generation subproblem while section 4 describes the techniques we use in the branchand-price algorithm. We analyze the numerical results in section 5 before offering some brief concluding remarks and suggestions for future research.

Problem definition and mathematical models

In this section, we formally define the problem we are studying. For this purpose, we disregard the industrial application that motivates this work and focus on a vertex covering problem restricted to subtrees.

Formal problem definition

We represent the geographical area to be serviced using a weighted directed graph G = (V, A, c) in which each vertex i ∈ V represents a waste collection point and each arc (i, j) ∈ A represents a connection from i to j of length c i,j . Let n = |V |. Any vertex i ∈ V is also associated with two values w 1 i ∈ N and w 2 i ∈ N corresponding, respectively, to the time required to service that point and the volume of waste to be collected therein. The fleet is composed of K vehicles, each of them limited to a maximum total service time W 1 ∈ N and a capacity W 2 ∈ N of waste that can be collected.

Each district that we want to generate is defined by a pair (U, j), where U ⊆ V is a set of vertices, and j ∈ U is the district's center. For a given district centered at j, let d i,j be the distance of the shortest path in G from j to vertex i and let π j (i) be the vertex preceding i in this path. If there is no path from j to i in G, then d i,j = +∞ and we assume, by convention, that π j (i) = j. These shortest paths define a tree T j = (V, A j ) rooted in j with A j = {(π j (i), i), i ∈ V \ {j}}. The radius of a district (U, j) is defined by max i∈U d i,j . In addition, for a district (U, j), we have that i ∈ U =⇒ π j (i) ∈ U .

In the following, given T j = (V, A j ) a tree rooted in j ∈ V and i ∈ V \ {j}, we define T j (i) the subtree of T j rooted in i and composed of all the descendants of i. In addition, given a set of vertices U ⊆ V , we define T j [U ] = (U, A j ∩ (U × U )) the subtree of T j rooted in j and containing every vertex of U . We say that set U induces a valid subtree of

T j if j ∈ U and T j [U ] is a connected component of T j .
Definition 1 (valid district, radius). Given a tree T j = (V, A j ) rooted in j, three vectors d •,j , w 1 , w 2 in R n + , and two real numbers W 1 and W 2 , a valid district centered at j is a pair (U, j) such that T j [U ] is a valid subtree and i∈U w i ≤ W for = 1, 2. The radius of a valid district (U, j) is defined as r(U, j) = max i∈U d i,j .

Problem 1 (minimum cost cover with valid districts). Given V a set of vertices, two vectors w 1 and w 2 in N n , two integer values W 1 ≥ 1 and W 2 ≥ 1, an integer K, and for each vertex j ∈ V , a tree T j of shortest paths from j to any other vertex in V , and a vector d •,j ∈ R n of shortest path distances from j to any other vertex in V , find a cover of V using K valid districts with distinct centers and minimizing the sum of the districts' radiuses.

The problem of covering a graph with valid districts is NP-hard, since the bin packing problem is reducible to this covering problem with trees of depth 1, null distances and capacity W 2 equal to +∞.

Note that an item may belong to two different districts (and thus counts in two knapsack constraints) in a solution, which may be counter-intuitive. However, this is mandatory in order to keep the precedence constraints relevant. It also helps finding solutions in which few items are covered more than once.

Compact formulations

The problem of vertex covering using valid districts of minimum cost can be formulated as an integer linear program. Our first model considers natural decision variables (center selection, assignments, and radiuses). Let x ∈ {0, 1} n×n be a matrix of decision variables. For i = j, x i,j is equal to 1 if i is assigned to center j, 0 otherwise. Each variable x j,j is equal to 1 if j is the center of some district, and 0 otherwise. Let also r ∈ R n be a vector of variables. For each vertex j ∈ V , r j is equal to the radius of the district centered at j if j is a center, and 0 otherwise. The problem can then be modeled as follows:

(P0) =                                                              min j∈V r j s.t. j∈V xj,j ≤ K, j∈V xi,j ≥ 1 ∀i ∈ V, i∈V w i xi,j ≤ W = 1, 2, ∀j ∈ V, xi,j ≤ x π j (i),j ∀i ∈ V, i = j, ∀j ∈ V, r j ≥ di,j xi,j ∀i ∈ V, ∀j ∈ V, r j ∈ R + ∀j ∈ V, xi,j ∈ {0, 1} ∀i ∈ V, ∀j ∈ V. (1)
(

The objective function is equal to the sum of the radiuses. Note that although this is not explicitly enforced, when j is not a center, r j = 0 in any optimal solution. Constraint (1) is an upper bound for the number of districts used ; constraints (2) ensure that each vertex belongs to at least one district ; and constraints (3) guarantee that no district can exceed any of the two maximum capacities. Each constraint (4) asks that if i is assigned to cluster j, then its predecessor is assigned to j as well. Constraints [START_REF] Elizondo-Amaya | A dual bounding scheme for a territory design problem[END_REF] guarantee that each variable r j carries the corresponding radius value.

We now propose a better model for the problem, with different variables. For that purpose, we use a model similar to that of [START_REF] Elloumi | A new formulation and resolution method for the p-center problem[END_REF] for the p-center problem. For each vertex j ∈ V , let σ j be a permutation of the vertices in V in order of increasing distance from center j. For a given j ∈ V and i ∈ {1, . . . , n} \ {j}, let ∆ i,j be the difference between the distances (from center j) of vertex i and its predecessor in σ j . It can be defined as follows: ∆ j,j = 0 and ∆ σj (k),j = d σj (k),j -d σj (k-1),j for k = 2, . . . , n. We can now distribute the maximum distance in a district among all the vertices it contains. To this end, we introduce a new matrix y ∈ {0, 1} n×n of binary variables. For any valid pair (i, j), we have that y i,j is equal to 1 if j is a center and the maximum distance in district j is greater than or equal to d i,j , and 0 otherwise. Thus, if j is chosen as a center, and vertex σ j (i) defines the radius of district j, we have y j,j = y σj (2),j = . . . = y σj (i),j = 1 and y σ(i+1),j = . . . = y σ(n),j = 0.

The problem can then be modeled as the following MIP.

(P) =                                                                              min j∈V i∈V \{j} ∆i,j yi,j s.t. j∈V xj,j ≤ K, j∈V xi,j ≥ 1 ∀i ∈ V, i∈V w i xi,j ≤ W = 1, 2, ∀j ∈ V, xi,j ≤ x π j (i),j ∀i ∈ V, i = j, ∀j ∈ V, y σ j (k),j ≤ y σ j (k-1),j ∀k ∈ {2, . . . , n}, ∀j ∈ V, xi,j ≤ yi,j ∀i ∈ V, ∀j ∈ V, xi,j ∈ {0, 1} ∀i ∈ V, ∀j ∈ V, yi,j ∈ {0, 1} ∀i ∈ V, ∀j ∈ V. (8) (9) (10) (11) (12) (13) (14)
(15)

The model described hereinabove differs from (1)-( 7) by its objective function, which is computed from decision variables y, and constraints [START_REF] Elloumi | A new formulation and resolution method for the p-center problem[END_REF] and [START_REF] Ceselli | An optimization algorithm for a penalized knapsack problem[END_REF]. Constraints [START_REF] Elloumi | A new formulation and resolution method for the p-center problem[END_REF] ensure that variables y are consistent, and constraints [START_REF] Ceselli | An optimization algorithm for a penalized knapsack problem[END_REF] state that if i is in the district of center j, then the radius of district j is greater than or equal to d i,j .

Extended formulation

In this section, we show how to derive an extended formulation of (P) using Dantzig-Wolfe decomposition.

Given j ∈ V , we define G j the set of extreme points of the convex hull of the set of integer solutions to system (10)-( 15) associated with j. Note that a solution associated with center j can refer to an empty district, therefore not covering center j. We can rewrite compact model (P) by substituting vectors x and y with an integer combination of the elements of G j . Given g ∈ G j , let x g i,j (resp. y g i,j ) be the value of variable x i,j (resp. y i,j ) in g. Let λ j ∈ {0, 1} |G j | be a vector of binary variables such that λ g j is equal to the number of times extreme point g is used in the solution. The new model can be defined as follows.

(P M ) =                                                min j∈V i∈V ∆ i,j   g∈G j y g i,j λ g j   s.t.
j∈V g∈G j

x g j,j λ g j ≤ K, (κ)

j∈V g∈G j x g i,j λ g j ≥ 1 ∀i ∈ V, (φ)
g∈G j λ g j = 1 ∀j ∈ V, (α)
λ g j ∈ {0, 1} ∀g ∈ G j , ∀j ∈ V. ( 16
) ( 17
) ( 18
) ( 19
)
This problem is called the master problem. It has a polynomial number of constraints and an exponential number of variables. A classical way to solve this problem is to use a branch-and-price method. It relies on solving iteratively the continuous relaxation of (P M ). For this purpose, we first solve a restriction of the model to a subset of variables (called the restricted master problem). Then, we need to solve a subproblem to find a non-generated variable corresponding to a negative reduced cost. That subproblem can be decomposed into n subproblems (one for each center candidate). In practice, each subproblem consists in computing the district that corresponds to the smallest reduced cost, with regards to its center candidate. If none of the generated districts correspond to a negative reduced cost, then the current solution for the continuous relaxation of (P M ) is optimal. Otherwise, we add all n generated columns to the restricted master problem (even if the corresponding reduced cost is non-negative).

Let us associate dual variables κ ∈ R -, φ i ∈ R + for all i ∈ V and α j ∈ R for all j ∈ V respectively to constraints [START_REF] Della Croce | New exact approaches and approximation results for the penalized knapsack problem[END_REF], [START_REF] Ibaraki | A dynamic programming method for single machine scheduling[END_REF] and [START_REF] Wentges | Weighted dantzig-wolfe decomposition for linear mixed-integer programming[END_REF]. We then get, for each center candidate, the following pricing problem.

(P j ) =                                                        min i∈V ∆i,j yi,j -κ xj,j - i∈V φi xi,j -αj ⇔ -αj -max κ xj,j + i∈V φi xi,j - i∈V ∆i,j yi,j s.t. i∈V w i xi,j ≤ W = 1, 2, xi,j ≤ x π j (i),j ∀i ∈ V, i = j, y σ j (k),j ≤ y σ j (k-1),j ∀k ∈ {2, • • • , |V |}, xi,j ≤ yi,j ∀i ∈ V, xi,j ∈ {0, 1} ∀i ∈ V, yi,j ∈ {0, 1} ∀i ∈ V. ( 20
) (21) (22) (23) (24) (25)
Considering φ i as the profit gained from covering a vertex with some district, this problem amounts to searching the valid district that minimizes its radius while maximizing the profit induced by the covered vertices.

In order to get an integer solution of (P M ), we apply a branch-and-price method. The details of this method are presented in a subsequent section. In the following section, we focus on solving the continuous relaxation of model (P M ). The restricted master problem is solved using a linear programming solver. The main difficulty is to solve the column generation subproblems ; therefore, we propose several algorithms that can solve them efficiently.

Methods for solving the column generation subproblem

In this section, we focus on the pricing problem introduced in the previous section. This problem can be decomposed into n pricing subproblems, each of them using a distinct center vertex j in V . Thus, for each vertex j in V , we want to generate an extreme point g in G j such that the reduced cost of λ g j , the associated decision variable in (P M ), is minimized. To simplify notations, we will drop index j henceforth, since we are always considering a subproblem centered at some vertex j. We will also assume, w.l.o.g., that the subproblem is a maximization problem. Consequently, the profit associated with each node is non-negative in the pricing subproblem. However, when applying the Lagrangian relaxation of the capacity constraints, negative profits may also arise. The problem can be formally stated as the following.

Problem 2 (valid district maximizing the sum of profits minus the radius). Given a tree T = (V, A) rooted in j, two vectors w 1 and w 2 in N n , two vectors d, and p in R n + , and two integers W 1 ≥ 1 and W 2 ≥ 1, find a valid district (U, j) maximizing i∈U p i -r(U, j). This problem is a generalization of two variations of the knapsack problem that have been studied in the literature. The first is the so-called Tree Knapsack Problem, defined by a single knapsack constraint and no radiuses (or, equivalently, radiuses equal to 0). In [START_REF] Cho | The critical-item, upper bounds, and a branch-and-bound algorithm for the tree knapsack problem[END_REF], a branch-and-bound algorithm is proposed to solve this problem while a dynamic programming approach is discussed in [START_REF] Cho | A depth-first dynamic programming algorithm for the tree knapsack problem[END_REF]. Both methods make use of a depth-first search (DFS) ordering of the items. The second variation is the so-called Penalized Knapsack Problem, which considers a single knapsack constraint and no precedence constraints. In this context, the radiuses are called penalties. In [START_REF] Ceselli | An optimization algorithm for a penalized knapsack problem[END_REF], this problem is solved using a dynamic program wherein the items are arranged in order of increasing penalty. However, the precedence constraints change the structure of our problem, preventing us from using a similar ordering.

We designed two solvers for the subproblem wherein the center vertex is given. The first is based on enumerating the relevant r(U, j) values. We then iteratively solve a variation of the knapsack problem with precedence constraints. The second solver directly computes both the set of covered vertices and the district radius simultaneously. Note that, since the Knapsack Problem is weakly NP-hard and the complexity of the dynamic program described in section 3.2 is in pseudo-polynomial time, Problem 2 is weakly NP-hard.

Subproblem solving by radius enumeration and

branch-and-bound One way to solve the subproblems is to enumerate all potential radiuses in order to obtain a variation of a classical problem from the literature called Tree Knapsack Problem (TKP).

Problem 3 (2D knapsack problem with tree-precedence constraints (TKP-2D)). Given a tree T = (V, A) rooted in j, a vector p in R n + , two vectors w 1 and w 2 in N n , and two integers W 1 ≥ 1 and W 2 ≥ 1, find a valid district (U, j) that maximizes i∈U p i .

There are, at most, n relevant radius values for each center. For each such value r, we can preprocess the instance data to remove each vertex whose distance from the center is greater than r. Given center j and radius r, the subproblem amounts to solving TKP-2D(T j ; r) -r + α j , where TKP-2D(T j ; r) is defined as follows.

(TKP-2D(T j ; r)) =                            max i∈V p i x i s.t. i∈V w i x i ≤ W = 1, 2, x i ≤ x πj (i) ∀i ∈ V, i = j, x i = 0 ∀i ∈ V, d i > r x i ∈ {0, 1} ∀i ∈ V (26) (27) (28) (29)
The problem described hereinabove has already been studied in one dimension in [START_REF] Cho | A depth-first dynamic programming algorithm for the tree knapsack problem[END_REF][START_REF] Cho | The critical-item, upper bounds, and a branch-and-bound algorithm for the tree knapsack problem[END_REF]. A dynamic programming algorithm [START_REF] Cho | A depth-first dynamic programming algorithm for the tree knapsack problem[END_REF] and a branch-and-bound algorithm [START_REF] Cho | The critical-item, upper bounds, and a branch-and-bound algorithm for the tree knapsack problem[END_REF] have been proposed to solve it. These algorithms have been extended to two knapsack constraints in [START_REF] Clautiaux | Clustering problem for waste collection[END_REF].

The branch-and-bound method that we propose is an extension of the algorithm from [START_REF] Cho | The critical-item, upper bounds, and a branch-and-bound algorithm for the tree knapsack problem[END_REF], initially addressing the 1D knapsack problem with tree-precedence constraints.

In [START_REF] Cho | The critical-item, upper bounds, and a branch-and-bound algorithm for the tree knapsack problem[END_REF], the branching scheme of the branch-and-bound algorithm consists in selecting a vertex i in V and constructing two sub-problems: a first subproblem in which node i does not belong to the solution, and a second subproblem in which node i does belong to the solution. The dual bound for each node of the branching tree is obtained by applying a Lagrangian relaxation of the knapsack constraint and performing a linear search to compute the best Lagrange multiplier. Let us first formalize the problem obtained after relaxing the capacity constraint.

Problem 4 (Maximum weight subtree). Given a tree T = (V, A) rooted in j, and a vector p in R n , find a valid subtree T [U ] that maximizes i∈U p i .

If μ ∈ R + is the Lagrange multiplier associated with the knapsack constraint, then the Lagrangian relaxation corresponds to the maximum weight subtree problem with pi = p i -μw i where w i is the weight of item i in the knapsack constraint.

This problem can be solved in linear time [START_REF] Cho | The critical-item, upper bounds, and a branch-and-bound algorithm for the tree knapsack problem[END_REF] using a dynamic program. A state of this program is defined by some vertex i ∈ V and corresponds to a maximum subtree of T (i). The dynamic program (DP) can then be written as

β(i) = max{0, pi + k:π(k)=i β(k)},
and the optimal solution is constructed from β(j). Using an inverse depth-first search (DFS) ordering to walk the vertices of T , all states of the dynamic program are computed iteratively.

To solve the Lagrangian problem (that is, to find the best multiplier μ), one performs a linear search on the value of μ. Note that the objective function of the Lagrangian problem is convex and piecewise linear. The linear search starts with a null Lagrange multiplier, and the maximum weight subtree associated with this multiplier. Increasing the value of μ will linearly decrease the value of the Lagrangian function (assuming the current solution is not already feasible) while keeping the same solution until, for some vertex i * of the solution, profit β(i * ) becomes zero (a break point of the piecewise linear function is reached). Vertex i * is called a critical item. The new solution is simply obtained by removing subtree T (i * ) from the previous solution. The procedure is repeated until a solution satisfying the knapsack constraint is obtained. The value of μ used to obtain the last breaking point is the optimal solution of the Lagrangian problem.

The sequence of critical items obtained along the solution of the Lagrangian problem is used to define the successive branching decisions in the branch-and-bound tree.

In the two dimensional case, we consider the same branching scheme and evaluation method. However, as the Lagrangian relaxation is obtained by penalizing two knapsack constraints, the linear search is no longer relevant to find the best Lagrange multipliers. Consequently, we use a descent algorithm to optimize them heuristically. In order to fully take advantage of the branching strategy used in [START_REF] Cho | The critical-item, upper bounds, and a branch-and-bound algorithm for the tree knapsack problem[END_REF], based on a sequence of critical items, the descent algorithm should iteratively remove subtrees rooted in critical items, similarly to the linear search in the one dimensional case.

Let us first formalize the problem obtained after relaxing both capacity constraints. Given multipliers µ ∈ R 2 , the problem is again equivalent to the maximum weight subtree problem with pi = p i -µ • w i for all i ∈ V and can be solved with the same dynamic program. Now, by restricting the descent algorithm to directions that only increase the Lagrange multipliers, we have that each successive solution is obtained by removing subtrees from the previous one. This property follows from the fact that the values pi and β(i) decrease proportionally to µ. The optimal solution remains identical until at least one of the DP state is such that β(i * ) = 0 for some vertex i * of the solution (which defines the length of the step made in the chosen direction), in which case we remove the subtree rooted in the critical vertex i * .

The sequence of critical vertices induced by the descent algorithm is used as the successive branching decisions in our branch-and-bound algorithm.

Finally, let us describe the descent direction that we use. The most straightforward direction corresponds to the violation of the capacity constraints by the current solution. However, preliminary tests showed that moving along only one capacity constraint at a time, namely the most violated one, leads to better performances.

When (TKP-2D(T j ; r)) is solved for a given value r, we obtain a primal bound for (P j ), which can be used to prune branch-and-bound nodes for subsequent radius values (sometimes even at the root node). We also obtain a dual bound for any radius r < r by keeping the same set of vertices and updating the objective function by replacing penalty r with r . The value obtained is a valid upper bound since it is the value of an optimal solution of a relaxation of the problem associated with radius r (some vertices of larger radius can still be used, which relaxes the problem).

Considering this, we first focus on the primal bound by evaluating the "expected" best radius first, i.e. the radius that is the most often associated with an optimal solution of subproblem (P j ), before evaluating the other radiuses in decreasing order to produce useful dual bounds.

Direct subproblem solving with dynamic programming

While our branch-and-bound algorithm is based on an enumeration of the relevant radiuses, the dynamic program (DP) presented hereafter solves a subproblem by computing both the set of covered vertices and the district radius simultaneously. First, we present a DP designed for the case where the radius is given, which is an extension of the one presented in [START_REF] Cho | A depth-first dynamic programming algorithm for the tree knapsack problem[END_REF] in two dimensions. Then, we propose a method to efficiently take the radius into account in the DP's recurrence, as well as several techniques that improve the algorithm's performances.

Dynamic program for the TKP-2D

This dynamic program was originally proposed in [START_REF] Cho | A depth-first dynamic programming algorithm for the tree knapsack problem[END_REF] in one dimension as a variant of [START_REF] Johnson | On knapsacks, partitions, and a new dynamic programming technique for trees[END_REF], then extended to two dimensions in [START_REF] Clautiaux | Clustering problem for waste collection[END_REF]. We present the dynamic program using our own formalism. In the following, we consider a tree T = (V, A) rooted in j. W.l.o.g., we assume that the vertices of V are indexed according to some depth-first search (DFS) ordering. Under this assumption, if item k is selected, the next choice is related to item k + 1 (the next in the DFS ordering). If item k is not selected, no vertices in T (k) can be selected. In this case, the next vertex to be considered is η

(k) = min{i : i > k, i ∈ T (k)}.
If k is the last vertex in the DFS ordering, or if all indices k > k correspond to descendants of k, then we assume, by convention, that η(k) = n + 1. Figures 1a and1b represent the relation between the precedence tree and the successorship structure in the dynamic program. Building on this precedence structure, we can define the states of the dynamic program by pairs (k, w), where k ∈ {1, . . . , n + 1} represents the index in the DFS order, and w ∈ {0, . . . , W 1 } × {0, . . . , W 2 } is a vector representing resource consumptions. Since profits are computed using dual values, they are not integer. Thus, only DP based on weights are viable.

In the following, for any two vectors u and v of the same dimensions, u ≤ v and u + v are componentwise operations. Let W = W 1 W 2 . Then, the dynamic program can be defined as follows.

ψ CS (k, w) =          0 if k = n + 1, max{p k + ψ CS (k + 1, w + w k ), if k ≤ n ψ CS (η(k), w)} and w + w k ≤ W , ψ CS (η(k), w) otherwise,
where ψ CS (k, w) represents the value of the best solution limited to vertices in {k, . . . , n + 1} and a resource consumption less than or equal to w. The optimal solution corresponds to ψ CS (1, 0 0 ).

Note that the capacity constraints are included in the structure of the graph. Therefore, solving this dynamic program is equivalent to finding a longest path in a directed acyclic graph, where the vertices are states of the program and the arcs represent the options available in the recurrence.

This directed acyclic graph has O(nW 1 W 2 ) vertices and arcs. A longest path from vertex (1, 0 0 ) can be computed in O(|V |W 1 W 2 ) time and space. The algorithm of [START_REF] Cho | A depth-first dynamic programming algorithm for the tree knapsack problem[END_REF] is a direct application of Bellman's algorithm on this graph, in the case of a single capacity constraint. On the one hand, using a DFS ordering is motivated by the precedence constraints, which can then be integrated in the recurrence relation. On the other hand, the best DPbased methods for the penalized knapsack problem [START_REF] Della Croce | New exact approaches and approximation results for the penalized knapsack problem[END_REF] consider the items in order of increasing penalty, so that the penalty itself does not require adding an explicit dimension to the state space. However, these two orderings are very likely to be contradictory. For example, consider an instance with four items {i 0 , i 1 , i 2 , i 3 } such that π(i 1 ) = π(i 2 ) = i 0 , π(i 3 ) = i 2 , and r(i 0 ) < r(i 2 ) < r(i 1 ) < r(i 3 ). Then, it is clear that there does not exist a DFS ordering such that the items are arranged in order of increasing radius (penalty). For this reason, the dynamic program in the next section will be based on a DFS ordering of the items.

Extension of the dynamic program to include district radiuses

In this section, we show how to extend the dynamic program defined hereinabove to include district radiuses in the computation of the solution cost. The dynamic program can then be directly used to solve (P j ). For the sake of clarity, let r i be the radius corresponding to vertex i.

The basic idea of the following DP is to modify ψ CS in such a way that it takes into account the current radius of the partial solution. Now, a state has the form (k, w, r),

where k and w are the same as in ψ CS , and r is the radius of the current partial solution. When item k is considered, the radius does not change if k is not selected. If k is selected, the new radius is udpated to value r k if r k > r.

In the latter case, the cost of the state will also be updated to account for the penalty entailed by the new radius.

Given x ∈ R, let (x) + = max{0, x}. If x is a vector, then (x) + applies componentwise. We assume, w.l.o.g., that the vertices are indexed according to some DFS ordering.

ψ R (k, w, r) =                  0 if k = n + 1, max p k + (r k -r) + + ψ R (k + 1, if k ≤ n w + w k , max{r, r k }), and w + w k ≤ W , ψ R (η(k), w, r) ψ R (η(k), w, r) otherwise.
In this dynamic program, ψ R (1, 0 0 , 0) corresponds to an optimal solution for problem (P j ). Since the considered graph has no oriented cycles, we use a label-setting algorithm to solve this problem, i. We now need to define a dominance relation that takes the radius into account. Let i and j be two labels such that k( i ) = k( j ). Label i dominates label j if and only if at least one of the following dominance relations is met. Clearly, i dominates j if i bears a higher profit than j while consuming less resource and being closer to the district center.

Dominance relation 1. Label i dominates label j if p( i ) ≥ p( j ), r( i ) ≤ r( j ), and w( i ) ≤ w( j )
The second dominance relation implies that i dominates j if i bears a better profit and balances its radius while consuming less resource than j . Dominance relation 2. Label i dominates label j if p( i ) ≥ p( j ), p( i ) -r( i ) ≥ p( j ) -r( j ), and w( i ) ≤ w( j ). Proof. If dominance relation 2 is verified and r( i ) ≤ r( j ), then dominance relation 1 is verified as well, and i dominates j .

Let us then assume that r( i ) > r( j ). i is a better label than j , since its net profit p( i ) -r( i ) is higher than that of j , while its resource consumption is lower. In order to justify maintaining j , we have to assume that it will be combined with a set of vertices I. Let ⊕ be the operation that adds a set of vertices to some label. Since k( i ) = k( j ) and w( i ) ≤ w( j ), we can also combine I and i , then compare the resulting new labels ( i ⊕ I) and ( j ⊕ I). Since the net profit of these two new labels is p( i ) + p(I) -max{r( i ), r(I)} ≥ p( j )+p(I)-max{r( j ), r(I)} , then ( i ⊕I) is still better than ( j ⊕ I). This proves that i dominates j . Dominance relation 2 is at least as strong as 1 and does not involve heavier computations. For this reason, it is natively included in our algorithm.

The validity of the dynamic program builds on the fact that it considers the vertices iteratively according to a DFS ordering based on the precedence graph. Any such ordering guarantees the consistency of indices k+1 and η(k), for any vertex k ∈ V . Note that there exists an exponential number of such DFS orderings. Among these orderings candidates, we break ties by selecting the one in which vertices are also arranged in order of decreasing radius induced by their subtree.

For specific instances, using the concept of subtree radiuses in a DFS ordering leads to better performances in practice for the DP. An explanation for this behavior is that the DP states are associated with their final radius earlier in the recurrence relation, thus generating labels with higher costs that can be eliminated earlier.

Lagrangian filtering method

The dynamic program described in the previous section induces a pseudo-polynomial number of states. The number of labels that are required to solve the program can be huge in some cases, since the number of dimensions is a limiting factor in the algorithm's ability to remove dominated labels. The objective of our filtering method is to detect state transitions that can never be associated with an optimal solution. For this purpose, we use a technique, derived from the SSDP method (see [START_REF] Ibaraki | A dynamic programming method for single machine scheduling[END_REF], for example), that solves a relaxation, filters out some transitions, and then reintroduces the constraints that were relaxed. In our implementation of this method, we use a Lagrangian relaxation of the capacity constraints.

Given ζ ∈ R 2 a Lagrange multiplier vector, we compute the Lagrangian profit pk of each item k, that is, pk = p k -ζ •w k . The recurrence can then be defined as follows.

ψ LR (k, r) = max pk + (r k -r) + + ψ LR (k + 1, max{r, r k }), ψ LR (η(k), r) (30)
Let us now consider the graph corresponding to DP (30). Each node of this graph is a state of the DP. The arcs correspond to the potential decisions from each state (selecting an item or not). Figure 2 represents an example of such a graph, where the arc costs are defined as follows. Each vertex (k, r) has two outgoing arcs. Arc a = (k, r), (k + 1, max{r, r k }) , represents selecting item k in the solution, and its cost pa is pk + (r k -r) + . Arc a = (k, r), (η(k), r) represents dismissing this item, which does not induce any additional cost, i.e. pa = 0.

The complexity of this DP is in O(|V | 2 ) in time and space, since each of the n vertices is associated with at most one state per relevant radius.

Let P B be some primal bound for the value of an optimal solution. For example, that bound could have been inferred from a feasible solution. Given a vector of Lagrange multipliers, we compute, for each vertex (i, r) in the state graph, ω -((i, r)) and ω + ((i, r)) defined, respectively, as the value of a longest path from vertex (1, 0) to vertex (i, r) and the value of the longest path from vertex (i, r) to vertex (n + 1, 0). Given arc a = ((i, r), (j, r )), if ω -((i, r)) + pa + ω + ((j, r )) < P B, then a is not associated with an optimal solution. Note that this bound is still valid in the original space with the capacity constraints. In practice, when we filter out the outgoing arc of vertex (k, r) that corresponds to selecting item k, we will not consider solutions constructed from label (p, k, w, r) for any values p and w. The same applies to outgoing arcs that correspond to disregarding a item.

Lagrange multipliers ζ are computed using a standard subgradient algorithm. Let pt ∈ R n be the vector of Lagrangian profits, that is, pt k = p k -ζ t • w k , and let x t be the optimal solution obtained after iteration t of the subgradient algorithm. The associated subgradient di-

(1, 0) (2, 0) (3, 2) (4, 2) (4, 4)
(5, 0) ( rection ∇ t is defined by the capacity constraints viola-

tion, ∇ t = w 1 • x t -W 1 w 2 • x t -W 2 +
and the step size s t follows Polyak's step size,

s t = 2 • 0.1( pt •x t ) ∇ t 2 2
, where 0.1( pt • x t ) is an approximation of the optimality gap between the current multipliers ζ t and the optimal multipliers ζ * . The algorithm terminates after two non-improving iterations, i.e. iterations where pt • x t ≤ pt • x t-1 .

Completion bounds: filtering out labels according to

their cost Each label (p, k, w, r) of the dynamic program corresponds to a partial solution that is feasible with regards to the capacity constraints and the precedence constraints of the pricing subproblem. In a dynamic program, the main way to remove labels is to check dominance relations. However, it is also possible to remove labels using dual bounds, as we would in a branch-and-bound algorithm. Given some primal bound P B and some label , filtering out label is valid if p( ) + DB( ) ≤ P B, where DB( ) corresponds to some dual bound on the value that can be gained from to some terminal state, thus completing label .

The strength of this test depends on the quality of dual bound DB( ). In order to save computation time, the method that we use recycles computations from the Lagrangian filtering method that was described previously. Let ζ ∈ R 2 + be a vector of Lagrangian multipliers. Given label defined by tuple (p( ), k( ), w( ), r( )), let be its projection, by relaxation of the capacity constraints, defined by triple (p( ) -ζ • w( ), k( ), r( )), and let v ζ ( ) be the largest value of a path from the vertex bearing to sink (n + 1, 0). A dual bound DB( ) on the best solution that can be constructed from is given by DB

( ) = v ζ ( ) + ζ • (W -w( )). Given a pair (k( ), r( ))
and a residual capacity W r = Ww( ), we could compute optimal Lagrange multipliers ζ * Wr . However, a less time-consuming approach consists in reusing ζ * the multipliers that were obtained at the end of the subgradient algorithm described in the previous section. This technique will be called "completion bound" in the computational experiments.

4. Speeding up the convergence of the branch-andprice algorithm

Generating initial columns

In essence, our column generation algorithm is based on the primal simplex algorithm. One of the bottlenecks of this algorithm is to find an initial feasible basis. In order to address this issue, one generally uses an artificial basis, made from a set of artificial columns that do not exist in the original linear program but whose cost is prohibitive, so that they cannot be part of an optimal basis with a positive coefficient. This phase might require a significant number of column generation iterations before removing all artificial columns from the basis. Thus, by generating a set of random columns that do not take reduced costs into account -since the latter are biased by the cost of the basic artificial columns -, we are able to reduce the time spent in phase one of the simplex algorithm.

In order to generate such columns, we apply the following greedy heuristic for every center vertex. We start from the empty solution, i.e. the solution that does not cover any vertex. If there are no vertex that can be appended to the current solution, stop the method. Otherwise, add a vertex that is 1) the closest to the center, 2) not covered yet, and 3) induces a new solution that still verifies both capacity constraints. This algorithm generates districts that cover vertices close to their centers. It can be adapted to generate sparse districts instead, which adds more variety to the set of initial columns.

Branch-and-price

In order to obtain an integer solution for the master problem (P M ), we use a branching algorithm. In each branching node, we solve the continuous relaxation of problem (P M ), to which we add new constraints representing the branching decisions.

The branching decisions are based on the variables x i,j from compact model (P). Given vertices i and j, adding constraint x i,j ≥ 1 (resp. x i,j ≤ 0) is equivalent to adding constraint g∈G j x g i,j λ g j ≥ 1 (resp. g∈G j x g i,j λ g j ≤ 0) to problem (P M ). As a consequence, a new dual variable, associated with this branching constraint, is added to the objective function of subproblem (P j ). Note that (P j ) might still generate columns that do not verify the branching constraints depending on the dual values of the constraints it violates. In this case, the column is added to the restricted master problem, and may only enter the basis with a null value. The algorithm branches on the variable x i,j whose value is still very fractional, i.e. close to 0.5, after solving the parent node. Finally, the branching nodes are explored according to a Best-first search ordering. Another way of branching would be to apply the branching decisions to the subproblem directly. This could be done by cleaning up the master problem and modifying the trees in the subproblems. We chose to branch in the master problem instead, since this method entails no modification in the subproblem (it just changes the cost of the arcs). As such, we do not have to address inconsistency issues such as branching decisions implying that for some i, j, k ∈ V with i predecessor of k in T j , x i,j = 0 while x k,j = 1, which would make subproblem (P j ) infeasible.

Column generation is known to suffer from convergence issues, known as the tailing-off effect. In order to address this weakness, we use the stabilisation technique from [START_REF] Wentges | Weighted dantzig-wolfe decomposition for linear mixed-integer programming[END_REF] that substitutes the dual solution for (P M ) by a smoothed dual solution. Let (κ, φ, α) t be the dual solution at iteration t, and let (κ, φ, α) be the dual solution associated with the best lower bound known for (P M ) (see below for the computation of this bound). Then, the smoothed solution that we use in the pricing subproblems is given by the following convex combination: (κ, φ, α) t = γ(κ, φ, α) + (1 -γ)(κ, φ, α) t . We use the implementation described in [START_REF] Pessoa | Automation and combination of linear-programming based stabilization techniques in column generation[END_REF], which automatically tunes parameter γ during the column generation process.

The so-called Lagrangian lower bound is based on the Lagrangian relaxation of constraints ( 8) and ( 9) in (P) where dual values κ and φ are respectively used as Lagrangian mutipliers for ( 8) and [START_REF] Silva De Assis | A redistricting problem applied to meter reading in power distribution networks[END_REF]. When this relaxation is applied, each center is now related to an independent supbroblem. Therefore, for any values κ and φ of proper sign, L(κ, φ) = Kκ + i∈V φ i + j∈V min g∈G j { j∈V i∈V ∆ i,j y g i,jκ j∈V x j,j -i∈V φ i j∈V x g i,j } is a valid lower bound for an optimal solution of (P). Now, let OP T (RM P ) be the optimum of the current reduced master problem and RC(SP j ) be the best reduced cost of a variable generated by subproblem j at the current iteration. By strong duality, OP T (RM P ) = Kκ * + i∈V φ * i + j∈V α * j , where κ * , φ * and α * are optimal dual values for RMP. The best reduced cost related to center j is RC(SP j ) = min g∈G j { j∈V i∈V ∆ i,j y g i,j -

κ * j∈V x g j,j -i∈V φ * i j∈V x g i,j -α * j }. Computing OP T (RM P ) + j∈V RC(SP j ), one obtains L(κ * , φ * ) + j∈V α * j - j∈V α * j = L(κ * , φ * ).
Therefore, OP T (RM P ) + j∈V RC(SP j ) is a valid lower bound for the problem.

Computational experiments

The purpose of our experiments is twofold: to empirically determine which oracle works best in the context of column generation, and to validate the hypothesis that the extended formulation leads to better performances than the original compact formulation.

Instances solved and hardware configuration

Two sets of instances are used to measure the performances of our methods. Instance set exeo is composed of 13 instances built from real-life data provided by a French company, which develops and publishes software solutions for route optimization and geolocation of household waste trucks. These instances range from 31 to 245 collection points. Instance set cvrp was randomly generated using 16 CVRP (Capacitated Vehicle Routing Problem) instances from the TSPLIB. Each instance from the TSPLIB was used to generate 25 new instances for our problem, resulting in a benchmark of 400 instances ranging from 7 to 262 collection points.

In order to generate several instances from one TSPLIB instance, we added some random noise to the arc lengths as follows. Given (i, j) an arc of length c i,j , we define the length c i,j = c i,j + r, with r a random real number generated according to a uniform distribution in [0, (c i,j ) 1.5 ]. In addition, we generated the second capacity constraint by keeping the same capacity (right-hand side) and shuffling the resource consumptions according to a random permutation.

Overall, the item resource consumptions can go as high as 4 100 and the district capacities range from 3 to 14 400, the average capacity being around 1 800.

The experiments were performed on PlaFRIM (Plateforme Fédérative pour la Recherche en Informatique et Mathématiques), using Dodeca-core Haswell Intel R Xeon R E5-2680 v3 @ 2,5 GHz processors and 128Go of RAM.

The linear programs and compact mixed-integer programs were solved using IBM ILOG CPLEX Optimization Studio 12.6.0 [START_REF]Ibm ilog cplex optimization studio cplex[END_REF]. We used the branch-and-price implementation from the BaPCod 2.0.0 framework [START_REF] Vanderbeck | Bapcod -a branch-and-price generic code[END_REF], which manages the branch-and-price search tree and provides a built-in stabilization.

Analysis of the branch-and-bound oracle

Our first numerical experiments focus on the branchand-bound method that we use to solve the column generation subproblem. Our objective is to find the best parameters for this algorithm. Specifically, there are two components that can impact the performances of the algorithm in a significant way: 1) the method used to solve the Lagrangian relaxation and 2) the order in which the radiuses are considered.

These results are described in Table 1. First, we compare two descent directions. Direction PROP corresponds to the capacity constraint violation projected onto R 2 + , while direction STEEP only increases the penalty associated with the most violated capacity constraint. Then, we combine the best of these two directions with a special ordering for the radiuses. In this ordering, we first solve the radius that is most frequently associated with an optimal solution for the pricing subproblem before proceeding with the remaining radiuses in order of decreasing distance. In the following, this last configuration is referred to as STEEP+BR. For each configuration and for each instance set, we report time the average time in seconds spent solving an instance to optimality (limited to two hours), nodes/cg the average number of branching nodes that were generated by the oracle during that time and solved the number of instances that were solved in less that two hours.

It transpires from these experiments that the direction we apply has a significant impact on the algorithm performances, since direction STEEP requires 43% less branching nodes in average compared to direction PROP. The improved radius ordering allows us to further reduce the number of branching nodes by another 80%. This reduction comes from the fact that this particular ordering allows us to quickly find a good primal bound for the subproblem that is then used to filter out other radiuses. In the next sections, the results for the branch-and-bound oracle correspond to configuration STEEP+BR.

Analysis of the dynamic programming oracle

In this section, we compare several parameters for the dynamic program in the context of column generation. The most important is the choice of the methods to eliminate states of the dynamic program, namely dominance relations, Lagrangian filtering, and completion bounds.

The default strategy is referred to as base. In this configuration, we use dominance relation 1 and we use neither Lagrangian filtering nor completions bounds. Furthermore, the items are ordered according to an arbitrary DFS ordering. Building on this configuration, we can enrich the algorithm by using the following options: dom corresponds to applying dominance relation 2 ; filt corresponds to the Lagrangian filtering method described in section 3.2.3 ; sort refers to a DFS ordering of the items that takes into account the radiuses of their subtrees ; options sgCB and stCB are associated with Lagrangian bounds computed from the Lagrangian relaxation of the capacity constraints. The difference between these last two is that sgCB uses the Lagrangian multipliers obtained at the end of the Lagrangian filtering method, while stCB sets one multiplier to 0 and computes the optimal multiplier for the other dimension. In strategy all, we enable every option.

We also compare the impact of removing one of these features from all: w/o dom use options filt, sort, sgCB and stCB but does not use dom.

The performances associated with each configuration are recorded in Table 2. For each instance set, we report time the average time in seconds spent solving an instance to optimality (limited to two hours), labels/cg the average number of labels of the dynamic program that were generated during that time and solved the number of instances that were solved in less that two hours.

The experiments confirm that all proposed techniques significantly improve the performances of the dynamic programming algorithm. In particular, two configurations are standing out: configuration all is best suited to solve instances from exeo set, while configuration w/o sort is more adapted to cvrp instances. In both configurations, the algorithm generates about 90% less labels than configuration base.

In the remainder, results for the dynamic programming oracle correspond to configurations all (referred to as DP1) and w/o sort (referred to as DP2).

Comparison of the subproblem oracles

We now compare both approaches for solving the column generation subproblem that we described previously: branch-and-bound and dynamic programming. To do so, we use two different experimental setups. In the first set of experiments, we compare the performances branchand-price algorithm when combined with one of the two oracles. However, this comparison is biased, as the subproblems solved by both oracles might be different: when the subproblem supports several equivalent solutions, the oracles might generate different columns, which changes the behavior of the branch-and-price algorithm. As a consequence, we propose a second way to compare the oracles. We first saved the pricing instances solved during the branch-and-price algorithm while using the branch-andbound oracle. Then, we solved every pricing instance independently using one of the oracles, outside of the branchand-price scheme. Table 3 and figures 3-14 summarize these computational experiments.

In Table 3, we report time the average time in seconds spent solving (P) using the branch-and-price algorithm, t MP the average time in seconds spent solving (P M )as opposed to solving the subproblems (P j ) -, and solved the number of instances that were solved in less that two hours.

When we focus on the branch-and-price performances, the branch-and-bound oracle is more efficient than the dynamic programming oracles, since it is 1.7 times faster on the exeo instances and can solve 21 additional cvrp instances. However, both methods exhibit performances that have the same order of magnitude for our test instances.

We now offer a more accurate comparison of the oracles' performances. On the one hand, figures 3-5 and 9-11 represent the number of pricing instances that were solved by each oracle under a given time limit. For example, in Figure 3, the curve corresponding to the branch-and-bound oracle shows that 90% of the pricing instances are solved in less than 0.1 time unit, while the remaining 10% are solved under 0.4 time units. In comparison, 20% of the pricing instances require at least 0.2 time units when we use a dynamic programming based oracle. Thus, the higher the curve, the more efficient the algorithm.

Based on these experiments, we can infer that the branch-and-bound oracle solves generally more instances than any dynamic programming oracle, given some time limit. Note, however, that the dynamic program seems more robust for solving pricing instances related to cvrp/eilB76 (Figure 10).

On the other hand, figures 6-8 and 12-14 represent the solving time ratios between a dynamic programming oracle and the branch-and-bound oracle. We can then see the proportion of instances for which the dynamic programming oracle is faster than the branch-and-bound method, as well as the associated factor of speed. For example, Figure 6 shows that the dynamic programming oracles are faster than the branch-and-bound oracle more than 75% of the time, and can be as much as 9 times faster than the branch-and-bound oracle. By contrast, the branch-andbound oracle is at least 5 times faster for solving about 20% of the pricing instances.

As in the previous experiments, the branch-and-bound oracle outperforms the dynamic programming method, except for the pricing instances based on exeo/075 (Figure 6) and cvrp/eilB76 (Figure 13). In Appendix A, we reported the diagrams associated with additional experiments using the same framework.

In the following experiments, we will only use the branch-and-bound oracle for solving the pricing subproblem.

Comparison of the compact formulations and the extended formulation

In this set of experiments, we compare the performances of our column generation method to those obtained by the compact formulations, which are solved using a commercial solver (CPLEX).

Table 4 compiles the performances of compact models (P 0 ) and (P) and the extended formulation. In this table, we report time tot the total time in seconds to solve an instance to optimality (or -if the algorithm did not terminate in less than 2 hours), time root the time spent solving the root of the branching tree, which corresponds to a continuous relaxation of the model. The optimality gap, referred to as gap opt, is ≈ 0% iff the algorithm did not terminate in less than 2 hours. If the algorithm did not even find any feasible solution, then this gap is equal to +∞. Column gap root indicates the quality of the continuous relaxation of each model compared to the best known primal solution. Note that we know the optimal value for each exeo instance except exeo/163 and exeo/171. Finally, nb bb nodes indicates the number of branching nodes that were evaluated, and time per node is the average time in seconds for each of these nodes. Regarding the cvrp instances, we only report the number of instances that were solved to optimality in less than 2 hours.

Clearly, the column generation algorithm can solve the instances more efficiently than the compact models combined with a generic solver. Even if evaluating a branching node is generally more time-consuming in the extended formulation, the quality of the continuous relaxation allows us to dramatically reduce the number of branching nodes that are required to solve the instances.

Note that, for instance exeo/245, the optimality gap is smaller at the end when we solve compact model (P) than when we solve the extended formulation. As a consequence, column generation methods seem less fit for this instance as the solving time for the restricted master problem is excessively high.

Conclusion

In this document, we studied a new covering problem where vertices of a graph are covered by rooted subtrees. We proposed several mathematical models. The most efficient model relies on column generation, for which we proposed different algorithms. To solve the pricing subproblem, the specialized branch-and-bound algorithm we developed is better on average than our state-of-the-art dynamic programming approach. However, there exists a significant number of instances of the subproblems for which the dynamic programming algorithm is faster.

In order to solve larger instances, valid inequalities could be added to the model in order to strenghten the dual bound in the branch-and-price algorithm. However this could significantly complicate the subproblem oracle, and a thorough study is needed to produce an efficient algorithm.

Several other variants of covering problems using trees are also interesting to consider, for examples problems where one would like to minimize the sum of distances between each pair of vertices of each cluster. A similar decomposition method can be applied, but the subproblem would be much more challenging to tackle. this appendix, we report the results of additional experiments comparing both pricing oracles. We separated these results from the main document since they are consistent with the ones we presented in subsection 5.4. In these experiments, we observed that the dynamic program can slightly outperforms the branch-and-bound oracle (see cvrp/eilB76 and cvrp/eilB101). In other cases, the oracles' performances are comparable, depending on what metric we use as a reference (see: exeo/174, cvrp/A76 and cvrp/C76 ; or exeo/75, exeo/102 and cvrp/gil262). In all remaining cases, the branch-and-bound clearly outperforms the dynamic programming approach (see: cvrp/D76 and exeo/109 ; or exeo/162 and exeo/87 ; or cvrp/att48 and cvrp/A101).

  A precedence graph (where node 8 is an artificial node bearing no profit nor weight) Successorship structure corresponding to the precedence graph in 1a. Arcs (k, k + 1) are drawn in black while arcs (k, η(k)) are grey.

Figure 1 :

 1 Figure 1: Relation between the precedence graph and the successorship structure.

  e. we consider the vertices of the graph in a topologic ordering, and keep a set of non-dominated partial solutions for each node. Nondominated solutions are represented by labels. A label is defined by tuple (p( ), k( ), w( ), r( )) containing the cost of the partial path from vertex (1, 0 0 ), the current position in the DFS ordering, the resource consumptions and the current district radius. The theoretical complexity of this new DP is O(|V | 2 W 1 W 2 ) in time and space.

Proposition 3 . 1 .

 31 Dominance relation 2 is correct.

Figure 2 :

 2 Figure 2: Graph representing the dynamic program induced by Lagrangian relaxation. States are defined by both a position in the DFS ordering and a radius. Bold arcs correspond to selecting the current item. Grey arcs correspond to not selecting the current item.

Figure 3 :

 3 Figure 3: exeo/075 walltimes

Figure 6 :Figure 12 :

 612 Figure 4: exeo/162 walltimes

Figure 14 :

 14 Figure 14: cvrp/eilD76 time ratios

Table 2 :Table 4 :

 24 cg solved time labels/cg solved time labels/cg solved time labels/cg solved time labels/cg solved time labels/cg solved time labels/cg solved cg solved time labels/cg solved time labels/cg solved time labels/cg solved time labels/cg solved time labels/cg solved time labels/cg solved Performances of the dynamic programming oracle time tot time root gap opt = (pb-db)/db gap root = (pb*-lp)Performances of the extended formulation Appendix A. Comparison of the branch-andbound oracle and dynamic programming oracle

Figure

  Figure A.15: exeo/087 walltimes

Table 1 :

 1 Performances of the branch-and-bound oracle

			PROP			STEEP			STEEP+BR
	Instance	time nodes/cg solved time nodes/cg solved time nodes/cg solved
	exeo/031	≈ 0	776	1/1	≈ 0	344	1/1	≈ 0	284	1/1
	exeo/032	≈ 0	46	1/1	≈ 0	46	1/1	≈ 0	19	1/1
	exeo/054	1	15	1/1	≈ 0	15	1/1	≈ 0	10	1/1
	exeo/075	67	51	1/1	70	51	1/1	65	31	1/1
	exeo/087	31	279	1/1	31	279	1/1	26	185	1/1
	exeo/102	50	109	1/1	50	109	1/1	48	71	1/1
	exeo/109	1 471	4 328	1/1 1 478	4 328	1/1	765	1 648	1/1
	exeo/125	7 200	1 860	0/1 7 200	1 859	0/1 7 200	720	0/1
	exeo/162	212	1 697	1/1	212	1 697	1/1	178	1 127	1/1
	exeo/163	7 200	888	0/1 7 200	893	0/1 7 200	418	0/1
	exeo/171	7 200	7 978	0/1 7 200	7 980	0/1 7 200	1 729	0/1
	exeo/174	7 200	8 000 257	0/1 1 288	8 963	1/1	635	1 998	1/1
	exeo/245	7 200	6 494 471	0/1 7 200	39 808	0/1 7 200	28 641	0/1
		time nodes/cg solved time nodes/cg solved time nodes/cg solved
	cvrp/att48	3 210	448 17/25 2 909	251 17/25 3 182	188	16/25
	cvrp/eil7	≈ 0	107	25/25	≈ 0	92	25/25	≈ 0	46	25/25
	cvrp/eil13	≈ 0	46	25/25	≈ 0	46	25/25	≈ 0	46	25/25
	cvrp/eil22	≈ 0	475	25/25	≈ 0	408	25/25	≈ 0	143	25/25
	cvrp/eil23	3	952	25/25	3	961	25/25	3	930	25/25
	cvrp/eil30	1	625	25/25	1	529	25/25	1	197	25/25
	cvrp/eil31	2	1 802	25/25	2	1 626	25/25	2	1 292	25/25
	cvrp/eil33	32	3 514	25/25	23	2 556	25/25	12	893	25/25
	cvrp/eil51	316	3 084	25/25	345	2 134	25/25	212	916	25/25
	cvrp/eilA101	6 810	22 557	1/25 7 018	10 782	1/25 7 200	4 016	1/25
	cvrp/eilA76	6 161	7 755	6/25 6 183	5 022	6/25 5 385	1 228	12/25
	cvrp/eilB101	6 970	15 405	2/25 6 915	11 631	2/25 6 717	3 117	2/25
	cvrp/eilB76	6 962	10 986	1/25 6 958	8 736	1/25 5 910	1 755	9/25
	cvrp/eilC76	5 148	6 438	9/25 4 914	3 407 11/25 4 651	1 524	11/25
	cvrp/eilD76	5 654	7 649	8/25 5 173	3 434	9/25 4 762	1 300	13/25
	cvrp/gil262	7 355	792 422	0/25 7 200	404 621	0/25 7 200	26 934	0/25
	Total (/400)			244			247			264

Acknowledgments

We would like to thank the anonymous referees for their valuable comments, which helped us to improve the quality of the manuscript.

Experiments presented in this paper were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI and IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux, CNRS and ANR in accordance to the programme d'investissements d'Avenir (see https://www.plafrim. fr/).

BB DP1 DP2

Instance time t MP solved time t MP solved time t MP solved exeo/031