COSYS, ESTAS Samuel Deleplanque Ifsttar

Martine Labbé

A Branch-Price-and-Cut Procedure for the Discrete Ordered Median Problem

The Discrete Ordered Median Problem (DOMP) is formulated as a set partitioning problem using an exponential number of variables. Each variable corresponds to a set of demand points allocated to the same facility with the information of the sorting position of their corresponding costs. We develop a column generation approach to solve the continuous relaxation of this model. Then, we apply a branch-price-and-cut algorithm to solve small to large sized instances of DOMP in competitive computational time.

Introduction

Logistics is a very active domain in contemporary Operations Research and Location Analysis is one of its main building blocks. In the last years, motivated by the need of applying more flexible models in Logistics, a new family of location models that takes into consideration the order of the costs to be minimized has emerged in the literature.

The ordered weighted average of a vector is equal to the weighted sum of its coordinates after ranking them by non decreasing order [START_REF] Nickel | Location Theory: A Unified Approach[END_REF]). When the ordered weighted average is used in a location problem, the vector contains the distances or allocation costs from clients to service facilities. Ordered median location problems were first introduced in networks and continuous spaces by [START_REF] Nickel | A unified approach to network location problems[END_REF] and [START_REF] Puerto | Geometrical properties of the symmetrical single facility location problem[END_REF], respectively. Later, they were extended to the discrete setting by [START_REF] Nickel | Discrete ordered weber problems[END_REF] and [START_REF] Boland | Exact procedures for solving the discrete ordered median problem[END_REF]. The Discrete Ordered Median Problem (DOMP) has been widely studied since the 90's and there is a number of different formulations, solution approaches and applications available in the literature. To cite a few, DOMP has been applied to discrete facility location in [START_REF] Boland | Exact procedures for solving the discrete ordered median problem[END_REF]; [START_REF] Marín | A flexible model and efficient solution strategies for discrete location problems[END_REF][START_REF] Marín | An extended covering model for flexible discrete and equity location problems[END_REF]; [START_REF] Nickel | Discrete ordered weber problems[END_REF]; [START_REF] Puerto | A new formulation of the capacitated discrete ordered median problem with {0, 1}-assignment[END_REF]; [START_REF] Puerto | Minimax regret single-facility ordered median location problems on networks[END_REF]; to location on networks in [START_REF] Nickel | A unified approach to network location problems[END_REF]; to hub network design problems in [START_REF] Puerto | Single-allocation ordered median hub location problems[END_REF][START_REF] Puerto | A specialized Branch & Bound & Cut for Single-Allocation Ordered Median Hub Location problems[END_REF][START_REF] Puerto | Ordered Median Hub Location Problems with Capacity Constraints[END_REF]; to determine values in cooperative game theory in [START_REF] Perea | Finding the nucleolus of any n-person cooperative game by a single linear program[END_REF]; to combinatorial optimization problems 1 INTRODUCTION 2 with ordering in [START_REF] Fernández | On Discrete Optimization with Ordering[END_REF][START_REF] Fernández | Ordered Weighted Average Combinatorial Optimization: Formulations and their properties[END_REF][START_REF] Fernández | Ordered Weighted Average Optimization in Multiobjective Spanning Tree Problems[END_REF], and to voting problems in [START_REF] Ponce | Mathematical programming formulations for the efficient solution of the k-sum approval voting problem[END_REF], etc. The reader is referred to the monographies by [START_REF] Nickel | Location Theory: A Unified Approach[END_REF] and [START_REF] Puerto | Ordered Median Location Problems[END_REF] for some other applications.

Given a set of clients, a set of candidate locations and assuming that the allocation costs of clients to facilities are known, DOMP consists in choosing p facility locations and assigning each client to a chosen facility with the smallest allocation cost in order to minimize the ordered weighted average of these costs. The ordered weighted average sorts the allocation costs in a non-decreasing sequence and then it performs the scalar product of this so-obtained sorted cost vector with a given vector of weights.

There are several valid formulations for DOMP that exploit specific features of the problem (see e.g. [START_REF] Boland | Exact procedures for solving the discrete ordered median problem[END_REF]; [START_REF] Marín | A flexible model and efficient solution strategies for discrete location problems[END_REF]; [START_REF] Labbé | A comparative study of formulations and solution methods for the discrete ordered p-median problem[END_REF] and the references therein). In [START_REF] Labbé | A comparative study of formulations and solution methods for the discrete ordered p-median problem[END_REF] a new formulation for DOMP has been proposed, based on a set packing approach, that is valid for general cost coefficients. This formulation gives rise to rather tight integrality gaps and was shown to be reasonably efficient to solve medium size instances when embedded in a branch-and-cut (B&C) scheme. In this paper we explore a different paradigm for solving DOMP based on an extended formulation using an exponential number of variables and that corresponds to a set partitioning model. Each variable represents a set of couples (client, position). These clients are served by the same facility and their position indicates the place of their allocation cost in the sorted list of allocation costs for the entire considered solution. To handle the exponential number of variables we use a column generation approach that is embedded in a branch-price-and-cut (B&P&C) algorithm. A recent similar approach to a different problem can be found in [START_REF] Doulabi | A constraint-programming-based branch-andprice-and-cut approach for operating room planning and scheduling[END_REF]. This scheme has never been applied to DOMP and it opens new avenues of research. Therefore, the contribution of this paper is to propose a new perspective in the resolution of DOMP based on formulations with an exponential number of variables and to develop an efficient B&P&C algorithm to handle them.

The remainder of this paper is organized as follows. Notation, models and algorithms are presented in Section 2. Subsection 2.2 introduces a new set partitioning formulation for DOMP. This formulation uses an exponential number of variables where each element of the partition is a set of clients that are assigned to the same facility together with their sorted positions. This formulation is theoretically compared with another valid formulation described in Subsection 2.1 and borrowed from [START_REF] Labbé | A comparative study of formulations and solution methods for the discrete ordered p-median problem[END_REF]. Subsection 2.3 describes the column generation algorithm that we have designed to overcome the large number of variables in the model. We prove that the pricing subproblem is solvable efficiently in polynomial time by using an ad hoc dynamic programming algorithm. We devote Section 3 to the implementation details of our B&P&C algorithm. We develop a GRASP heuristic, in Subsection 3.1, that is used to generate both a promising initial solution and a pool of variables to initialize the column generation routine. We also develop a stabilization routine, based on [START_REF] Pessoa | Exact Algorithm over an Arctime-Indexed Formulation for Parallel Machine Scheduling Problems[END_REF], that reduces considerably the number of iterations of the column generation approach in Subsection 3.2. In addition, Subsection 3.3 is devoted to an additional improvement, namely a preprocessing. The next two subsections, 3.4 and 3.5, present our branching strategies and some families of valid inequalities that will be added to the branch-and-price algorithm. In the last section, namely Section 4, we report on the final computational experiments. We evaluate the performance of the B&P&C algorithm and compare it to the compact formulation in Section 2.1. The paper ends with some concluding remarks.

Problem definition and formulations

Let I be a set of n points which at the same time represent clients and potential uncapacitated facility locations and let c ij denote the cost for serving client i's demand from facility j.

Given a set J of p open facilities, let c i (J) represent the cheapest cost for allocating client i to a facility in J, i.e. c i (J) := min j∈J c ij . Now let us sort the costs c i (J), i ∈ I by non-decreasing order of their values. The elements of the resulting vector of ordered costs are denoted by c (k) (J) and satisfy c (1) (J) ≤ • • • ≤ c (n) (J). We denote the set of all possible positions 1, . . . , n in this ordered vector by K.

Given vector λ = (λ k ) k∈K satisfying λ k ≥ 0, k ∈ K, the objective function of DOMP, is defined as

z(J) := k∈K λ k c (k) (J). ( 1
)
Recall that this objective function provides a very general paradigm to encompass standard and new location models. For instance, if λ 1 = . . . = λ n = 1 we obtain the median objective, if λ 1 = λ 2 = . . . = λ n-1 = 0, λ n = 1 we obtain the center objective, if (DOMP)

λ 1 = λ 2 = . . . = λ n-1 = α, λ n =
Several formulations of DOMP have been proposed in the literature using different types of variables. Among them we mention those based on a combination of the p-median and permutation polytopes [START_REF] Boland | Exact procedures for solving the discrete ordered median problem[END_REF]) or on covering approaches by using radius variables [START_REF] Puerto | A new formulation of the capacitated discrete ordered median problem with {0, 1}-assignment[END_REF], [START_REF] Marín | A flexible model and efficient solution strategies for discrete location problems[END_REF][START_REF] Marín | An extended covering model for flexible discrete and equity location problems[END_REF]).

An explicit formulation for DOMP: The Weak Order Constraints

In the following, we recall the Weak Order Constraints formulation, W OC, introduced in [START_REF] Labbé | A comparative study of formulations and solution methods for the discrete ordered p-median problem[END_REF], and that will be the starting point for the developments presented in this paper. This formulation uses two types of binary variables. Variable y j assumes value 1 if facility j ∈ I is open (i.e. j ∈ J) and 0 otherwise. Variable x k ij is equal to 1 if client i ∈ I is allocated to facility j ∈ I and the corresponding cost occupies position k ∈ K in the allocation cost ranking (i.e. c (k) (J) = c ij ). The choice of this formulation is motivated by its good performance in terms of integrality gap (see [START_REF] Labbé | A comparative study of formulations and solution methods for the discrete ordered p-median problem[END_REF]). However, it requests important memory space since it needs O(n 3 ) binary variables which may become prohibitive for moderate n.

We denote the rank of the allocation cost c ij by r ij , i.e. r ij = if c ij is the -th element in the list of the costs c ij , for all i, j ∈ I, sorted by order of non decreasing values and where ties are broken arbitrarily. For the sake of readability the reader is referred to Example 1 in Section 2.3. The formulation is as follows:

(WOC) min i∈I j∈I k∈K λ k c ij x k ij (2) s.t. j∈I k∈K x k ij = 1 i ∈ I (3) i∈I j∈I x k ij = 1 k ∈ K (4) k∈K x k ij ≤ y j i, j ∈ I (5) j∈I y j = p (6) i∈I j∈I      i ∈I j ∈I: r i j ≤r ij x k i j + i ∈I j ∈I: r i j ≥r ij x k-1 i j      ≤ n 2 k ∈ K, k = 1 (7) x k ij , y j ∈ {0, 1} i, j ∈ I, k ∈ K. ( 8
)
By means of (3) we ensure that each location is served by exactly one facility. In the same way, in each position there must be exactly one allocation cost (4). Constraints (5) translate the fact that a client can be allocated to a facility only if this facility is open and that the allocation cost of a client to a facility can be placed in at most one position. The equality constraint (6) implies that there are exactly p open facilities. Constraints (7), called weak order constraints, ensure that if client i is allocated to facility j and the corresponding costs c ij occupies the k-th position in the cost ranking of the solution then in the (k -1)-th position there must be a smaller allocation cost. This property is enforced by the coefficients of each variable in the inequality. In each constraint there are two different positions, k and k -1, so that, by (4), only two variables must take value one and all the others will be equal to zero. If we do not take into account the variables assuming the value zero and we assume that the variables with value one for positions k and k -1 correspond to allocation pairs in sorted position s and t, respectively, the inequality reduces to the following expression: (n 2 -(s -1))x k isjs + tx k-1 itjt ≤ n 2 , which is valid if and only if t < s. Finally, the variables are binary, see (8).

