1 Introduction

At times, very hard combinatorial optimization problems contain easy combinatorial subproblems after relaxing some of their constraints. A paradigmatic example is the Traveling Salesman Problem: after the elimination of its subtour elimination constraints it turns into the Linear Assignment Problem which is polynomially solvable. This pattern calls for developing techniques that take advantage of this situation to solve some combinatorial problems based on their constraint relaxation. This approach is not new and the reader is referred to Focacci et al. (1999, 2002a, 2002b) and the references therein for further details.

This behavior is not only observed in problems where formulations include an exponential number of constraints. Actually, it also occurs in many polynomial size formulations. One of these cases is the Discrete Ordered Median Problem (DOMP) modeled with the strong order constraints formulation as introduced in Labbé et al. (2017). If one removes the family of strong order constraints, whose acronym is SOC, the resulting problem is the standard p-median problem that is known to be combinatorially friendly (Hakimi 1964; Marín and Pelegrín 2019; ReVelle and Swain 1970). The aim of this work is to develop solution techniques for DOMP based on constraint relaxations.

The DOMP is a discrete location model that allows to generalize several classical discrete location problems (see, e.g., Nickel and Puerto 2005). For instance, the discrete p-center and p-median are particular cases of DOMP. Assume that we are given a set of clients, a set of candidate locations for facilities, and the allocation costs from each candidate facility to each client. The objective of DOMP is to locate p facilities in such a way that a certain weighted function of the allocation costs is minimized. These weights are not assigned to specific costs but to their sorted values. Namely, the weighted average sorts the allocation costs in a nondecreasing order and then, it performs the scalar product of this so obtained sorted cost vector by the vector of weights.

In the literature, one can find different applications of the ordered median operator. For instance, it has been applied to facility location (Aouad and Segev 2019; Domínguez and Marín 2020; Espejo et al. 2009; Kalcsics et al. 2010; Martínez-Merino et al. 2017; Tamir 2001), multicast communication (Fourour and Lebbah 2020), multiobjective Markov decision processes (Ogryczak et al. 2011), voting problems (Ponce et al. 2018), supervised classification (Marín et al. 2022), tomography reconstruction (Calvino et al. 2022), and network design (Puerto et al. 2013), among other situations.

DOMP was first introduced in Nickel (2001) as an integer nonlinear problem. Then, in Boland et al. (2006), this problem was modeled as a mixed integer linear program. Some works on this problem take advantage of some particular characteristics. Specifically, Marín et al. (2009) introduce an efficient covering formulation for DOMP considering free self-service, ties in the cost matrix and a non-negative weighted order vectors in the objective function. Futhermore, Marín et al. (2010) present a covering reformulation for weighted order vectors containing zeros and an extended model for vectors even with negative elements.

In Labbé et al. (2017), a new three-index formulation based on set packing constraints is proposed. These set packing constraints are known as strong order constraints or SOC, and the number of these constraints appearing in the formulation is \(\mathcal {O}(n^3)\). In addition, another new formulation, solved by an efficient branch-and-cut procedure that provides good results in terms of time, is introduced.

This second formulation is based on the aggregation of the SOC corresponding to the same position. The resulting order constraints are the weak order constraints (from now on WOC). This formulation includes SOC as valid inequalities. Both formulations present small integrality gaps.

Recently, in Deleplanque et al. (2020), a novel branch-price-and-cut algorithm has been proposed. This procedure is based on a formulation with an exponential number of variables that corresponds to a set partitioning model. The proposed approach allows to handle larger instances since it requires less memory to run the model.

In this paper, we want to explore different exact approaches to solve DOMP. The first one uses branch-and-cut techniques based on one of the most promising formulations, namely the formulation based on WOC proposed in Labbé et al. (2017), adding SOC as valid inequalities. Additionally, we compare the use of cut pools in the branch-and-cut with respect to the use of callbacks to implement an ad hoc separator proposed in Labbé et al. (2017). By setting up pools of cuts, all SOC are initially stored and then, solvers decide which cuts are included during the branch-and-cut process. In contrast, applying the callbacks using the separation algorithm introduced in Labbé et al. (2017), SOC are not initially stored and the implementation of the SOC separation is based on a sequential update of the left-hand side of the corresponding order constraints. This separation can be performed in \(\mathcal {O}(n^3)\).