W OC can be reinforced by adding the following valid inequalities:

i ∈I j ∈I:

r i j ≤r ij x k i j + i ∈I j ∈I: r i j ≥r ij x k-1 i j ≤ 1, i, j ∈ I, k ∈ K, k = 1. (9)
Observe that constraints (7) are the aggregation over i, j ∈ I of inequalities (9). These inequalities are the so called strong order constraints, See [START_REF] Labbé | A comparative study of formulations and solution methods for the discrete ordered p-median problem[END_REF] for a detailed explanation.

A set partitioning formulation

From a linear programming relaxation point of view the above formulation is not the strongest one but it provides a good compromise between the number of required constraints and the quality of its linear relaxation bound, see [START_REF] Labbé | A comparative study of formulations and solution methods for the discrete ordered p-median problem[END_REF]. Further, it allows to solve to optimality problems of moderate size. One of its drawbacks is the use of a cubic number of variables, which can be prohibitive for large n. A second important problem of most known formulations for DOMP is their high degree of symmetry in case of allocation costs c ij or weighted (λ k ) with many ties. The reasons above motivate the introduction of a new formulation based on a different rationale. We observe that a solution for DOMP is a partition of the clients together with their positions in the sorted vector of costs so that each subset of clients in the partition is allocated to the same facility.

Let us consider sets of couples (i, k) where the first component refers to a client i and the second to a position k, namely S = {(i, k) : for some i ∈ I, k ∈ K}. Further, we denote by P(I × K) the family of all sets S for which all first (resp. second) coordinates of its couples are different.

Associated with each set S and facility j, we define a variable y j S equal to 1 if the set S is part of a feasible solution ((i, k) ∈ S iff x k ij = 1); and 0 otherwise. Let S be the set of couples whose first coordinate corresponds to the clients allocated to a given facility j of a feasible solution. The positions of these clients in the solution, i.e. the second coordinates of couples in S must be compatible with the ranking of all the allocation costs involved in the solution. Hence, they must, in particular, be compatible with the ranking of the costs c ij of the clients i allocated to j. This implies that for facility j ∈ J we only need to consider couple subsets S belonging to S(j) = {S ∈ P(I × K) : c ij ≤ c i j for all (i, k), (i , k ) ∈ S and k < k }.

Since in any feasible solution, each client i must be allocated to a unique facility j and its allocation cost must occupy a unique position k in the sorted list, the following relationship holds:

x k ij = S∈S(j):(i,k)∈S y j S , i, j ∈ I, k ∈ K. (10)
Next, we can evaluate the cost c j S induced by the set S provided that its clients are assigned to facility j in a feasible solution:

c j S = (i,k)∈S λ k c ij . (11)
To simplify the presentation in the following we denote by (i, •) any couple whose first entry is i regardless of the value of the second entry. Analogously, (•, k) denotes any couple whose second entry is k regardless of the value of the first entry.

The following valid formulation uses variables y j S and constitutes our Master Problem (M P ).

(MP) min j∈I S∈S(j) c j S y j S (12) s.t.

j∈I S∈S(j):

(i,•)∈S

y j S = 1 i ∈ I (13) j∈I S∈S(j): (•,k)∈S y j S = 1 k ∈ K (14) S∈S(j) y j S ≤ 1 j ∈ I ( 15
)
j∈I S∈S(j)

y j S ≤ p (16) n i=1 n j=1         S∈S(j): (i ,k)∈S r i j ≤r ij y j S + S∈S(j): (i ,k-1)∈S r i j ≥r ij y j S         ≤ n 2 k ∈ K, k = 1 (17) y j S ∈ {0, 1} S ∈ S(j), j ∈ I. (18)
The objective function (12) accounts for the sorted weighted cost of any feasible solution. Constraints (13) ensure that each client appears in exactly one set S. Constraints (14) ensure that each position is taken by exactly one client appearing in one set S. Constraints (15) guarantees that each facility j serves at most one set S of clients. Inequality ( 16) states that at most p facilities will be opened. By the following family of inequalities (17) we enforce the correct sorting of the costs in any feasible solution. Finally, the variables are binary.

One can relate M P and W OC. First, remark that for a given facility j, there is at most one cost c ij which occupies a given position k. Hence, the following constraints are valid for W OC:

i∈I x k ij ≤ y j , j ∈ I, k ∈ K. ( 19
)
Let W OC+ denote the formulation given by ( 2)-( 8) and ( 19) and consider the Dantzig-Wolfe reformulation of W OC+ in which constraints (5), ( 19) and ( 8) constitute the subproblem. The subproblem can be decomposed by facility.

On the one hand, the feasible points of the subproblem of a facility j correspond one to one to the sets S ∈ P(I × K). Hence, this Dantzig-Wolfe reformulation of W OC+ is given by the above master problem in which we consider variables y j S for all S ∈ P(I × K) (instead of only S ∈ S(j)). More precisely, the variables of W OC+ are related to the variables y j S through the following two equations:

x k ij = S∈P(I×K):(i,k)∈S y j S i, j ∈ I, k ∈ K and y j = S∈P(I×K)
y j S j ∈ I.

Moreover, constraints (13) correspond to constraints (3), constraints ( 14) to (4), constraints ( 16) to (6), and constraints (17) to (7). Finally, constraints (15) constitute the "convexity" constraints for the subproblems.

On the other hand, it is easy to see that the polyhedron of each subproblem, defined by constraints ( 5) and ( 19) together with x k ij ≥ 0 and y j ≤ 1, is integer. This implies that the linear relaxations of W OC+ and M P in which all sets S ∈ P(I × K) are considered provide the same bound. By restricting the subsets S to be considered for each facility j to belong to S(j), our formulation M P provides thus a stronger model. The computational experiments presented in Subsection 4.3 show that there exist instances for which the linear relaxation of M P provides a strictly better (higher) lower bound than the linear relaxation of W OC.

Formulation M P can be strengthened by adding valid inequalities borrowed from W OC. Indeed, one can translate valid inequalities (9) in terms of the y j S variables so that they can be used in the set partition formulation of DOMP. The translation of (9) results in:

S∈S(j): (i ,k)∈S r i j ≤r ij y j S + S∈S(j): (i ,k-1)∈S r i j ≥r ij y j S ≤ 1, i, j ∈ I, k ∈ K, k = 1. ( 20
)

Column generation to solve the linear relaxation of M P (LRMP)

Since the number of variables in M P is too large to be handled directly, in this section, we describe a column generation approach to solve it. Let (α, β, γ, δ, ) be the dual variables associated, respectively, to constraints ( 13), ( 14), ( 15), ( 16) and ( 17). The dual problem DP of LRMP is

(DP) max i∈I α i + k∈K β k - j∈I γ j -pδ - k∈K: k =1 n 2 k (21) s.t. i∈I: (i,•)∈S α i + k∈K: (•,k)∈S β k -γ j -δ - i ∈I j ∈I          (i,k)∈S: r i j ≥r ij k =1 k + (i,k)∈S: r i j ≤r ij k =n k+1          ≤ c j S j ∈ I, S ∈ S(j) (22) δ, γ j , k ≥ 0 j ∈ I, k ∈ K, k = 1. ( 23
)
In order to apply the column generation procedure, let us assume that we are given a set of columns that define a restricted Master Problem, and denote its linear relaxation by ReLRMP. This problem is solved to optimality and (α * , β * , γ * , δ * , * ) represents its optimal dual solution. See Example 1. The reduced cost, c j S , of column y j S , namely c j S = c j S -z j S is given by:

c j S = c j S + γ * j + δ * + i ∈I j ∈I        (i,k)∈S: r i j ≥r ij k =1 * k + (i,k)∈S: r i j ≤r ij k =n * k+1        - i∈I: (i,•)∈S α * i - k∈K: (•,k)∈S β * k . (24)
If c j S ≥ 0 for all j, S ∈ S(j) the current solution of ReLRMP is also optimal for the LRMP and the column generation procedure stops.

Otherwise, one has identified one (some) new column(s) to be added to the current reduced master problem to proceed further. In each iteration, ReLRMP and its reduced costs provide lower and upper bounds for the LRMP. Indeed it holds that (Desrosiers and Lübecke (2005))

z ReLRM P + p • min j∈I,S∈S(j) c j S ≤ z LRM P ≤ z ReLRM P , (25)
z ReLRM P + j∈I min S∈S(j) c j S ≤ z LRM P ≤ z ReLRM P , (26)
where z ReLRM P and z LRM P denote the optimal value of ReLRM P and LRM P respectively.

Example 1. Consider the following vector λ = (4, 2, 1), cost matrix C and precedence matrix R:

C =   1 3 6 3 1 8 6 8 1   , R =   1 4 6 5 2 8 7 9 3   .
For n = 3, there are 33 different sets of couples (i, k) in S. 1,1),(2,3),(3,2)} S 30 = {(1,2),(2,1),(3,3)} S 31 = {(1,2),(2,3),(3,1)} S 32 = {(1,3),(2,1),(3,2)} S 33 = {(1,3),(2,1),(3,2)}.

S 1 = {(1, 1)} S 2 = {(1, 2)} S 3 = {(1, 3)} S 4 = {(2, 1)} S 5 = {(2, 2)} S 6 = {(2, 3)} S 7 = {(3, 1)} S 8 = {(3, 2)} S 9 = {(3, 3)} S 10 = {(1, 1), (2, 2)} S 11 = {(1, 1), (2, 3)} S 12 = {(1, 1), (3, 2)} S 13 = {(1, 1), (3, 3)} S 14 = {(1, 2), (2, 1)} S 15 = {(1, 2), (2, 3)} S 16 = {(1, 2), (3, 1)} S 17 = {(1, 2), (3, 3)} S 18 = {(1, 3), (2, 1)} S 19 = {(1, 3), (2, 2)} S 20 = {(1, 3), (3, 1)} S 21 = {(1, 3), (3, 2)} S 22 = {(2, 1), (3, 2)} S 23 = {(2, 1), (3, 3)} S 24 = {(2, 2), (3, 1)} S 25 = {(2, 2), (3, 3)} S 26 = {(2, 3), (3, 1)} S 27 = {(2, 3), (3, 2)} S 28 = {(1, 1), (2, 2), (3, 3)} S 29 = {(
The sets S(j) are the following: We consider as initial pool of columns the variables y 1 18 and y 3 8 . With this set of variables, the ReLRMP is

S(1) = {S 1 , S 2 ,
(ReLRMP) min +2y 2 5 +10y 1 13 s.t. +y 1 13 ≥ 1 i = 1 +y 2 5 ≥ 1 i = 2 +y 1 13 ≥ 1 i = 3 +y 1 13 ≥ 1 k = 1 +y 2 5 ≥ 1 k = 2 +y 1 13 ≥ 1 k = 3 -y 1 13 ≥ -1 j = 1 -y 2 5 ≥ -1 j = 2 ≥ -1 j = 3 -y 2 5 -y 1 13 ≥ -2 -8y 2 5 -y 1 13 ≥ -9 k = 2 -2y 2 5 -3y 1 13 ≥ -9 k = 3 y ≥ 0 .
Actually, we are interested in its dual problem:

(DP) max +α 1 +α 2 +α 3 +β 1 +β 2 +β 3 -γ 1 -γ 2 -γ 3 -2δ -9 2 -9 3 s.t. +α 2 +β 2 -γ 2 -δ -8 2 -2 3 ≤ 2 (y ) +α 1 +α 3 +β 1 +β 3 -γ 1 -δ -2 -3 3 ≤ 10 (y ) α, β, γ, δ, ≥ 0 . Solving (DP), the solution is α 2 = 2, β 3 = 10, α 1 = α 3 = β 1 = β 2 = δ = 2 = 3 =
0 and the value of the objective function is f = 12.

Solving the pricing subproblem

Although any column y j S with negative reduced cost may be added to ReLRMP, we will follow a strategy that identifies the most negative reduced cost for each facility j. This approach may give rise to several candidate columns (multiple pricing, see [START_REF] Chvátal | Linear Programming (W. H. Freeman and Company)[END_REF]), which is advantageous for this procedure.

In order to do that, for each facility j ∈ I, we solve a subproblem to find the column y j S , S ∈ S(j), with minimum reduced cost. This set S must be such that there is at most one couple (i, •) for each client i and one couple(•, k) for each position k. Furthermore, the set S must enjoy that the allocation costs of its couples are compatible. We solve this problem by the following dynamic programming algorithm.

Let d k ij be the contribution of the pair (i, k) to the reduced cost of any column y j S such that (i, k) ∈ S. Depending on the values of k, d k ij is given by

d k ij =                            λ k c ij + i ∈I j ∈I: r i j ≤r ij k+1 -α i -β k if k = 1, λ k c ij + i ∈I n j ∈I: r i j ≥r ij k + i ∈I j ∈I: r i j ≤r ij k+1 -α i -β k if k = 2, . . . , n -1, λ k c ij + i ∈I j ∈I: r i j ≥r ij k -α i -β k , if k = n. (27)
Then, for a facility j, the problem of finding the variable y k S with minimum reduced costs can be formulated as: min

S∈S(j) c j S = γ * j + δ * + (i,k)∈S d k ij . (28)
Now for each facility j, we define a matrix D j as follows:

D j =      d 1 i 1 j d 2 i 1 j • • • d n i 1 j d 1 i 2 j . . . . . . d 1 inj d n inj      (29)
where i 1 , i 2 , . . . , i n is a permutation of the indices i = 1, . . . , n such that

c i 1 j ≤ c i 2 j ≤ • • • ≤ c inj .
Example 2. (Continuing from Example 1). We illustrate the procedure that computes the elements d k ij for all i, k = 1, . . . , n of the matrix D 1 . (j=1)

d 1 11 = λ 1 c 11 + r 11 2 -α 1 -β 1 = 4 d 2 11 = λ 2 c 11 + (n 2 -r 11 + 1) 2 + r 11 3 -α 1 -β 2 = 2 d 3 11 = λ 3 c 11 + +(n 2 -r 11 + 1) 3 -α 1 -β 3 = -9 d 1 21 = λ 1 c 21 + r 21 2 -α 2 -β 1 = 10 d 2 21 = λ 2 c 21 + (n 2 -r 21 + 1) 2 + r 21 3 -α 2 -β 2 = 4 d 3 21 = λ 3 c 21 + +(n 2 -r 21 + 1) 3 -α 2 -β 3 = -9 d 1 31 = λ 1 c 31 + r 31 2 -α 3 -β 1 = 24 d 2 31 = λ 2 c 31 + (n 2 -r 31 + 1) 2 + r 21 3 -α 3 -β 2 = 12 d 3 31 = λ 3 c 31 + +(n 2 -r 31 + 1) 3 -α 3 -β 3 = -4
Since r 11 < r 21 < r 31 the valid permutation is (1, 2, 3). This implies that

D 1 =   4 2 -9 10 4 -9 24 12 -4   i = 1 i = 2 i = 3
By using D j we obtain that a set S belongs to S(j) if and only if for every (i 1 , k 1 ) and

(i 2 , k 2 ) ∈ S such that i 1 < i 2 : k 1 < k 2 .
Our dynamic programming algorithm to obtain the minimum reduced cost for each j ∈ J builds upon this observation by constructing a solution to a reduced version of (28) in which only the first i l clients and the first k positions are considered.

For each couple (i l , k), we define a function:

g j (i l , k) = min{c j S : S ∈ S(j) and for all (i l , k ) ∈ S : i l ≤ i l and k ≤ k} (30)

and we denote an optimal solution of this restricted optimization problem by S j (i l , k).

Hence, the optimal value of problem ( 28) is equal to g j (i n , n) + δ + γ j and a corresponding optimal solution by S j (i n , n).

Our recursive procedure computes g j (i l , k) and S j (i l , k) for increasing values of l and k and exploits the following feasibility conditions on S:

(i) for each client i (resp. position k), at most one couple containing i (resp. position k) belongs to S.

(ii) if (i l 1 , k 1 ) and (i l 2 , k 2 ) ∈ S and k 1 < k 2 then r i l 1 j < r i l 2 j . More precisely, if (i l , k) belongs to S j (i l , k) then, from (i) it follows that g j (i l , k) = g j (i l-1 , k -1) + d k i l j . Otherwise, S j (i l , k) may contain a couple (i l , k ) with k ≤ k -1 or a couple (i l , k) with l ≤ l -1 but not both for otherwise condition (ii) would be violated. Hence, in this case, g j (i l , k) = min{g j (i l-1 , k -1), g j (i l , k -1), g j (i l-1 , k)}.
Combining the two cases, we obtain the following recurrence relation for l, k = 2, . . . , n:

g j (i l , k) = min{g j (i l-1 , k -1) + d k i l j , g j (i l-1 , k -1), g j (i l , k -1), g j (i l-1 , k)}. (31)
Obviously, if at the end of the procedure, g j (i n , n)+δ +γ j is negative the variable y j S j (in,n)

is a good candidate to be chosen in the next iteration of the column generation scheme. If we solve this problem for all j, we get c j R = min S c j S and if c j R < 0, we can activate (at least) y j R . Next, we solve a new reduced master problem ReLRMP with this (these) new activated variable(s).

Remark 1. Computing each matrix D j can be done in O(n 2 ). Next, obtaining g j (i n , n) requires the evaluation of the function g j (i, k) for all i ∈ I and k ∈ K. According to the above algorithm, the evaluation of each g j (i, k) is done in constant time. Solving the pricing subproblem amounts to evaluate g j (i n , n) for all j ∈ I. Therefore, the entire pricing subproblem can be solved in O(n 3 ) time.

Example 3. (Continuing from Example 2). We show the computation of the g j (i n , n) and S j (i n , n) for j = 1.

g 1 (i 1 , 1) = min{0, 4} = 0, S 1 (i 1 , 1) = ∅ g 1 (i 1 , 2) = min{2, 0} = 0, S 1 (i 1 , 2) = ∅ g 1 (i 1 , 3) = min{-9, 0} = -9, S 1 (i 1 , 3) = {(1, 3)} g 1 (i 2 , 1) = min{10, 0} = 0, S 1 (i 2 , 1) = ∅ g 1 (i 3 , 1) = min{24, 0} = 0, S 1 (i 3 , 1) = ∅ g 1 (i 2 , 2) = min{0 + 4, 0, 0, 0}, S 1 (i 2 , 2) = ∅ g 1 (i 3 , 2) = min{0 + 12, 0, 0, 0}, S 1 (i 3 , 2) = ∅ g 1 (i 2 , 3) = min{0 -9, 0, -9, 0}, S 1 (i 2 , 3) = {(1, 3)} g 1 (i 3 , 3) = min{0 -4, 0, -9, 0}, S 1 (i 3 , 3) = {(1, 3)}
We have obtained g 1 (i 3 , 3) = -9 and S 1 (i 3 , 3) = S 3 is the potential set to be used, since its reduced cost is negative. The corresponding reduced cost c 1 3 = g 1 (i 3 , 3) + δ + γ 1 = -9 + 0 + 0 = -9 < 0. Hence, we active variable y 1 3 .

Algorithm 1 Pricing Subproblem Algorithm.

1: g j (i1, 1) = min{0, d 1 i 1 j }; 2: if g j (i1, 1) = d 1 i 1 j < 0 then 3:

S j (i1, 1) = {(i1, 1)}; 4: else 5:

S j (i1, 1) = ∅; 6: end if 7: for k = 2, . . . , n do 8:

g j (i1, k) = min{d k i 1 j , g j (i1, k -1)}; 9: if g j (i1, k) = g j (i1, k -1) then 10: S j (i1, k) = S j (i1, k -1); 11: else 12: S j (i1, k) = {(i1, k)}; 13:
end if 14: end for 15: for l = 2, . . . , n do 16:

g j (i l , 1) = min{d 1 i l j , g j (i l-1 , 1)}; 17:

if g j (i l , 1) = g j (i l-1 , 1) then 18:

S j (i l , 1) = S j (i l-1 , k); 19: else 20:

S j (i l , 1) = {(i l , 1)}; 21:

end if 22: end for 23: for k, l = 2, . . . , n do 24:

g j (i l , k) = min{g j (i l-1 , k -1) + d k i l j , g j (i l-1 , k -1), g j (i l , k -1), g j (i l-1 , k)}; 25: if g j (i l , k) = g j (i l-1 , k -1) then 26: S j (i l , k) = S j (i l-1 , k -1); 27: else if g j (i l , k) = g j (i l , k -1) then 28: S j (i l , k) = S j (i l , k -1); 29:
else if g j (i l , k) = g j (i l-1 , k) then 30:

S j (i l , k) = S j (i l-1 , k); 31: else 32:

S j (i l , k) = S j (i l-1 , k -1) ∪ {(i l , k)}; 33:

end if 34: end for Next, the process continues with the following facilities, i.e. j = 2, 3. In this example the optimal solution can be certified after four complete iterations of the above process.

Dealing with infeasibility

One important issue when implementing a column generation procedure to solve a linear optimization problem is how to deal with infeasibility. This is specially crucial if the procedure is used within a branch-and-bound scheme to solve the linear relaxation of the problem at every node of the branching tree. To handle this, we resort to the so called Farkas pricing. This method was used previously, to the best of our knowledge, in [START_REF] Günlük | A Branch-and-Price Algorithm and New Test Problems for Spectrum Auctions[END_REF] and [START_REF] Ceselli | Optimizing the Cargo Express Service of Swiss Federal Railways[END_REF]. The term Farkas pricing was coined in [START_REF] Gamrath | Generic Branch-Cut-and-Price[END_REF].

According to Farkas' Lemma [START_REF] Farkas | A Fourier-féle mechanikai elv alkalmazásai[END_REF]), a reduced master problem is infeasible if and only if its associated dual problem is unbounded. Thus, to recover feasibility in the ReLRMP we have to revoke the certificate of unboundedness in the dual problem. This can be done by adding constraints to it. Since we are only interested in recovering feasibility in ReLRMP, one can proceed in the same way that for the usual pricing, but with null coefficients in the objective function of the primal. In this way, the Farkas dual problem is

max i∈I αi + k∈K β k - j∈I γj -pδ - k∈K: k =1 n 2 k (32) s.t. i∈I: (i,•)∈S αi + n k∈K: (•,k)∈S β k -γj -δ - i ∈I j ∈I        (i,k)∈S: r i j ≥r ij k =1 k + (i,k)∈S: r i j ≤r ij k =n k+1        ≤ 0 j ∈ I, S ∈ S(j) (33) δ, γj k ≥ 0 j ∈ I, k ∈ K, k = 1. ( 34
)
To identify new variables that make the reduced master problem feasible we use our dynamic programming approach in which we replace the column costs c j S by zeros. Farkas pricing is an important element in our approach because it allows to start the column generation algorithm with an empty pool of columns, although this is not advisable. Furthermore, Farkas pricing will be crucial in the branching phase to recover feasibility (whenever possible) in those nodes of the branching tree where it is lost after fixing variables.

A branch-price-and-cut implementation

In this section, we precise several components of the implementation of our set partitioning formulation based on a column generation approach. B&P&C is a branch-and-cut scheme that solves the linear relaxation at each node of the branching tree with the column generation algorithm previously described and may apply cuts to improve the obtained lower bound. (The reader is referred to [START_REF] Doulabi | A constraint-programming-based branch-andprice-and-cut approach for operating room planning and scheduling[END_REF] for another recent implementation of a B&P&C.)

Unless otherwise specified, to calibrate the best choice of the different parameters used in our B&P&C, we have performed a preliminary computational study based on a set of 60 instances with sizes n = 20, 30 and with a time limit of 1800 sec. Those are the smallest instances that we will eventually use in Section 4.

Upper bound for the Master Problem: A GRASP heuristic and an initialization stage

A heuristic algorithm that generates a good feasible solution for M P will provide a promising pool of initial columns as well as a good upper bound. GRASP [START_REF] Feo | A probabilistic heuristic for a computationally difficult set covering problem[END_REF], [START_REF] Feo | Greedy ramdomized adaptive search procedures[END_REF]) is a well-known heuristic technique that usually exhibits good performance in short computing time. In our case, it consists in a multistart greedy algorithm to construct a set of p facilities from a randomly generated set of facilities with smaller cardinality. Following [START_REF] Puerto | A modified variable neighborhood search for the discrete ordered median problem[END_REF] we have chosen, in a greedy manner, an initial set of p/2 facilities. Next, we improve this initial solution by performing a fixed number of iterations of a local search procedure.

The greedy algorithm adds iteratively a new facility to the current set of open facilities, choosing the one with the maximum improvement of the objective value. The local search consists in an interchange heuristic between open and closed facilities. The pseudocode of the GRASP used to solve the problem is described in Algorithm 1.

Algorithm 2 GRASP for DOMP.

1: Input(n, p, C, λ, n1, n2, q); 2: for n1 replications do 3:

PartialSolution ← ConstructRandomizedPartialSolution(q); 4:

Solution ← ConstructGreedySolution(PartialSolution); 5:

for n2 iterations do 6:

Solution ← LocalSearch(Solution); 7:

BestSolution ← UpdateSolution(Solution, BestSolution); 8:

end for 9: end for First of all, we would like to point out the remarkable behavior of the GRASP heuristic for this problem. In order to illustrate the appropriateness of our heuristic we have solved to optimality a number of instances of the problem (using the MIP formulation) to be compared with those given by our GRASP. In all instances, up to a size of n = 400, the solution provided by GRASP is always as good as the one obtained by the any of our MIP formulations with a CPU time limit of 7,200 seconds, see Section 4.

Moreover, it is not only advisable to use the GRASP heuristic because it provides a very good upper bound thus helping the exploration of the searching tree by pruning many branches of the branch-and-bound tree, but in addition, the construction phase of the heuristic also provides a very promising pool of initial columns for the B&P&C, in combination with the technique described in the following.

Since we are solving the linear relaxation of our master problem, LRMP, without generating its entire set of variables, using the primal simplex algorithm, the goal of the initialization phase is to find an initial set of columns that allows solving the M P by performing a small number of iterations in the column generation routine. We create variables using a modification of the local search routine of the GRASP algorithm. Every time that we find a promising feasible solution in the heuristic, we create the variables that define that solution (CreateSetVariables(J)). Algorithm 3 presents the pseudocode of this process.

Function CreateSetVariables(J) determines the costs involved in the solution, i.e. the minimum for each client among the open facilities. Then those costs are ordered to determine the position of each client. Once we know the couples (i, k) assigned for each open facility, the corresponding variables are added to the pool. According to R, the ranks of these costs are r 11 = 1 < r 33 = 3 < r 21 = 5. Thus, we get the couples (1, 1), (3, 2) and (2, 3). This means that client 1 goes to facility 1 in position 1, client 3 goes to facility 3 in position 2 and client 2 goes to facility 1 in position 3. Therefore, the variables y 1 {(1,1),(2,3)} and y 3 {(3,2)} are added to the pool.

Algorithm 3 Initial columns.

1: Input(|J| = p); 2: z = z(J); CreateSetVariables(J); 3: for n2 iterations,j1 ∈ J,j2 ∈ J do 4:

if z((J \ {j1}) ∪ {j2}) < z then 5: z = z((J \ {j1}) ∪ {j2}); J = (J \ {j1}) ∪ {j2}; CreateSetVariables(J); 6: end if 7: end for In order to test the usefulness of GRASP in solving problem instances, Table 1 reports results for the 60 instances of sizes n = 20, 30 enabling or not the use of the GRASP. It shows average results of CPU time (Time(s)), percent gap at termination, i.e. 100(z U B -z LB )/z LB (GAP(%)), and number of unsolved problems (in parentheses), number of nodes (#nodes) and number of variables (|V ars|). According to Table 1 it is clearly advisable to use the upper bound provided by the GRASP heuristic: it reduces the number of nodes, thus improving the size of the branch-and-bound tree.

GRASP

In Table 2, the same information as in Table 1 is reported but only for the instances solved to optimality within the time limit. One can observe from this table that enabling the use of GRASP reduces the CPU time and number of nodes of the B&B tree and at the same time reduces the overall number of variables required by the B&P&C. In addition, by using the GRASP heuristic, B&P&C is able to solve 6 more instances. For those instances for which B&P&C does not certify optimality, GRASP provides an upper bound that leads to an average gap of 0.89 %. Finally, without the use of GRASP, in many cases, no feasible solutions are found within the time limit and thus no % gap ("-") can be reported.

Our results show that by using the GRASP heuristic, 2.03% of the final number of variables are generated when applying Algorithm 3. The combination of the incumbent solution (given by GRASP) and that initial pool of variables leads to solve the considered instances faster, requiring less number of nodes and variables to certify optimality.

Figure 1 reports the performance profile of GAP versus number of solved instances within a time limit of 1800 seconds, for the 60 instances with sizes n = 20, 30. The blue line reports results using GRASP and the orange one without it. It is interesting to point out that when GRASP is enabled the B&P&C is able to solve to optimality 30 instances and the GAP of the remaining never goes beyond 6.43%. On the other hand, if GRASP is disabled then B&P&C solves only 24 instances. In addition, it is capable to obtain a feasible solution for only 2 more instances whereas in the remaining 34 instances the gap is greater than 100% (no feasible solution is found).

Stabilization

When using a column generation procedure, the vector of dual variables may be quite different from an iteration to the next resulting in a slow convergence. For this reason, the stabilization is sometimes a critical step in order to reduce the number of variables and iterations needed to solve each reduced master problem [START_REF] Du Merle | Stabilized column generation[END_REF]). We follow the stabilization procedure of [START_REF] Pessoa | Exact Algorithm over an Arctime-Indexed Formulation for Parallel Machine Scheduling Problems[END_REF] which depends on only one parameter. The idea consists in using a vector of dual variables which is a convex combination of the previous vector and the current solution of the dual problem.

Let π = (α, β, γ, δ, ) be a generic vector of dual multipliers, π be the best known vector of dual multipliers (found so far) and π ReM P be the current solution of the dual problem. Let c j S (π) be the reduced cost of y j S computed with the dual variable π and LB(π) the lower bound provided by the same vector of dual multipliers, namely π. Finally, let z D (π) be the value of the dual objective function of ReLRMP for the dual vector π, see ( 25). The stabilization algorithm that we have implemented is described by the following pseudocode:

Algorithm 4 Stabilization in ReLRM P .

1: ∆ = ∆init; π = 0; LB(π) = 0; GAP = 1; 2: while GAP > do 3: Solve ReLRMP, obtaining zReLRMP and πReLRMP ; πst = ∆πReLRMP + (1 -∆)π; 4:

for j = 1, . . . , n do 5:

Solve the pricing using πst, obtaining S; 6: if c j S (πReLRMP ) < 0 then Add variable y j S ; end if 7:

end for 8:

LB(πst) = z(π t st ) + GAP = z ReLRM P -LB(π)

z ReLRM P ; 17: if GAP < 1 -∆ then ∆ = 1 -GAP ; end if 18: end while
In words, the algorithm performs a while loop where in each iteration it makes a convex combination of the current vector of dual multipliers and the best vector of multipliers found so far. This loop ends whenever both vectors of multipliers are close enough based on the gap between the incumbent lower bound and the actual value of the reduced master problem. It is important to realize that the coefficient (importance), ∆, given in the convex combination to π ReLRM P (the current solution of ReLRMP) increases with the number of iterations of the algorithm since ∆ = 1 -GAP and GAP decreases with the number of iterations. Eventually in the very last iterations of the stabilization algorithm we will use the actual vector of dual multipliers since π st ≈ π ReLRM P .

In order to check the efficiency of the stabilization and to determine the best value for parameter ∆ init , LRMP has been solved to optimality for 270 instances from n = 20 to n = 100. In our implementation, we have chosen ∆ = 0.4 based on the computational study shown in Figure 2. As one can observe in this figure, the best performance profile is obtained by ∆ = 0.4 (blue dashed line) because it is the configuration that solves the instances in less time. It is worth mentioning that LRMP can be solved in one third of the time required for solving the problem without stabilization. Detailed results of the linear relaxation for ∆ = 0.4 are reported in Table 4. We report in Figure 3 the evolution of the lower and upper bounds with respect to number of iterations for a single instance. Stabilization results in a better behavior: the dual bound is not infinite at iteration 0 and it does not improve for some iterations. The reason is because we start with a feasible solution of the problem.

∆ init = 0.2 ∆ init = 0.4 ∆ init = 0.6 ∆ init = 0.8 No Stabilization
The control over the dual variables significantly improves the necessary number of iterations and the number of variables used to certify optimality. Note that this improvement becomes more important when M P is solved using a branch-and-bound procedure because the number of variables should be small at every node.

Preprocessing

In order to improve the performance of the algorithm we use two different preprocessings to set some variables to zero. Our approach is based on Claims 1 and 2 in [START_REF] Labbé | A comparative study of formulations and solution methods for the discrete ordered p-median problem[END_REF].

The reader may observe that although those results allow to fix some x k ij variables to zero, this variable-fixing can be translated to the new setting by using the relation (10) between the variables in formulations W OC and M P .

Therefore, the above results imply that those variables y j S such that (i, k) ∈ S and x k ij = 0 will not be considered to be added to the ReLRMP. This can be simply enforced by setting the corresponding d k ij = 0 in every pricing subproblem.

Branching strategies

Branching on the original variables is a common option when the master problem involves set partition constraints. See for instance [START_REF] Johnson | Modeling and Strong Linear Programs for Mixed Integer Programming[END_REF]. In spite of that, we have also considered other branching strategies as using the set partitioning variables or the Ryan and Foster branching, [START_REF] Ryan | An integer programming approach to scheduling[END_REF]; [START_REF] Barnhart | Branch-and-price: column generation for solving huge integer programs[END_REF]. However, these two alternatives were discarded because our pricing subproblem is polynomially solvable when we branch on the original variables whereas using any of the other branching strategies mentioned above, makes it N P-hard.

Recall that x k ij = S∈S(j):(i,k)∈S y j S , thus, a way to branch using a fractional solution can be derived directly from the integrality conditions on the original variables. Proposition 1. If x k ij ∈ {0, 1} for i, j ∈ I, k ∈ K, then y j S ∈ {0, 1} for j ∈ I, S ∈ S(j).

Proof. Suppose on the contrary there exists a variable with fractional value y j S . Since x k ij are binary for all i, j, k (in particular for i 1 , j , k 1 where (i 1 , k 1 ) is a pair of S ), there must be another fractional variable y j S such that (i 1 , k 1 ) ∈ S . Note that S = S since the column generation procedure never generates duplicate variables. Hence, there is a pair (i 2 , k 2 ) such that either (i 2 , k 2 ) ∈ S or (i 2 , k 2 ) ∈ S but not both. Therefore, we obtain the following relationship

1 ≥ S∈S(j ): (i 1 ,k 1 )∈S y j S > S∈S(j ): (i 2 ,k 2 )∈S y j S > 0. (35)
The first inequality comes directly from the formulation. The second inequality is strict because the term S∈S(j ):(i 2 ,k 2 )∈S y j S has at least one fractional variable less than the term S∈S(j ):(i 1 ,k 1 )∈S y j S . The third inequality is strict because of the choice of (i 2 , k 2 ). Finally, a contradiction is found because x j i 2 k 2 is not binary. The reader may note that this branching can be seen as a SOS1 branching [START_REF] Beale | Special facilities in a general mathematical programming system for non-convex problems using ordered sets of variables[END_REF]) since at most one of the above y j S variables can assume the value 1. The way to implement this branching in the pricing subproblem is to set locally (in the current node) to zero the y j S variables which are in conflict with the condition implied by the branch x k ij = 0 or x k ij = 1. In the case x k ij = 0 we set y j S = 0 for all sets S containing couples (i, k) ∈ S. Analogously, in the case x k ij = 1 we set y j S = 0 for all sets S containing (i, k) ∈ S such that j

= j , (i , k) ∈ S such that i = i or (i, k ) ∈ S such that k = k .
This condition can be transferred to the pricing subproblem modifying the d k ij coefficients accordingly. Specifically, this transformation is done as follows:

• If x k ij = 0 then d k ij = 0. • If x k ij = 1 then      d k ij = 0, j ∈ I : j = j. d k i j = 0, j , i ∈ I : i = i. d k ij = 0, j ∈ I, k ∈ K : k = k.
Moreover, it is also well-known that branching on SOS constraints (original variables) gives rise to more balanced branching trees (see e.g. Chapter 7 of [START_REF] Wolsey | Integer programming[END_REF]) than branching on the variables of M P .

Among the fractional original variables one has to decide which will be the next variable to branch on. One of the easiest techniques for this choice is to consider the most fractional variable. This is not difficult to implement but it is not better than choosing randomly [START_REF] Achterberg | Branching rules revisited[END_REF]). Alternative techniques are pseudocost branching [START_REF] Benichou | Experiments in mixed-integer programming[END_REF]) or strong branching [START_REF] Applegate | Finding cuts in TSP[END_REF]) although they are rather costly.