The second method is based on a constraint relaxation on the formulation using SOC to model the order. This procedure, to solve DOMP, starts with a relaxed formulation where all SOC are removed and feasibility is enforced adding model constraints from the SOC family in the searching tree. Although the number of SOC is polynomial, \(\mathcal {O}(n^3)\), this number of constraints becomes too large to be handled when n increases. Consequently, it is interesting to study this approach since we could improve the time and memory performance by only including the necessary constraints in the solution process. Again, in this case, we compare the branch-and-cut procedure through callbacks with respect to the branch-and-cut based on constraint pools.

The contributions of this paper can be summarized as follows:

  1. (1)

    Comparing a branch-and-cut approach to solve DOMP based on the so called WOC formulation with a constraint relaxation approach for DOMP based on removing SOC.

  2. (2)

    Comparing the performance of the branch-and-cut and the constraint relaxation approach when using callbacks based on specific tailor made separation oracles with respect to the use of fixed pools of cuts/constraints.

  3. (3)

    Reporting intensive computational tests which show the limits of the different considered solution methods.

The remainder of this work is organized as follows. In Sect. 2, we introduce the notation and description of DOMP. Besides, we recall two formulations for DOMP that will be used along the paper (\(\text {DOMP}_{\text {WOC}}\) and \(\text {DOMP}_{\text {SOC}}\)). In Sect. 3, we present two solution methods for the problem. First, we describe the branch-and-cut procedure for \(\text {DOMP}_{\text {WOC}}\) introduced in Labbé et al. (2017). Then, we propose a novel row generation procedure for \(\text {DOMP}_{\text {SOC}}\). Section 4 is devoted to the analysis of both solution methods. In addition, we present a comparison between the results of using pools of cuts/contraints and using callbacks for the branch-and-cut and the row generation techniques. Finally, in Sect. 5, we include some conclusions and future research lines.

2 Problem definition and formulations

This section is devoted to recall the definition and some formulations of DOMP that will be instrumental in our discussion. We shall follow the following notation. We denote by \(I=\{1,\ldots ,n\}\) the set of n clients and, at the same time and without loss of generality, the set of n potential facility locations. Facilities are assumed to be uncapacitated, i.e., they can supply as many clients as desired. Besides, \(c_{ij}\) denotes the cost for serving client i from facility j, for \(i,j\in I\).

Given a set J composed by p open facilities, \(c_i(J)\) represents the cost of allocating client i to the facility set J, i.e., \(c_i(J)=\min _{j\in J}c_{ij}\). In addition, if the vector of costs \(c_i(J)_{i=1,\dots ,n}\) for \(J\subset I\) is sorted in non-decreasing order, we denote by \(c^{(k)}(J)\) the allocation cost in position \(k\in K=I\) of this sorted vector. Thus, it holds that \(c^{(1)}(J)\le c^{(2)}(J)\le \ldots \le c^{(n)}(J)\).

The aim of DOMP is to determine a subset of p facilities \(J\subset I\) to open, and to assign each client to an open facility in order to minimize the ordered median objective function. Given a vector \(\lambda =(\lambda ^k)_{k\in K}\) such that \(\lambda ^k\ge 0\), \(k\in K\), the objective function of DOMP can be expressed as \(\sum _{k\in K}\lambda ^k c^{(k)}(J).\) Consequently, the definition of DOMP is

$$\begin{aligned} \min _{J\subset I:|J|=p}\sum _{k\in K}\lambda ^k c^{(k)}(J) \quad \quad \quad \text {(DOMP)}. \end{aligned}$$

Observe that this formulation generalizes several standard discrete location problems. For instance: if \(\lambda ^1=\lambda ^2=\ldots =\lambda ^n=1\), this model is the p-median problem; if \(\lambda ^1=\lambda ^2=\lambda ^{n-1}=0,\lambda ^n=1\), one gets the p-center; if \(\lambda ^1=\lambda ^2=\ldots =\lambda ^{n-k}=0,\lambda ^{n-k+1}=\ldots =\lambda ^n=1\), the resulting problem is the k-centrum; etc.

DOMP is known to be \(\mathcal{NP}\mathcal{}\)-complete, see Nickel and Puerto (2005), and as mentioned in the introduction, different formulations have been proposed to deal with this problem. In this paper, we will elaborate on two of the most recent and promising formulations presented in Labbé et al. (2017). On the one hand, that paper introduces a three-index formulation where order is modeled by a family of set packing constraints, SOC. On the other hand, it also presents an aggregated version of that formulation where the order is ensured by a different set of constraints called WOC. We will build on those two formulations, thus in order to be self-contained, in the following subsection we provide full details of them.