This issue has motivated us to propose another rule to select the variable to branch on, based on the improvement of the bounds in each of the new created nodes. We use the following indices corresponding to the down and up branches of the variable x k ij :

ς k,- ij = λ k c ij x k ij and ς k,+ ij = λ k c ij 1 -x k ij . (36)
They account, respectively, for the unitary contribution to the objective function due to fixing the variable x k ij either to zero (down branching) or to one (up branching). Branching down stimulates the improvement of the lower bound, whereas branching up helps the problem to find integer solutions.

We have tested several strategies that make use of the indices, ς, defined above.

Strategy 1: arg min{θς k,- ij + (1 -θ)ς k,+ ij : 0 < x k ij < 1}. Strategy 2: arg min{min{ς k,- ij , ς k,+ ij } : 0 < x k ij < 1}. Strategy 3: arg min{max{ς k,- ij , ς k,+ ij } : 0 < x k ij < 1}.
Based on our computational experience (see Figure 4), we have concluded that the best strategy to choose the following variable to branch on corresponds to strategy 1 with θ = 0.5.

Each node of the branching tree can be fathomed before it is fully processed comparing lower bounds, as given by ( 25) and ( 26), with the current incumbent solution. This strategy implies reducing the number of calls to the pricing subproblem and as a result savings in the number of variables added to the restricted master problem.

Valid inequalities

Clearly, the addition of valid inequalities (20) to M P modifies the structure of the master problem and thus the pricing must be modified accordingly. Let us denote by ζ k ij the dual variable associated with valid inequality (20) for indices i, j, k. After some calculation, one obtains the following expression of the reduced costs of variable y j S :

c j S = c j S + γ * j + δ * + i ∈I j ∈I        (i,k)∈S: r i j ≥rij k =1 ( * k + ζ k * i j ) + (i,k)∈S: r i j ≤rij k =n ( * k+1 + ζ (k+1) * i j )        - i∈I: (i,•)∈S α * i - k∈K: (•,k)∈S β * k . (37)
Furthermore, solving the pricing subproblem to find a new column or to certify optimality of the column generation algorithm requires to adapt the dynamic programming algorithm that computes the g(i l , k) terms using the new dual multipliers. This implies to modify the D j matrices. Once again, after some calculations the modified d k ij elements are now given by:

d k ij =                          λ k c ij + i ∈I j ∈I: r i j ≤r ij ( k+1 + ζ k+1 i j ) -α i -β k if k = 1, λ k c ij + i ∈I j ∈I: r i j ≥r ij ( k + ζ k i j ) + i ∈I j ∈I: r i j ≤r ij ( k+1 + ζ k+1 i j ) -α i -β k if k = 2, . . . , n -1, λ k c ij + i ∈I j ∈I: r i j ≥r ij ( k + ζ k i j ) -α i -β k , if k = n.
(38) These new elements allow us to apply the adapted column generation algorithm to solve LRMP, reinforced with valid inequalities (20).

To justify the use of the mentioned cuts we have done some preliminary computational experiments with instances of sizes n = 50 and 60. Table 3 compares the behavior of the standard branch-and-price without cuts, (B&P(M P )), against the strategy with cuts, B&P&C(M P ).

From Table 3, we conclude that it is always better to add cuts because the final gap is always smaller with this strategy. This solution scheme has been implemented and the results are reported in the next section.

Computational Experiments

The B&P&C implementation of the formulation M P has been experimentally compared with the B&C implementation of the formulation W OC on the instances detailed below.

The B&P&C algorithm considered in these experiments is based on the description in the previous section.

The computer used for these tests has an Intel Core i7 CPU clocked at 2.8 GHz with 4GB of RAM. Each implementation has a maximum of 7,200 seconds (2 hours) to solve each individual instance.

Both implementations are using the SCIP 4.0.1's API (see [START_REF] Gamrath | The SCIP optimization suite 3.2[END_REF]) and both are calling the LP solver of IBM ILOG CPLEX 12.7.

Instances

Since no standard libraries of instances for DOMP are available in public repositories we generate our own instances with the pseudorandom number generator from the C random library. In this work we consider that the sets of clients and potential facilities coincide, thereby we will refer both as points.

We consider 20 sets of 30 instances. Each set has a different number of points such that n ∈ {20, 30, . . . , 90, 100, 120, 140, . . . , 280, 300, 400}. For a given n, we generate one subset of 10 instances for each value of p, where p ∈ { (n/4) , (n/3) , (n/2) }.

For a given n, we first randomly generate the Cartesian coordinates of the points in the square [0, 400] 2 . Then, we calculate the cost for each pair of points with the Euclidean distance between the two related nodes in the square. We round each distance to the nearest integer to build the cost matrices. We also fix the values of the matrix diagonal to the smallest admissible cost to avoid free self service.

Finally, we randomly generate the weight vector λ such that λ k ∈ [n/4, n] for k ∈ K. All these instances, with n up to 400, are available at https://gom.ulb.ac.be/gom/ wp-content/uploads/2018/12/DOMP_Repository.zip. Detailed information about the instances generation can be found in [START_REF] Chvátal | Linear Programming (W. H. Freeman and Company)[END_REF].

M P vs W OC linear relaxations

We assess experimentally the linear relaxation of M P by comparing with W OC on all the instances generated. For these experiments, neither cuts nor preprocessing have been applied.

In Table 4, we report averages of the numerical results of the linear relaxation for both formulations. The value GapLP (%) represents the gap percentage between the optimal integer value z * (alternatively the best known solution) and the linear relaxation optimal value z * LP : GapLP (%) = 100(z * -z * LP )/z * LP . The computational times in the column T ime(s) are given in seconds. Table 4 also includes average number of variables (|V ars|) and required memory (M emory(M B)). We highlight the small number of variables that are generated to certify optimality with the column generation approach applied to M P , besides the time and memory saving, which is likewise significant. Table 4: Numerical results on linear relaxation for W OC and M P .

As expected, the integrality gap of formulation M P is smaller than the one of W OC. Moreover, formulation M P also outperforms W OC in the number of required variables (see Figure 5) which results in much smaller memory requirements (see Figure 6). Indeed, the implementation of W OC fails to solve, already for sizes of n = 100, the linear relaxation of all instances by lack of RAM memory; whereas with the same parameter configuration, formulation M P is relatively far from experiencing that problem. Figure 6 shows the performance profile of the memory requirement of both formulations. As one can see M P outperforms W OC with respect to this factor for all instance sizes.

B&P&C(MP) vs B&C (WOC)

We now compare the B&P&C implementation of M P with the B&C implementation of W OC. The former follows the procedure explained in Section 3 and the latter consists of W OC formulation with (9) as valid inequalities.