2.1 Strong order constraints formulation

For the formulation based on strong order constraints, the next families of variables are required:

$$\begin{aligned} y_{j}= & {} {\left\{ \begin{array}{ll} 1,&{}\text{ if } \text{ facility } j \text{ is } \text{ open, }\\ 0,&{}\text{ otherwise, }\\ \end{array}\right. } \text{ for } j\in I,\\ x_{ij}^k= & {} {\left\{ \begin{array}{ll} 1,&{}\text{ if } \text{ client } i \text{ is } \text{ allocated } \text{ to } \text{ facility } j \text{ and } \text{ the }\\ &{}\text{ associated } \text{ cost } \text{ is } \text{ in } \text{ position } k \text{ of } \text{ the }\\ &{} \text{ sorted } \text{ sequence } \text{ of } \text{ allocation } \text{ costs, }\\ 0,&{}\text{ otherwise },\\ \end{array}\right. }\text{ for } i,j\in I, k\in K. \end{aligned}$$

Besides, we denote the rank of the allocation cost \(c_{ij}\) by \(r_{ij}\), i.e., \(r_{ij}=\ell\) if \(c_{ij}\) is the \(\ell\)-th element in the list of the costs \(c_{ij}\), for all \(i,j \in I\), sorted in a non decreasing sequence and where ties are broken arbitrarily. Then, \(\text {DOMP}_{\text {SOC}}\) formulation is the following.

$$\begin{aligned} {(\text {DOMP}_{\text {SOC}})} \min&\displaystyle \sum _{i\in I}\sum _{j\in I}\sum _{k\in K}\lambda ^kc_{ij}x_{ij}^k \end{aligned}$$
(1)
$$\begin{aligned} \text{ s.t. }&\displaystyle \sum _{j\in I}\sum _{k\in K}x_{ij}^k&=1,&i \in I, \end{aligned}$$
(2)
$$\begin{aligned}&\displaystyle \sum _{i\in I}\sum _{j\in I}x_{ij}^k&= 1,&k\in K, \end{aligned}$$
(3)
$$\begin{aligned}&\displaystyle \sum _{k\in K}x_{ij}^k&\le y_j,&i,j\in I, \end{aligned}$$
(4)
$$\begin{aligned}&\displaystyle \sum _{j \in I}y_{j}&= p, \end{aligned}$$
(5)
$$\begin{aligned}&\sum _{i'\in I}\sum _{\begin{array}{c} j'\in I:\\ r_{i'j'}\le r_{ij} \end{array}}x_{i'j'}^k + \sum _{i'\in I}\sum _{\begin{array}{c} j'\in I:\\ r_{i'j'}\ge r_{ij} \end{array}}x_{i'j'}^{k-1}&\le 1,&i,j \in I,\; k\in K, k \ne 1,\quad \quad \quad \quad {\text {(SOC)}} \end{aligned}$$
$$\begin{aligned}&x_{ij}^k,\; y_j&\in \{0,1\},&i,j \in I,k\in K. \end{aligned}$$
(6)

Constraints (2) ensure that each client is served by just one facility in one position. Similarly, constraints (3) are necessary to guarantee that only one allocation cost is in each sorted position. Constraints (4) ensure that each client is allocated to an open facility and that the allocation cost of a client can only be in at most one position. Constraint (5) restricts that exactly p facilities must be open. Constraints (SOC) are the so-called strong order constraints which ensure the correct sorting of the allocation costs. These constraints are set packing constraints, i.e., at most one of the variables of the l. h. s. could take value one. The incompatibility of two or more variables taking value one is due to (4) and the fact that a variable \(x_{ij}^k\) cannot take value one if \(x_{i'j'}^{k-1}=1\) when \(r_{ij}<r_{i'j'}\), for \(i,j,i',j' \in I, k\in K, k \ne 1\). We refer the reader to Labbé et al. (2017) for a more detailed explanation. Finally, constraints (6) are the domain of definition of the variables. The reader may note that removing (SOC) the formulation results in the p-median.

2.2 Weak order constraints formulation

Despite the good mathematical properties of \(\text {DOMP}_{\text {SOC}}\), as the lower bound given by its linear relaxation, the number of (SOC) is \(\mathcal {O}(n^3)\). Consequently, when the number of clients increases, the number of (SOC) becomes too large to be handled by a solver. For this reason, Labbé et al. (2017) introduce an alternative family of constraints to ensure the order of costs. These new constraints are based on the aggregation of (SOC) corresponding to the same position. (The reader is referred to Labbé et al. (2017) for further details.) This alternative formulation results in the following.

$$\begin{aligned} {(\text {DOMP}_{\text {WOC}})} \min&\displaystyle \sum _{i\in I}\sum _{j\in I}\sum _{k\in K}\lambda ^kc_{ij}x_{ij}^k\nonumber \\ \text{ s.t. }&(2) - (6),\nonumber \\&\sum _{i\in I}\sum _{j\in I}\left( \sum _{i'\in I}\sum _{\begin{array}{c} j'\in I:\\ r_{i'j'}\le r_{ij} \end{array}}x_{i'j'}^k + \sum _{i'\in I}\sum _{\begin{array}{c} j'\in I:\\ r_{i'j'}\ge r_{ij} \end{array}}x_{i'j'}^{k-1} \right)&\le n^2,\nonumber \\&k\in K, k \ne 1.\quad \quad \quad {\text {(WOC)}} \end{aligned}$$

Constraints labeled by (WOC) are known as weak order constraints. They ensure that if facility j serves client i and its cost \(c_{ij}\) is in position k of the sorted cost vector of the solution, then there must be a smaller or equal allocation cost in position \(k-1\). This is due to the coefficients corresponding to each variable in the constraints. In each inequality, there are represented two positions (\(k-1\) and k). By constraints (3), only two variables must take value 1, and the remaining ones take value 0. Assuming that the variables with value one for position k and \(k-1\) correspond to positions s and t of the sorted costs, respectively, the inequality can be expressed as follows:

$$\begin{aligned} (n^2-(s-1))x_{i_s j_s}^k+tx_{i_t j_t}^{k-1}\le n^2, \end{aligned}$$

with \(i_s,j_s,i_t,j_t\in I\) such that \(c_{i_sj_s}\) and \(c_{i_tj_t}\) are the s-th and t-th smallest allocation costs in matrix \((c_{ij})_{n\times n}\). This is valid if and only if \(t<s\).

In Labbé et al. (2017), it is shown that \(\text {DOMP}_{\text {SOC}}\) formulation provides a relevant improvement of the integrality gap with respect to \(\text {DOMP}_{\text {WOC}}\) formulation. In other words, for most of the instances, fractional solutions satisfying (WOC) could be cut including (SOC). Thus, it is recommended to use (SOC) as valid inequalities of \(\text {DOMP}_{\text {WOC}}\).

3 Solution methods

Both formulations, \(\text {DOMP}_{\text {SOC}}\) and \(\text {DOMP}_{\text {WOC}}\), can be solved by using standard MIP solvers as CPLEX, Gurobi, Xpress, or SoPlex. However, the good performance of these formulations is rather limited for large sizes of the problem as we will see in Sect. 4. The reader should observe that already for \(n=100\) clients the number of (SOC) is almost \(10^6\).

To improve the performance of \(\text {DOMP}_{\text {WOC}}\), Labbé et al. (2017) propose a branch-and-cut procedure which starts by solving the linear program relaxation of \(\text {DOMP}_{\text {WOC}}\), and then it includes (SOC) as valid inequalities whenever necessary. This procedure is described in Sect. 3.1.

In this paper, we propose an alternative solution method which exploits the good properties of \(\text {DOMP}_{\text {SOC}}\) avoiding the use of the complete family of strong order constraints. This solution method consists in a row generation procedure which initially considers \(\text {DOMP}_{\text {SOC}}\) without (SOC), and then it iteratively includes these order constraints. As far as we know, this row generation method has not been considered before for DOMP. Section 3.2 is devoted to the description of this procedure.

3.1 Branch-and-cut for \(\text {DOMP}_{\text {WOC}}\)

The branch-and-cut procedure has become an efficient method for solving large instances of models where the number of constraints is intractable for solvers. For instance, it has been successfully applied to the matching problem with blossoms (Edmonds and Johnson 1973; Grötschel and Holland 1985; Letchford et al. 2004), problems related to trees (Fernández et al. 2017; Magnanti and Wolsey 1995), clustering (Benati et al. 2017), and the orienteering arc routing problem (Archetti et al. 2016), to name a few. In Labbé et al. (2017), a branch-and-cut method for (SOC) in \(\text {DOMP}_{\text {WOC}}\) formulation is proposed.

This branch-and-cut procedure could be handled by two different perspectives. On the one hand, most of the current solvers have some options to define a fixed pool of cuts that are added automatically when necessary in the cut generation. This means that, changing some parameters in the solver, we can remove some families of constraints from the formulations (thus avoiding to include them initially), which are later added when necessary in the solution process. This automatic feature of solvers is interesting when the constraint separation must be done by enumeration due to the efficient implementation of the solvers in this case.

In our particular case, the use of cut pools in the branch-and-cut method for (SOC) in formulation \(\text {DOMP}_{\text {WOC}}\) seems to require \(\mathcal {O}(n^6)\) operations, since there are \(\mathcal {O}(n^3)\) (SOC) and \(\mathcal {O}(n^3)\) x-variables in the formulation. Actually, it is \(\mathcal {O}(n^5)\) since each constraint has only \(\mathcal {O}(n^2)\) variables to check. Nevertheless, based on the knowledge of the structure of the problem, an efficient separation method of (SOC) constraints can be developed. This ad hoc separation can be included by callbacks allowing a more efficient implementation of the branch-and-cut.

Focusing in DOMP, Labbé et al. (2017) propose an algorithm to separate (SOC) with complexity \(\mathcal {O}(n^3)\). It is a remarkable quadratic improvement with respect to the pure enumerative approach. Algorithm 1 shows in detail this separation procedure, based on the calculation of left-hand sides of the possible cuts adding and subtracting two values in each iteration.

In Sect. 4, we will develop a complete analysis of the branch-and-cut method using pools of cuts and the branch-and-cut method using callbacks. We will determine whether or not it is advisable to use a separation algorithm that takes advantage of the knowledge of the problem through callbacks.

Remark 1

According to previous experiences (Labbé et al. 2017; Deleplanque et al. 2020), when valid inequalities (SOC) are embedded in a branch-and-cut procedure over \(\text {DOMP}_{\text {WOC}}\) formulation, they should be added at the root node, but not deeper, in order to find a compromise between the integrality gap and the size of the problem. Hence, for our computational study, we will use this cut-and-branch procedure to check the performance of the solution method proposed in this section.

3.2 Row generation procedure for \(\text {DOMP}_{\text {SOC}}\)

Since \(\text {DOMP}_{\text {WOC}}\) formulation presents coefficients that are not zero-one, in this work we explore the use of a (SOC) relaxation of \(\text {DOMP}_{\text {SOC}}\) adding iteratively these constraints whenever they are necessary. Therefore, the solution method of this section starts with a formulation which has a zero-one coefficient matrix to later add set packing constraints. Hence, we provide a well-behaved (from the solvers point of view) formulation of DOMP without using a huge number of constraints.

The initial formulation which is considered in this row generation procedure is the following.

$$\begin{aligned} {(\text {DOMP}_{\text {relax}})} \min \quad&\displaystyle \sum _{i\in I}\sum _{j\in I}\sum _{k\in K}\lambda ^kc_{ij}x_{ij}^k\nonumber \\ \text{ s.t. }\quad&(2) - (6). \end{aligned}$$

Observe that this formulation corresponds to the DOMP model without imposing the order constraints. Therefore, we are dealing with a relaxation of \(\text {DOMP}_{\text {SOC}}\). In this case, the proposed relaxation results in the p-median problem.

For each obtained solution in the branch-and-bound of \(\text {DOMP}_{\text {relax}}\), (SOC) are checked by using Algorithm 1 and added to the model when necessary. Consequently, this row generation method ensures the order by only using a moderate number of (SOC). This allows to handle bigger instances to be solved in a reasonable computing time as we will see in Sect. 4.

figure a

3.2.1 Bounds in constraint relaxations

In constraint relaxations, any integer solution has to be checked to be valid according to the problem definition. However, there are different alternatives for continuing the branch-and-bound tree exploration when a fractional solution arises. One of them is checking all model constraints to improve the lower bound. Another one is to branch in a particular fractional variable.

One issue to be taken into account is the way in which a subset of (SOC), that a solution does not verify, is selected to be included in the formulation in order to improve the lower bound without increasing too much the formulation size. To develop this idea, we introduce the following proposition.

Proposition 1

Given an integer solution \((\widehat{x},\widehat{y})\) for \(\text {DOMP}_{\text {relax}}\), if \(\widehat{x}\) verifies that

$$\begin{aligned} \sum _{i'\in I}\sum _{\begin{array}{c} j'\in I:\\ r_{i'j'}\le r_{ij} \end{array}}x_{i'j'}^k + \sum _{i'\in I}\sum _{\begin{array}{c} j'\in I:\\ r_{i'j'}\ge r_{ij} \end{array}}x_{i'j'}^{k-1} \le b,\quad i,j \in I,\; k\in K, k \ne 1, \end{aligned}$$
(7)

for \(b\in [1,2)\), then \(\widehat{x}\) satisfies all (SOC).

Proof

If \(b=1\), then (7) are equal to (SOC). Thus, the result follows trivially. Assume that \(b>1\). In this case, since \(\widehat{x}\) is integer and \(b<2\), then the left hand side of (7) must be at most one. Consequently, \(\widehat{x}\) satisfies (SOC). \(\square\)

As a result of Proposition 1, when an integer solution is obtained in the branch-and-bound tree of \(\text {DOMP}_{\text {relax}}\), the lhs described in Algorithm 1 could be compared with \(b\in [1,2)\) instead of comparing it with 1. Therefore, if \(lhs>b\) in Algorithm 1, then the corresponding (SOC) are included.

However, varying this b value could affect the number of added cuts when a fractional solution is found and thus, the number of explored nodes in the branch-and bound tree. When b is close to 2, then the number of identified constraints (7) which are not verified by the solution is smaller and they are the most violated cuts. Consequently, for big values of b, the number of added cuts will be reduced. Nonetheless, since the number of added cuts is smaller, the number of explored nodes in the branch-and-bound is expected to be bigger.

Remark 2

Following Remark 1, we separate fractional solutions only at root node. Hence, for deeper nodes, Algorithm 1 is called only when integer solutions are found obtaining upper bounds. Beyond these concerns, lower and upper bounds get closer within the branch-and-bound tree as usually.

In order to experimentally check how the value of b could impact in times, cuts, and nodes of the row generation proposed in this section, we present Table 1. We show the results for the instances of sizes \(n=20\), 30, and 40 that will be detailed in Sect. 4. Particularly, in Table 1, first column shows the number of clients; the second column reports the number of open facilities; the third set of columns shows the computing time for each b value; the fourth set of columns represents the number of added cuts in the row generation procedure; and finally, the last group of columns reports the number of explored nodes in the branch-and-bound tree. In all cases, each row reports the average value of ten instances. We have tested the results for \(b=1\), \(b=1.1\), and \(b=1.3\).

In Table 1, we can observe that \(b=1\) reports better results than \(b=1.1\) and \(b=1.3\) when \(n=40\). Note that for the largest instances, although the number of added cuts in the branch-and-bound process is bigger, the number of nodes and the computing times are smaller since the gaps at the root node are smaller. Consequently, for the computational results reported in Sect. 4, the choice of \(b=1\) is used in the row generation procedure.

Table 1 Time, cuts, and nodes for different r. h. s. in the row generation procedure using Algorithm 1 to separate (7)
Full size table

4 Computational experiments

This section is devoted to the analysis of the solution methods introduced in this paper. The goals of this computational study are the following: (1) checking the differences among the results of formulations \(\text {DOMP}_{\text {SOC}}\) and \(\text {DOMP}_{\text {WOC}}\) and determining their limitations; (2) comparing two approaches to implement the branch-and-cut and the row generation algorithms for DOMP. The difference between these two approaches relies on the fact that the first one uses a fixed constraint/cut pool handled by the solver and the second one applies the separator described in Algorithm 1 implementing a callback; (3) comparing the different results between the two solution methods defined in Sect. 3.

The instances used in this computational study were introduced for the first time in Deleplanque et al. (2020). These instances were created to test different weighted order vectors \(\lambda\) beyond the classical ones, namely p-median, p-center, k-centrum problems, etc. The weighted vector \(\lambda\) was randomly generated such that \(\lambda ^k \in \left[ \frac{n}{4}, n\right]\) for \(k \in K\). Furthermore, in that data set, there are small- to large-sized instances to perform an exhaustive computational study. The reader can find the mentioned instances in (Deleplanque et al. 2022).

The models were coded in C and solved with SCIP v.6.0.2 (Achterberg 2009; Gleixner et al. 2018) using as optimization solver CPLEX 20.1.0 on a Mac OS Catalina with a Core Intel Xeon W clocked at 3.2 GHz and 96 GB of RAM memory.

In the computational experience, for all the considered formulations and solutions methods, we have included a preprocessing phase. In this stage, we are able to reduce the number of necessary variables to define the problem in terms of optimality. We refer the reader to Labbé et al. (2017) for more details. Besides, we have given an incumbent solution provided by a GRASP heuristic (Deleplanque et al. 2020). This solution let us provide a good upper bound from the beginning of the corresponding solution method.

Table 2 contains the results within two hours of 90 instances up to 40 clients, namely ten instances of each configuration of n and p for two different formulations: \(\text {DOMP}_{\text {WOC}}\) and \(\text {DOMP}_{\text {SOC}}\). This table and the following ones show the average results for these ten instances: the average CPU time (Time), the number of instances not solved in the time limit (#Unsolved), the gap at the root node (GAProot(%)), the gap at termination (GAP(%)), the number of variables after preprocessing (Vars), the number of constraints (OrigCons), the number of cuts/constraints added in the procedure (Cuts), the number of nodes (Nodes), and the required memory (Memory (MB)). Observe that, for instances with \(n=30\) and \(n=40\), \(\text {DOMP}_{\text {SOC}}\) formulation provides better results than \(\text {DOMP}_{\text {WOC}}\). \(\text {DOMP}_{\text {SOC}}\) presents a better linear relaxation value (see gap at the root node), and it needs less nodes at the branch-and-bound tree. Consequently, \(\text {DOMP}_{\text {SOC}}\) requires less solution time. In addition, we can conclude that the large number of constraints implies an increment of memory for \(\text {DOMP}_{\text {SOC}}\) which makes this formulation too heavy for sizes of n greater than 40. Then, together with the fact that \(\text {DOMP}_{\text {WOC}}\) cannot solve any of the instances for \(n=40\), we study other alternatives which add (SOC) iteratively until certifying optimality. Thus, in the following tables, we report the results of the methods proposed in Sect. 3.

Table 2 Results of formulations (\(\text {DOMP}_{\text {WOC}}\)) and (\(\text {DOMP}_{\text {SOC}}\)) for \(n\in \{ 20,30,40\}\)
Full size table

Nowadays, commercial solvers include options to code valid inequalities in a branch-and-cut procedure. In this context, these valid inequalites are usually known as user cuts, while in constraint relaxations, these model constraints are known as lazy constraints. The coding of the cuts/constraints can be done just giving them as input of the linear program. These automatic strategies need to encode all the cuts/constraints in advance within a fixed pool, with an ensuing waste of computer memory. Also, the management of these potential cuts/constraints follows a pre-implemented strategy without taking advantage of the particularities of the formulation beyond the solver pattern recognition based on developers’ experience. On the other hand, the use of an oracle (in our case, Algorithm 1) to add the cuts/constraints on the fly implementing callbacks could save memory on the whole process (see, e.g., Ackooij and de Oliveira 2014; Blado and Toriello 2021; de Oliveira and Sagastizábal 2014; Mazzi et al. 2021; Wolf et al. 2014, and the references therein). Besides, it allows to control when to check and to add those cuts/constraints that is an advantage by itself. For instance, in the row generation solution method, we check model constraints (SOC) at the root node (regardless the solution is fractional or integer) and in any node with integer solution. We refer the reader to (CPLEX 2022; SCIP 2022a, b) for a detailed discussion.

Table 3 B&C: \(\text {DOMP}_{\text {WOC}}\) with (SOC) as valid inequalities
Full size table
Table 4 Row generation, i.e, SOC (model constraints) added iteratively
Full size table

Table 3 presents the results, in two hours of time limit, of the branch-and-cut introduced in Sect. 3.1, i.e., (\(\text {DOMP}_{\text {WOC}}\)) with (SOC) as valid inequalities. Here, we follow two different strategies: we code (SOC) defining a fixed user cut pool in SCIP and the solver decides when to check and to add them (Pool); or we check the constraints using Algorithm 1 which adds them by an user callback when needed at the root node (Callback). The results exhibit that the method using Algorithm 1 shows better solving times than the automatic approach. Besides, the method using the callback can solve nine more instances than the automatic one. Note that these better results are explained by the efficiency of our separation algorithm. Regarding the required memory, observe that memory used by the automatic method increases quickly with n. Therefore, the application of this method does not seem to be useful for bigger instances since they could not be loaded: it requires around 50 GB of RAM memory already for \(n=60\).

Table 4 reports the results, in two hours of computing time, of the row generation method, i.e., \(\text {DOMP}_{\text {relax}}\) with (SOC) as model constraints not included from the beginning. Two approaches to carry out the row generation are considered: the automatic use of (SOC) defining a fixed lazy constraint pool (Pool) and the application of Algorithm 1 to add (SOC) when necessary (Callback). In this table, we have included the same columns as in Table 3. Note that the callback approach provides the best computing time results and only 16 instances remain unsolved after the time limit. The automatic method requires less cuts and nodes than Callback. However, the required memory increases faster when using the automatic approach because it needs to encode all the original (SOC) constraints which are \(\mathcal {O}(n^3)\). Consequently, the performance of the row generation following the automatic separation shows, in general, worse results than using Algorithm 1.

From Tables 3 and 4, we can conclude that the performance of the automatic branch-and-cut and the automatic row generation are quite limited in comparison with the use of the separation presented in Algorithm 1. The reason of the better performance of the callback approach is that Algorithm 1 exploits the knowledge of the problem. Therefore, for the study of larger instances, we focus on the branch-and-cut and row generation approaches using Algorithm 1.

Table 5 shows a comparison between the results provided by the callback versions of the branch-and-cut method (B&C) and the row generation technique (RowGen). For these experiments, we establish a time limit of five hours. Overall, the row generation technique outperforms the branch-and-cut in terms of computational times. Moreover, the row generation is able to solve more instances than the branch-and-cut in five hours.

Table 5 B&C vs row generation in 5 hours of time limit
Full size table
Table 6 B&P &C, B &C, and row generation in 24 hours of time limit
Full size table

Note that, for some huge instances, the problem does not get even the root node bounds. For those instances, \(\texttt {GAProot(\%)}\) and \(\texttt {GAP(\%)}\) report the same value and thus, \(\texttt {GAProot(\%)}\) cannot be analyzed as the gap at the root node, but as the gap at the root node at the time limit.

Regarding \(n=60\) instances, the branch-and-cut is able to solve 20 out of 30 instances in five hours, whereas this solution method certifies optimality for only 13 instances in two hours (see Table 3). These three extra hours also let the row generation algorithm to solve 24 instances, seven more than the same algorithm in two hours (see Table 4).

For most of the instances, the integrality gap at termination provided by the row generation procedure is smaller than the one obtained by the branch-and-cut, even with less cuts added. This gives us the idea that the added cuts are more accurate when (WOC) family is not included in the formulation. However, for \(n=100\), the gap at termination is smaller for the branch-and-cut procedure since the added cuts cannot improve the lower bound given by the linear relaxation of the program in both solution methods within the time limit.

To analyze the differences between the branch-and-cut and the row generation, starting from \(\text {DOMP}_{\text {WOC}}\) and \(\text {DOMP}_{\text {relax}}\), respectively, we give the solver up to 24 hours of time limit. Thereby, the influence of the linear relaxation bound is not so decisive. Furthermore, the branch-and-price-and-cut (B&P &C) described in Deleplanque et al. (2020) has been tested for those instances in the same computer and with the same time limit. The reader could see, in Table 6, how the gaps at termination are smaller on average for the row generation procedure. In fact, they are reduced by half and for the row generation are less than 2%. This approach needs less cuts to have a reasonable bound at the root node what let it branch faster to generate more nodes and improve the bounds. Thus, whereas the solution method detailed in Sect. 3.1 solves four instances out of 60 for these large-sized instances, the solution method proposed in Sect. 3.2 is able to solve eight instances to optimality and on top of that, the row generation procedure also reduces the gap at termination for those instances which are not solved. The B&P &C procedure uses less variables but the gap at termination is worse because it is not able to solve even the root node. However, this approach would still be useful for bigger instances since it requires much less memory.

5 Conclusions

In this work, we have introduced a row generation solution method which has improved the best known performances for DOMP regarding the medium-sized instances of the data set described in Deleplanque et al. (2020). In adittion, we have improved the best known solution for an instance with \(n=90\) and \(p=45\) (domp90p45v5.domp) that can be found in the mentioned dataset.

For large-sized instances, the lower bound provided by formulations which include (WOC) makes them also a good alternative. For these instances, the lower bound given by the linear relaxation can be barely improved within the time limit. Moreover, comparing with the integrality gap given by the branch-and-price algorithm (Deleplanque et al. 2020), one should note that the column generation of its master problem gives theoretically better lower bounds.

Taking into account that the limits of our row generation algorithm come from the huge number of variables, a combination of row and column generation seems to be a promising approach to be considered as future research line.