The results are reported in Table 5. In that table, we denote by T ime(s) the average computational time (in seconds) required by each method to obtain an optimal solution for a given set of 10 instances defined by number of clients (n) and number of open facilities (p). With |V ars| we refer to the average of the numbers of variables used by M P or W OC. We also denote by |N odes| and |Cuts| the average of the number of nodes explored and the average of the number of cuts used, respectively, in the corresponding methodology. The column #unsolved(T /M ) in the case of B&C(WOC) reports the number of unsolved instances out of the 10 in each group. It distinguishes between those instances not solved by exceeding the maximum running time (T ) or the memory limits (M ). Observe that in the similar column within the blocks B&P&C(MP) no distinction is shown since the memory limit is never reached and instances may be not solved only because of time limitation. Finally, we also include the gap at termination as Gap(%) = 100(z U B -z LB )/z LB , where z LB and z U B are the lower and upper bound, respectively.

Analyzing further the results in Table 5 we conclude that on average B&C(WOC) uses less variables than B&P&C(MP). This allows us to solve larger sized instances that were not affordable for the original W OC. We also observe that the number of required cuts for B&P&C(MP) is smaller than for B&C(WOC). This could be explained by the tightness the linear relaxation of DOMP problems of sizes n ≥ 100. This fact shows the usefulness of our new approach.

Conclusions

This paper presents a first branch-price-and-cut, B&P&C(MP), algorithm for solving DOMP. This approach is based on an extended formulation using an exponential number of variables coming from a set partitioning model. Elements in the partitions are couples containing information about a client and its sorted position in the sorted sequence of allocation costs. To address the solution of this formulation we develop a column generation algorithm and we prove that the pricing routine is polynomially solvable by a dynamic programming algorithm. We embed the column generation algorithm within a branch-and-price framework. Furthermore, we adapt preprocessing and incorporate families of valid inequalities that improve its performance. Extensive computational results compare the performance of our B&P&C(MP) against the most recent algorithm in the literature for DOMP, B&C(WOC), showing that for the largest considered instances B&P&C(MP) performs better and it requires less memory to upload and run the models. The methodology presented in this paper is able to solve sized instances for DOMP that had never been solved in the literature.

  1, where α = [0, 1], we obtain a convex combination of median and center objectives (centdian), etc. The p-facility Discrete Ordered Median Problem looks for the subset J of p facilities to open in order to minimize the ordered median function: min J⊆I:|J|=p z(J).

Example 4 .

 4 (Continuing from Example 1) We illustrate the use of the function CreateSet-Variables(J) with the following set J = {1, 3} (open facilities). The allocation costs for this set J of open facilities are c 11 = 1, c 21 = 3, c 33 = 1.

Figure 1 :

 1 Figure 1: Performance profile graph with GRASP enabled or disabled after 1800 seconds, GAP / # of instances.

Figure 2 :

 2 Figure 2: Performance profile graph with different combination of ∆ init , #solved instances / n.

Figure 3 :

 3 Figure 3: Bound's behavior at the root node in a particular instance on successive iterations.

Figure 4 :

 4 Figure 4: Performance profile graph of #solved instances using different branching strategies.

Figure 5 :

 5 Figure 5: Graph of Number of Variables versus size n for W OC and M P .

Figure 6 :

 6 Figure 6: Graph of Memory usage (MB versus size n for W OC and M P .

  S 3 , S 4 , S 5 , S 6 , S 7 , S 8 , S 9 , S 10 , S 11 , S 12 , S 13 , S 15 , S 17 , S 22 , S 23 , S 25 , S 28 }, S(2) = {S 1 , S 2 , S 3 , S 4 , S 5 , S 6 , S 7 , S 8 , S 9 , S 12 , S 13 , S 14 , S 17 , S 18 , S 19 , S 22 , S 23 , S 25 , S 30 }, S(3) = {S 1 , S 2 , S 3 , S 4 , S 5 , S 6 , S 7 , S 8 , S 9 , S 10 , S 11 , S 15 , S 16 , S 20 , S 21 , S 24 , S 26 , S 27 , S 31 }.

Table 3 :

 3 Numerical results with and without cuts.

Table 5 :

 5 |V ars| |N odes| |Cuts| #unsolved(T /M ) Gap(%) T ime(s) |V ars| |N odes| |Cuts| #unsolved Gap(%) Numerical results for B&C(WOC) and B&P&C(MP).

				B&C(WOC)					B&P&C(MP)
	n p T ime(s) 20 10 4.48	4211	38	689	0/0	0.00	58.70	4272	111	373	0.00
	30 15	131.89	13952	19197	2519	0/0	0.00 3670.21 16369	392	1078	0.45
	40 20 6202.85	32820	605812	4727	8/0	2.54 6751.27 14591	41	1819	3.33
	50 25 6575.59	63776	355560 10131	9/0	1.35 6860.93 13407	9	2709	1.87
	60 30 6707.38 109804	85723 15676	8/2	1.82 7200.20 16838	11	2864	4.47
	70 35 2474.98 173955	835 19238	2/8	7.72 7200.31 17758	2	4150	4.95
	80 40 3428.13 259186	1 12406	0/10	3.38 7201.36 18902	2	4680	2.13
	90 45 6243.49 368560	1 12157	7/3	4.24 7200.38 20028	7	4127	2.37
	20 6	11.50	5706	440	1249	0/0	0.00	952.76 10959	97	615	0.00
	30 10 1578.21	18245	305595	3056	1/0	0.12 6270.41 17502	81	1503	1.61
	40 13 7061.36	43664	628962	6559	8/2	2.38 7200.72 11186	2	3073	3.69
	50 16 7116.54	85630	284028 10423	9/1	1.14 7200.49 14977	6	3066	2.89
	60 20 3306.54 144983	20330 19887	2/8	3.02 7200.77 16691	1	4192	4.52
	70 23 2119.13 231680	1 23603	0/10	6.20 7200.97 19307	2	4365	4.51
	80 26 2886.25 346926	1 25187	0/10	5.59 7201.35 21675	1	5449	3.69
	90 30 5214.89 488316	1 32406	0/10	4.62 7201.66 24507	1	6116	3.17
	20 5	16.54	6054	1215	1537	0/0	0.00 1989.43 14699	97	731	0.09
	30 7	1807.41	20643	198424	4789	1/1	0.65 6840.10 14934	29	2093	1.85
	40 10 7050.93	48065	602685	7939	7/3	1.68 7200.90 10730	1	3455	3.95
	50 12 7200.00	94784	270959 12579	10/0	0.91 7201.39 14807	1	3526	2.83
	60 15 2768.88 161807	1 18081	0/8	2.90 7201.08 17131	2	3489	5.04
	70 17 1842.00 259406	1 16115	0/10	6.04 7201.78 19454	1	4649	4.51
	80 20 2902.00 383199	1 27129	0/10	6.95 7201.52 25278	3	4320	4.89
	90 22 5999.16 549561	1 46216	0/10	6.82 7201.75 27418	1	5549	5.50

For instances with n = 400 the time limit was set to

hours.

Acknowledgements

This research was supported by Spanish/FEDER grants number MTM2016-74983-C02-01. The research of the second and third authors was supported by the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy Office. We thank the SCIP team (Gamrath et al. (2016)) for the helpful advices.

23

of B&P&C(MP) with respect to B&C(WOC). After adding cuts B&P&C(MP) is able to solve the problem with a smaller branch-and-bound tree. The number of instances solved to optimality, for small size instances up to n = 40, is slightly better for B&C(WOC). As the size increases this number is similar in both cases. Gaps at termination, after 7200 seconds, are always smaller than 8% for B&C(WOC) and smaller than 6% for B&P&C(MP), being the later clearly better from instances of n = 70. Since B&C(WOC) is not able to handle any instance with n = 100 (reporting out of memory flags) we continue our study in Table 6: Numerical results for B&P&C(MP) for bigger instances.

Table 6 contains the results within the time limit of 2 hours for bigger instances of DOMP. This table has the same layout as Table 5 except that we replace column #unsolved by M emory(M B). This new column shows the average required memory to solve the corresponding set of instances. In that table extensive computational experiments are reported for instances up to 400 points. We would like to remark that the increase of the complexity, with respect to the instance sizes, of |V ars|, |Cuts|, M emory(M B) and Gap(%) is moderate (almost linear) which allows one to handle DOMP problems of larger size. Moreover the Gap(%) are similar to those reported in Table 5.

To conclude, the results show that the overall performance of B&P&C(MP) in solving DOMP is systematically better than the branch-and-cut formulation B&C(WOC) for instances for n ≥ 70. In addition, it is worth noting that B&C(WOC) is not even able to solve


Jules Waldhart

Mamoun Gharbi

Rachid Alami

A Novel Software Combining Task and Motion Planning for Human-Robot Interaction

Task planning and motion planning softwares operate on very different representations, making it hard to link them. We propose a software filling that gap in a generic way, mainly by choosing the best way to physically perform a task according to a higher level plan and taking explicitly into account human comfort and preferences.

The literature on task planning abounds in methods to synthesize to symbolic plan that fulfill an objective (e.g. HTN task planners [START_REF] Ghallab | Automated Planning: Theory and Practice[END_REF]). An other domain of the literature proposes methods to solve motion planning problems (such as RRT or PRM [START_REF] Lavalle | Planning Algorithms[END_REF]).

In other terms we know, thanks to task planning, which action to take, and thanks to motion planning we compute how to perform motions. The purpose of the software presented in this paper is to fill the gap between these two planners, providing a generic and easy to use interface enabling task planners to perform tests on action feasibility, and to make relevant geometric choices such as how to grasp an object, where to place it, or where to stop navigation in order to talk to a person.

These choices are made to assess the geometric feasibility of the plan, but the role played by our software, named Geometric Task Planning (GTP), is even more important when humans are involved in the task: one of its main purpose is to choose how to perform a task, and depending on human placement some ways to achieve a task might be more pertinent than others. Figure 1 depicts a set of options computed for a situation where the robot has to pick an object close to a human. According to the task requirements, the software chooses the best solution (configurations and trajectories) that minimizes the disturbance toward the human (figure 1 d), making the overall process human aware.

To the best of our knowledge, there is no other software dedicated to fill the gap in a generic way as GTP does. Most of the available software are specific to the task planner used [START_REF] Kaelbling | Integrated task and motion planning in belief space[END_REF][START_REF] Bidot | Geometric backtracking for combined task and motion planning in robotic systems[END_REF]. GTP design makes it easy to use, and it provides multiple interfaces using several middlewares.

The next section presents the geometric reasoning and planning software, including its algorithms and the features

Copyright c 2016, Association for the Advancement of Artificial Intelligence (www.aaai.org). All rights reserved. related to human-aware planning. We then briefly discuss an illustrative example and discuss future extensions.

Geometric Reasoning and Planning

The GTP software takes requests as inputs, in addition to an initial geometric world state (the positions of all objects, obstacles and agents (robots or humans) in the environment).

The inputs contains the type of action to perform, the involved agents and objects alongside with their role. The actions available for now are: pick, place, placeReachable, give, navigateTo. Here is an input example:

action: placeReachable, mainAgent:

robot, targetAgent:

Jane, mainObject:

box, supportObject: table2

In that example, GTP will compute a plan where the robot places the object named "box" on the "table2", in a place where it is reachable by a person called Jane. Here, the support object is optional, so GTP would find a suitable support to place the box in the reach of Jane. Other keywords are available for different usage, such as "target" for navigation.

The solution provided by GTP is a labelled sequence of trajectories performing the action specified. For example, the output for a pick action is a sequence of four trajectories: an "approach" trajectory where the robot gripper goes close to the object to pick, the second trajectory is a straight line to reach a grasping position, the "closing" of the end effector is a third trajectory, and finally an "escape" trajectory to lift the object from the surface it is lying on. The next subsections gives a sketch of the algorithms used in GTP.

Main Algorithm

The algorithm is based on four steps: Choices Making a choice concerning the action: for example in a pick, choosing the grasp to use, for a place, the final object position and so on. Configuration computation Finding the configuration to use, based on inverse kinematics and on the choices made in the previous step. Motion planning Computing all the trajectories based on the initial state and on the configurations found in the previous step. We use a cost-based motion planner with human-related costs [START_REF] Jaillet | Sampling-Based Path Planning on Configuration-Space Costmaps[END_REF][START_REF] Mainprice | Planning human-aware motions using a sampling-based costmap planner[END_REF]). Cost computation Computing the cost of the action, e.g. reflecting its safety, legibility, human comfort. The first iterations are done without step 3 (Motion Planning), which samples alternative sets of key configurations to achieve the task. These alternatives are sorted in term of cost. We then compute the motion of the best alternative and iterate over the following if no valid path is found. We finally compute the full task cost, including the path. When a certain number of choices are explored and no solution is found, the action is set as infeasible.

Details

Costs used The costs are mainly related to human preferences and comfort. They are inspired by [START_REF] Mainprice | Planning human-aware motions using a sampling-based costmap planner[END_REF][START_REF] Sisbot | A human-aware manipulation planner[END_REF]. The main effect is to keep the robot as far as possible from the human. Virtually any cost can be used.

Relations and Affordances GTP is able to compute a set of relations between objects such as is on or is in, but also relations between an agent (human or robot) and an objects (affordances) such as is reachable by or is visible by. These relations can be transformed into symbolic facts that can be used at a higher level to determine the effects of actions.

Assumptions We plan in a fixed world, and assume the human is not moving in the plan. We are working on a way to adapt on-line the motion to human movement, in a taskdependent way, but this is not presented here.

Usage and interface

GTP has already been successfully used for different purposes, on the simulation and on the real robot1 [START_REF] Gharbi | Combining Symbolic and Geometric Planning to synthesize human-aware plans: toward more efficient combined search[END_REF][START_REF] Devin | An implemented theory of mind to improve human-robot shared plans execution[END_REF]. It has two interfaces: the first one is a GUI available through the software move3d [START_REF] Simeon | Move3d: A generic platform for path planning[END_REF] and the second interface is a ROS interface where the same functionalities are made accessible.

Conclusion

This software provides to task planners an easy access to a number of functionalities that a service robot must have such as basic manipulation and navigation, and can integrate other, like handovers [START_REF] Mainprice | Sharing effort in planning human-robot handover tasks[END_REF][START_REF] Waldhart | Planning handovers involving humans and robots in constrained environment[END_REF]. GTP can be used for human-robot interaction, looking for optimal solutions in a given costspace. It was designed to interact closely with the Hierarchical agent-based planner (HATP) as depicted in [START_REF] Gharbi | Combining Symbolic and Geometric Planning to synthesize human-aware plans: toward more efficient combined search[END_REF]. As a result, GTP will soon be made available to use as stand alone and can be linked through ROS in order to execute the trajectories.

Figure 1 :

 1 Figure 1: Options computed for a pick action, with initial situation (a) and 3 grasping configuration from the worst (b) to the best (d), according to human-related cost. Solution c is the best when forcing the use of the right arm.

Videos are available at http://frama.link/2V4ZCOW7 and https://youtu.be/2iBElhKi-AE

Acknowledgements

This work is supported by the European Union's H2020 programme under grant agreement No. 688147, MuMMER http://mummer-project.eu/


Mary Ellen Foster
email: maryellen.foster@glasgow.ac.uk

Rachid Alami
email: rachid.alami@laas.fr

Olli Gestranius
email: olli.gestranius@ideapark.fi

Oliver Lemon
email: o.lemon@hw.ac.uk

Marketta Niemelä
email: marketta.niemela@vtt.fi

Jean-Marc Odobez
email: odobez@idiap.ch

Amit Kumar Pandey
email: akpandey@aldebaran.com

The MuMMER Project: Engaging Human-Robot Interaction in Real-World Public Spaces

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Developing an artificial agent capable of coexisting and interacting independently, naturally, and safely with humans in an unconstrained real-world setting has been the dream of robot developers since the very earliest days. In popular culture and science fiction, the prototypical image of a "robot" is precisely this: an artificial human that is able to engage fully in all aspects of face-to-face conversation. As modern robot hardware becomes safer, more sophisticated, and more generally available, there is a clear and significant consumer demand for robots that are able to coexist in everyday, real-world human environments in this way. For example, demand for SoftBank Robotics' Pepper robot has been so high in Japan that, since July 2015, each monthly run of 1000 units has sold out in under a minute [17]; also, the Jibo robot broke records on the crowdfunding site Indiegogo in July 2014 when it raised US$1million in less than a week [START_REF] Annear | Makers of the world's 'first family robot' just set a new crowdfunding record[END_REF]. This form of socially intelligent HRI [START_REF] Dautenhahn | Socially intelligent robots: dimensions of human-robot interaction[END_REF] has therefore received a particular focus in recent years.

However, while the hardware capabilities of such robots are increasing rapidly, the software development has not kept pace: even with the most recent technological developments, the most advanced such robots have generally supported limited, scripted interactions, often relying on a human operator to help with input processing and/or appropriate behaviour selection.

In the new European project MuMMER ("MultiModal Mall Entertainment Robot"), we are developing an entertaining humanoid robot that is able to operate autonomously and naturally in the different dynamic environments of a public shopping mall. The overall concept underlying MuMMER is that an interactive robot deployed in a public space such as a shopping mall should be entertaining: this will increase the acceptability of such a robot to the people it interacts with, therefore improving their overall experience in the mall. To ensure that the robot scenarios and interactions are engaging and entertaining, we will involve a wide number of potential users through continuous situated co-design studies, and will use the results of these studies to inform the design and operation of all components of the robot system.

In a number of recent studies, humanoid robots have been deployed in a range of public spaces, mainly in Japan. For example, Kanda et al. [START_REF] Kanda | A communication robot in a shopping mall[END_REF] carried out a five-week field trial of a robot in a shopping centre, involving 2343 interactions with 235 tracked participants (as well as many other members of the general public who did not choose to take part in the tracking study). The robot was semi-autonomous: it used a human operator to carry out speech recognition, to monitor and override its behaviour selection, and to provide additional domain knowledge when needed. In post-experiment questionnaires, the robot was highly rated on most subjective measures, with many positive comments. Outside Japan, other successful locations for public robot deployment have included museums [START_REF] Karkaletsis | Natural interaction with personality and dialogue enabled robots[END_REF], city centres [START_REF] Bauer | The autonomous city explorer: Towards natural human-robot interaction in urban environments[END_REF][START_REF] Weiss | Robots asking for directions: The willingness of passers-by to support robots[END_REF], care homes [START_REF] Hebesberger | Lessons learned from the deployment of a long-term autonomous robot as companion in physical therapy for older adults with dementia: A mixed methods study[END_REF], and airports [START_REF] Triebel | Spencer: A socially aware service robot for passenger guidance and help in busy airports[END_REF] .

Note, however, that the interactive behaviour of the robots in these previous deployments has been limited in various ways. For example, most of the robots used in these previous deployments either were not mobile, were remotely teleoperated, or were able to move only within an area delimited with environmental sensors such as floor sensors. Also, these robots are generally semi-autonomous, incorporating a human operator for operations such as speech recognition and overall behaviour selection. In general, recent experiments with robots in public spaces have tended to concentrate either purely on human-aware navigation (e.g. [START_REF] Clodic | Rackham: An interactive robot-guide[END_REF][START_REF] Kruse | Human-aware robot navigation: A survey[END_REF]), or purely on natural-language interaction (e.g. [START_REF] Barker | Towards human-robot speech communication in everyday environments[END_REF]), rather than on the combination of the two.

Objectives

In MuMMER, we will develop a humanoid robot capable of performing its interactive tasks independently and naturally in a dynamic environment and with several users simultaneously. To meet this goal, we have defined the following set of objectives which will be used to guide the work across the project. These objectives encompass both state-of-the-art technical requirements, as well as our intentions to ensure that our work in MuMMER is grounded in real-world needs and applications throughout the project.

Objective 1: Development of an interaction robot for entertainment applications.

The tasks of interaction and user experience design for entertainment are at the core of MuMMER. These highly challenging objectives are multiparty and multimodal at numerous levels (sensing, synthesis, conversation). We aim to achieve these objectives through the seamless integration of user state and mood monitoring, design of entertainment-oriented artificial social signals, and interaction management. Objective 2: Fusion of interactive applications through a co-design process. Interactive applications tend to be specific to one or two simple tasks. However, in MuMMER, we aim to develop an interaction robot that blends different tasks (information providing, guidance, entertainment) within a single entity, relying on real user needs and interests. The goal is to ensure that the developed robotic application will be user-driven in terms of HRI, is socially and ethically accepted, and is also interesting to commercial end-users. Objective 3: Situated-perception with on-board sensors. Rather than instrumenting the scene or people, in MuMMER we aim to design an autonomous robot able to operate robustly and naturally in diverse environments by relying only on onboard sensors. To address the resulting sensing limitations (e.g., reduced field of view and visibility, noisy acoustics), we plan to design a truly interaction-aware perception approach combining perception algorithms aware of the robot gestures, and active sensing, along with speech synthesis and dialogue to make users implicitly or explicitly aware of perception uncertainties. Objective 4: Automatic learning of engaging interaction strategies. We aim to develop robots which can interact naturally with humans in a social setting, being entertaining as well as helpful, and this requires learning from human interaction data. We therefore aim to develop new techniques for automatic learning of interaction strategies for socially appropriate, engaging, and robust robot action selection, based on uncertain information about the 'social' state of the interaction (e.g., where the user(s) are, whether users are engaged in the interaction, what speech has been recognised, what information is in the common-ground, and the estimated perspectives of the different users). Such capabilities will create a step-change in the naturalness and acceptance of human interaction with robots. Objective 5: Human-aware and situated interactive tasks and motion synthesis. We aim to investigate further, develop, and demonstrate a set of decisional capabilities that will allow such a robot to effectively "close the loop" between situated dialogue and robot behaviour. The robot will be able to estimate the visual perspective of its human interaction partners, their positions and postures, and based on the current interaction and dialogue context, to synthesize the pertinent human-aware motion and interaction. This will be achieved based on a novel combination of perspective-taking, situation assessment, and human-aware task and motion planning in order to make a substantial step forward and produce more complex and validated dialogue and interaction-aware behaviours. Objective 6: Development of new business opportunities. One of the most challenging outcomes of the MuMMER project will be the development of novel business

Approach

In the preceding section, we set out the objectives of the MuMMER project in areas including technical development as well as the broader context in which this work is situated. In this section, we describe our concrete approach to meeting these objectives: beginning with a description of the selected deployment locations and discussing how scenarios will be co-designed with a range of stakeholders; next, summarising how the necessary technical advances will be achieved and integrated on a sophisticated humanoid robot; and finally, outlining our plans for real-world deployment and user acceptance studies of the integrated robot system throughout the project.

Stakeholders, scenarios, and co-design

Throughout the duration of the project, we will carry out continuous data collection, deployment, and evaluation activities in the Ideapark shopping centre in Lempäälä, Finland. Ideapark is the biggest commercial city in Scandinavia, with 7.3 million visitors in 2013. It contains almost 200 stores, restaurants and cafes within 100 000 square metres. We are targeting two main classes of interaction scenario in MuMMER. In the first set of scenarios, the robot will be situated in the main shopping area of the mall (Figure 1a): in addition to supporting a set of helping tasks, such as providing guidance, information, or collecting customer feedback about the events going on, the robot will also engage fully with users, proposing events to attend or performing entertainment activities such as telling jokes and riddles, organizing contests, or distributing items such as balloons. In the other scenarios, the robot will be situated in the Pii Poo centre8 (pictured in Figure 1b), which is designed as an accessible space providing inclusive cultural activities for all ages and abilities. In this case, the robot will provide entertainment activities to children or families, focussed on social games such as quizzes, memory games, and the like.

The details of the scenarios will be developed through a continuous co-design process, which will ensure that multiple relevant perspectives from consumer target groups as well as from commercial and other stakeholders will be identified, analysed and integrated into the design of the robot, the human-robot interaction, and the applications for future businesses. In particular, the co-design process will work along with technical developments to defined a set of relevant scenarios of increasing complexity (e.g. among the entertaining possibilities), increasing the chances that the robot will perform according to the goals set, while being representative of future business applications identified with the stakeholders.

Co-design will include workshops, seminars and scenario co-creation from both consumer and business perspectives. In addition, co-design will support more focused study methods to answer specific questions or tests, such as questionnaires, interviews, and observations. Scenarios and requirements will be co-created with the stakeholders to guide the technical development work. Also the success metrics for the development and assessment of human-interactive mobile robots for consumer markets will be co-created to ensure their validity with regard to the purpose (the validity will be tested in the field trials). The work on success metrics will be supported by market analyses and market studies. Co-design engages people in the technical development, and MuMMER will work to integrate their contributions to the technology development and reach mutual learning of different stakeholders throughout the development. Co-design increases acceptance and interest of consumers toward such robots; this is a major research question and goal in the project.

Technical development

From a technical point of view, endowing a robot with the ability to support the desired interaction scenarios outlined above will require the integration of a number of state-of-the-art components on a sophisticated robot platform. The primary target robot platform for MuMMER is SoftBank Robotics' Pepper9 robot (Figure 2). Pepper is a 1.20m tall, wheeled humanoid robot, with 17 joints for graceful movements, 3 omni-directional wheels to move around smoothly, and more than 12 hours of energy for non-stop activities, along with the ability to go back to the recharging station if required. It has a 3D camera to sense and recognize humans and their movements at a distance of up to 3 meters. The software is built on the same architecture as the existing Nao robot, and can be programmed in a very efficient and open way. For the initial stages of the MuMMER project, which focus on scenario development, data collection, and preliminary development, we will use the default Pepper hardware which is currently available commercially. However, based on the particular needs of MuMMER, an adapted Pepper prototype will be developed during the first period of the project, including additional sensors if necessary and feasible; this will be the robot that is used for the final autonomous MuMMER system in the later stages. The high-level architecture of the MuMMER system is presented in Figure 3. As shown in this diagram, the technical developments will fall into four main classes: audiovisual processing, social signal processing, interaction management, and interactive task and motion control. In the remainder of this section, we summarise how we will develop state-of-the-art components to address each of these tasks.

Audiovisual processing

To improve the performance of the perception modules, we will investigate several research directions. First, we will extend on our previous real-time visual-only probabilistic multiparty head and pose tracking framework used in [START_REF] Sheikhi | Leveraging the robot dialog state for visual focus of attention recognition[END_REF]. In particular, to handle non-stationary sensors and situations with multiple people occluding each other and not facing the camera, we will leverage depth data to gain performance, e.g. by combining RGB-D based upper-body detectors (in addition to body part classifiers [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF]) with very fast segmentation based tracking methods [START_REF] Bibby | Robust real-time visual tracking using pixel-wise posteriors[END_REF][START_REF] Bibby | Real-time tracking of multiple occluding objects using level sets[END_REF] on the depth data stream. Second, we will develop tighter links between perception and action, and design situated algorithms that have a better knowledge of the situation and are able to exploit it. The situational context can originate from different sources (interaction state, perception state and quality, action state, etc.) and can be used in different ways: by influencing recognition and decisions, actively triggering specific perception algorithms or behaviuors (to quickly verify with a head gesture that a person that is currently outside of the field of view is still present), or via model parameter adaptation in the longer term.

Social signal processing Social signal processing involves both the recognition of human social signals, as well as the generation and synthesis of social signals by the robot; in MuMMER, we will advance the state-of-the-art in both of these areas. For social signal recognition, we will extend algorithms for detecting and classifying user engagement and attention to apply also to the target domain of multi-party, real-world, entertaining human-robot interaction. For social signal generation, we will investigate whether "non-natural" social signals such as lights and acoustic effects can actually be socially useful in this context. We will also extend the investigation of humour from the verbal (i.e., what makes a text funny?) to the behavioural cues that can make others laugh (i.e., why two people telling the same joke are not necessarily equally funny?) [START_REF] Hurley | Inside Jokes[END_REF]. Finally, we will will investigate how synthetic personality traits can influence the relationship between a robot and a possibly large audience (rather than a single interaction partner).

Onboard sensors

Interaction management

Multimodal Fusion for Social State Estimation

Synthesis of Social Signals

Social signal processing Fig. 3: Technical architecture of the MuMMER system; arrows show the main data flow Interaction management In this project, we will advance the state-of-the-art in interaction management in two directions. First, we will apply current state-of-the-art statistical models (see e.g. [START_REF] Keizer | Machine learning of social multi-party human-robot interaction[END_REF][START_REF]Data-driven Methods for Adaptive Spoken Dialogue Systems: Computational Learning for Conversational Interfaces[END_REF]) to the new scenarios and tasks that arise in the context of engaging, entertaining, socially appropriate human-robot dialogue interaction. Second, we will scale up machine learning techniques to support robust human interaction with a robot in noisy, populated spaces. To support the MuMMER scenarios and tasks, the range of system actions will be extended to include communicative gestures and behaviours such as moving towards a particular person, as well as pro-active behaviour, such as initiating an engagement with a person to attract their attention. In order to enable the robot to make informed decisions what to do next, given this extended range of possible system behaviours, the social state model will need to be extended with the additional relevant features. As a consequence, the extended state and action model will pose additional challenges of scalability of the machine learning techniques for state tracking and action selection that must be addressed.

Interactive tasks and motion control

Most of the prior work in the area of human-aware motion planning are focused mainly from a planning perspective, to be able to exhibit safe and comfortable motions. Further, their target scenario is not the kind of crowded space we are aiming for, with potentially many people, standing, moving, interacting and many a times blocking the path. In this context, another important dimension becomes more prominent: considering the social and cultural favours and interaction the robot and the human can exchange in the environment. This aspect has not been sufficiently exploited in the context of navigation planning. For example the robot can reason from the human's perspective, analyze the feasibility that the human can provide the way for the robot and requests the human to do so through verbal or physical interaction. Similarly, the human should be able to request for the same and expect the similar natural response from the robot. Such reasoning will greatly affect the way the state of the art navigation planning and decision making system works. This requires a fusion of rich reasoning about humans, their abilities, affordances, perspective taking, interaction and navigation planning. With this, we aim to elevate the human-aware and socially-aware path planning towards more intelligent planning with social interaction possibilities. Hence, the aim is to achieve interactive-aware social motion and behaviour synthesis capabilities in a real life environment.

Real-world deployment and acceptance studies

A crucial aspect of the MuMMER project plan is that the robot will be regularly deployed in Ideapark throughout the whole course of the project. For the first 3.5 years of the project, the interim deployments will be regular (i.e., every few months) but rather short-term (maximum 1-2 weeks); these visits will serve to support the co-design activities as well as to provide formative evaluation of the technical components of the system and pilot testing of the overall robot deployment process.

However, at the end of the project, we will carry out a continuous, long-term user acceptance study in Ideapark lasting several months, where the robot is present in the shopping centre continuously (ideally, daily) throughout the entire period: this should allow a wide range of behaviours to be examined, and also to greatly reduce the novelty effect of the robot by providing more concrete data on how it is actually used in context across time. The results of this long-term study will serve to finalise the use cases and success metrics developed during the project, and will also provide a concrete demonstration of the utility, flexibility, and generalisability of the robot system developed in MuMMER.

As mentioned at the start of this paper, this sort of long-term, public-space deployment has only recently begun to be addressed (e.g., [START_REF] Hebesberger | Lessons learned from the deployment of a long-term autonomous robot as companion in physical therapy for older adults with dementia: A mixed methods study[END_REF][START_REF] Triebel | Spencer: A socially aware service robot for passenger guidance and help in busy airports[END_REF]), and not very widely outside of Japan. The results of the MuMMER long-term deployment should therefore add significantly to our knowledge of how such robots are received and accepted in public spaces around Europe.

Current work

The project is currently in its first year, and efforts have begun towards all of the technical tasks outlines in Section 3. In particular, the first co-design workshops have been carried out with management and personnel from Ideapark, which has resulted in initial scenario proposals (including the idea of taking a 'selfie' with the robot), and the first consumer workshop will be carried out shortly. Consultations with retailers will begin with interviews in the autumn, and the first Pepper deployment in the mall will also take place in the autumn. We have also begun to lay the groundwork for the final acceptance study by developing a questionnaire which will be administered regularly throughout the project to monitor the attitudes of Ideapark customers towards the robot: the first questionnaire will be carried out before Pepper is taken to the mall in the autumn to provide a baseline measure of these attitudes.

On the technical side, development has also begun on all of the components outlined above in Section 3, and the infrastructure to allow the components to work together on Pepper is also being prepared. All partners have also begun development and testing with the actual Pepper hardware platform. We anticipate being in a position to test an initial integrated system in the mall at the end of the project's first year, as a step towards the eventual long-term deployment of the full system at the end of the MuMMER project's fourth year.

  Fig. 1: Deployment locations

Fig. 2 :

 2 Fig. 2: Pepper robot (Images © SoftBank)

http://www.kulttuuripiipoo.fi/content/fi/1/20106/In+English.html

https://www.ald.softbankrobotics.com/en/cool-robots/pepper

Acknowledgements

This research has been partially funded by the European Union's Horizon 2020 research and innovation programme under grant agreement no. 688147 (MuMMER, mummer-project.eu).