# Models and Methods for Two-Level Uncapacitated Facility Location problems

Bernard Gendron, Paul-Virak Khuong, Frédéric Semet

► To cite this version:

Bernard Gendron, Paul-Virak Khuong, Frédéric Semet. Models and Methods for Two-Level Uncapacitated Facility Location problems. Springer. Contributions to Combinatorial Optimization and Applications, 358, Springer Nature Switzerland, pp.59-75, 2024, International Series in Operations Research & Management Science, 10.1007/978-3-031-57603-4_4. hal-04384267

HAL Id: hal-04384267

https://inria.hal.science/hal-04384267v1

Submitted on 10 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Models and Methods for Two-Level Uncapacitated Facility Location problems

## Bernard Gendron

CIRRELT (Interuniversity Research Centre on
Enterprise Networks, Logistics and Transportation)
Département d'informatique
et de recherche opérationnelle
Université de Montréal
C.P. 6128, succ. Centre-ville
Montreal, Quebec H3C 3J7

## Paul-Virak Khuong

CIRRELT (Interuniversity Research Centre on
Enterprise Networks, Logistics and Transportation)
Département d'informatique
et de recherche opérationnelle
Université de Montréal
C.P. 6128, succ. Centre-ville
Montreal, Quebec H3C 3J7
pvk@pvk.ca

## Frédéric Semet

Univ. Lille, CNRS, Inria, Centrale Lille, UMR 9189 CRISTAL
F-59000 Lille, France
frederic.semet@centralelille.fr

June 2023
Abstract

We consider the two-level uncapacitated facility location problems with single-assignment
constraints (TUFLP-S), a problem that arises in industrial applications in freight trans-
portation and telecommunications. First, we recall the real-life application that moti-
vates this research work. Then, we present an integer programming formulation for this
variant with modular costs of the TUFLP-S, and we describe an adaptation of vari-
able neighborhood search to the classical TUFLP-S and the variant with modular costs.
Focusing on the TUFLP-S, we introduce mathematical programming models based on
reformulation techniques and on the relaxation of the integrality of some of the variables
associated with location decisions. We then present a Lagrangian relaxation approach
for the TUFLP-S, based on solving a single-level uncapacitated facility location problem
(UFLP) as the Lagrangian subproblem. We also describe a matheuristic based on a
mixed-integer programming-based large neighborhood search. Synthetic computational
results on instances derived from the industrial application as well as on large, hard,
artificial instances are reported throughout the paper.
# 1 Introduction

This article aims to provide a synthesis of the research work carried out by B. Gendron and his co-authors on the two-level uncapacitated facility location problem (TUFLP). The motivation to work on this line of research came from a project with a multi-channel retailer in the early 2000s [Gendron et al., 2002; Gendron and Semet, 2009]. The TUFLP is an extension of the uncapacitated facility location problem (UFLP) [Krarup and Pruzan, 1983]. Solving an UFLP consists of selecting a set of depots from potential locations to minimize the sum of the fixed costs associated with each depot and of the transportation costs from any depot to each customer. In the TUFLP, the single set of locations is substituted with two tiers of locations (depots and satellites), and the path to each customer must start at a depot and pass through a satellite. In other words, the only arcs are between depots and satellites and between satellites and customers. A classical variant of the TUFLP is the two-level uncapacitated facility location problem with single assignment constraints (TUFLP-S). The TUFLP-S imposes the additional restriction that each satellite can be connected to at most one depot. These constraints appear in a number of applications, including transportation [Gendron and Semet, 2009], and telecommunications [Chardaire et al., 1999]. For a large class of TUFLP instances for which the unique assignment constraints are not explicitly enforced, there exists an optimal solution that satisfies these constraints, due to the structure of the objective function [Chardaire et al., 1999].

In this article, we address two versions of the TUFLP-S: the *classical* TUFLP-S and the TUFLP-S with modular costs. In the classical TUFLP-S, the objective function includes fixed costs associated with the depots and the satellites, fixed costs for establishing connections between depots and satellites, and transportation costs from any depot to each customer, i.e., each path of the form depot-satellite-customer has a corresponding transportation cost. In the TUFLP-S with modular costs suggested to us by the industrial application, we have modular costs at depots and satellites associated with the number of product units handled and modular costs between depots and satellites
associated with the number of vehicles used.

Multi-echelon location analysis is a central research theme in supply chain manage-
ment reviewed by Melo et al. [2009] and Ortiz-Astorquiza et al. [2018]. Most real-life
applications lead to more challenging problems than the TUFLP, since they include
multi-period and multi-commodity aspects or a complex cost function to optimize, see,
e.g. Amrani et al. [2011], [Gendron and Semet, 2009], Huang et al. [2010], Thanh et al.
[2010].

If we restrict ourselves to TUFLPs, there is a wide body of literature as illustrated by
the surveys of Klose and Drexl [2005] and by Şahin and Süral [2007]. Two main formula-
tions are considered: arc-based [Kratica et al., 2014; Marín, 2007; Pirkul and Jayaraman,
1996, 1998] and path-based models [Barros and Labbé, 1994; Gao and Robinson Jr, 1992;
Kaufman et al., 1977; Ro and Tcha, 1984]. Most authors propose exact solution meth-
ods based on strengthening the models with valid inequalities and facets [Aardal et al.,
1996; Baïou and Barahona, 2014; Chardaire et al., 1999; Landete and Marín, 2009], on
Benders decomposition [de Oliveira et al., 2020] and on Lagrangian relaxation (Barros
[1998]; Marín and Pelegrín [1999]; Pirkul and Jayaraman [1996, 1998]. Some researchers
compare relaxations for the arc-based and path-based formulations [Bloemhof-Ruwaard
et al., 1994, 1996; Chardaire et al., 1999; Marín and Pelegrín, 1999]. Considering the
approximate solution of TUFLPs, heuristic methods based on Lagrangian or linear pro-
gramming relaxations have been developed Barros et al. [1998]; Pirkul and Jayaraman
[1996, 1998], while other heuristic approaches rely on greedy strategies or on simple neigh-
borhood search methods based on add, drop or exchange moves [Barros and Labbé, 1994;
Mitropoulos et al., 2009; Narula and Ogbu, 1979; Zhang, 2006]. To the best of our knowl-
edge, there are very few attempts to solve TUFLPs with metaheuristics. Chardaire et al.
[1999] describe a simulated annealing procedure that produces good solutions, while Ig-
nacio et al. [2008] present a tabu search algorithm that is shown to be efficient compared
to solving the mathematical program using CPLEX.

The remainder of this article is organized as follows. First, in Section 2, we recall the
real-life problem at the origin of this research work introduced in [Gendron et al., 2002],
and give a path-based integer programming model which was proposed in [Gendron and Semet, 2009]. Next, the two heuristic approaches to solve this problem described in Gendron et al. [2015] are presented in Section 3. Since this problem is a variant of the TUFLP-S, several formulations for this problem proposed in [Gendron et al., 2017] are presented and compared in Section 4, followed by the description of solution methods developed (see Gendron et al. [2016]) in Section 5.

## 2 A real-life TUFLP-S with modular costs

The motivation to study TUFLP originated in a project for a company that sells a wide range of products (clothing, electronics, and appliances) through its website, mail-order catalogs, and physical stores [Gendron et al., 2002; Gendron and Semet, 2009]. One of the main logistical challenges it faces is to adapt its distribution system to demands that vary daily to meet a level of service (typically, a maximum delivery time of 24 hours) while minimizing the total cost. As a B2C company, consolidation is a major concern. Indeed, most orders contain few items, resulting in small to medium packages.

The company addresses this challenge through a multi-echelon distribution network with single-sourcing (see Figure 1). Its operations are based in a few fixed primary facilities, which are central warehouses owned by the company. All transportation activities are subcontracted to carriers. A fleet of large vehicles delivers parcels from these primary facilities to secondary facilities, which are cross-docking terminals owned by carriers. Medium-size vehicles depart from these terminals and transfer parcels to smaller cross-docking locations (tertiary facilities), where the items are sorted into pre-assigned tours for final delivery to customers.

The distribution network design problem is to determine the secondary and tertiary facilities to open and the number of trucks to connect the primary and secondary facilities and the secondary and tertiary facilities. Indeed, the company owns the primary facilities, and their locations are known. Moreover, it is assumed that the routes to serve the customers from a satellite are pre-computed.

Since there are no fixed costs associated with primary facilities and all products are
available at each of them, each secondary facility can be assigned to its closed primary
facility without losing optimality. Thus, the problem can be modeled as a TUFLP-S with
modular costs. The transportation costs for individual connections are step functions
that take into account the distance traveled and the capacity of vehicles. A fixed cost
is associated with the use of each secondary facility. Secondary facilities incur costs for
each batch of parcels sorted, which are also represented by step functions. To keep the

<figure><img src="image_2.png" /><figcaption>Figure 1: Initial distribution network</figcaption></figure>

vocabulary consistent with the following sections, secondary facilities will be referred to
as depots and tertiary facilities as satellites. Gendron and Semet [2009] propose an arc-
based and a path-based formulation for this real-life variant of the TUFLP-S. We present
here only the path-based model, which outperforms the arc-based formulation. Indeed,
even though the model size is larger, the linear relaxations are stronger and the better
upper bounds are better when the path-based model is solved with CPLEX for a given
CPU time. Let us first define the different sets of nodes:

$I$ : set of potential sites to locate depots;

$J$ : set of potential sites to locate satellites;

$K$ : set of customers;

$I_j$ : set of potential sites to locate depots connected to satellite $j \in J$;

$J_i$ : set of potential sites to locate satellites connected to depot site $i \in I$;
$K_j$: set of customers connected to satellite $j \in J$;

The location and the transportation costs associated with the facilities and the customers
are denoted as follows:

$f_i$: fixed cost for using and operating depot $i \in I$;

$g_j$: cost per $Q$ product units for using and operating satellite $j \in J$;

$d_i$: transportation cost for using one vehicle to transport product units to depot $i \in I$ from its closest hub;

$e_{ij}$: transportation cost for using one vehicle from depot $i \in I$ to satellite $j \in J_i$;

$c_{jk}$: transportation cost between satellite $j \in J$ and customer $k \in K_j$.

The customer demands and the vehicle capacities are characterized by the following
parameters:

$n_k$: number of product units to deliver to customer $k \in K$;

$v_k$: volume of product units to deliver to customer $k \in K$;

$Q$: capacity (in number of product units) of one batch of products handled at any satellite;

$P$: volumetric capacity of a vehicle transporting product units to any depot, from its closest hub;

$R$: volumetric capacity of a vehicle transporting product units from any depot to any satellite.

The path-based model involves six families of decision variables. Three families are
associated with the design of the distribution system, while three families are required
to compute the total cost. First, $Y_i$ is a 0-1 variable which takes value 1 if depot $i \in I$
received product units from its closest hub; $T_{ij}$ is a 0-1 variable which takes value 1 if
product units are transported from depot $i \in I$ to satellite $j \in J$; $X_{ijk}$ is a 0-1 variable
which takes value 1 if product units are transported from depot $i \in I$ through satellite $j \in J_i$ to customer $k \in K_j$. Next, $V_i$ is an integer variable that represents the number of trucks used between depot $i \in I$ and its closest hub; $H_{ij}$ is an integer variable which represents the number of trucks used between depot $i \in I$ and satellite $j \in J_i$; $U_j$ is an integer variable which represents the number of batches of products handled at satellite $j \in J$.

$$ \min \sum_{i \in I} f_i Y_i + \sum_{i \in I} d_i V_i + \sum_{j \in J} g_j U_j + \sum_{i \in I} \sum_{j \in J_i} e_{ij} H_{ij} + \sum_{i \in I} \sum_{j \in J_i} \sum_{k \in K_j} c_{jk} X_{ijk} \quad (1) $$

$$ \sum_{i \in I} \sum_{j \in J_i} X_{ijk} = 1 \qquad \forall k \in K \qquad (2) $$

$$ \sum_{i \in I_j} T_{ij} \le 1 \qquad \forall j \in J \qquad (3) $$

$$ \sum_{j \in J_i} \sum_{k \in K_j} v_k X_{ijk} \le (\sum_{k \in K} v_k) Y_i \qquad \forall i \in I \qquad (4) $$

$$ X_{ijk} \le T_{ij} \qquad \forall i \in I, j \in J_i, k \in K_j \qquad (5) $$

$$ T_{ij} \le Y_i \qquad \forall i \in I, j \in J_i \qquad (6) $$

$$ \sum_{i \in I_j} \sum_{k \in K_j} n_k X_{ijk} \le Q U_j \qquad \forall j \in J \qquad (7) $$

$$ \sum_{j \in J_i} \sum_{k \in K_j} v_k X_{ijk} \le P V_i \qquad \forall i \in I \qquad (8) $$

$$ \sum_{k \in K_j} v_k X_{ijk} \le R H_{ij} \qquad \forall i \in I, j \in J_i \qquad (9) $$

$$ V_i \ge 0 \text{ and integer} \qquad \forall i \in I \qquad (10) $$

$$ U_j \ge 0 \text{ and integer} \qquad \forall j \in J \qquad (11) $$

$$ H_{ij} \ge 0 \text{ and integer} \qquad \forall i \in I, j \in J_i \qquad (12) $$

$$ Y_i \in \{0, 1\} \qquad \forall i \in I \qquad (13) $$

$$ T_{ij} \in \{0, 1\} \qquad \forall i \in I, j \in J_i, \qquad (14) $$

$$ X_{ijk} \in \{0, 1\} \qquad \forall i \in I, j \in J_i, k \in K_j \qquad (15) $$

In this model, the objective (1) represents the total cost for operating the distribution systems, i.e., the total handling cost and the total transportation costs between the facilities and customers. Constraints (2) enforce that each customer is served once. Constraints (3) ensure that each satellite is connected to at most one depot. Constraints (4), (5) and (6) are the linking constraints between the variables associated with the
design of the distribution system. Constraints (7) count the number of batches processed in each satellite. Constraints (8) and (9) count the number of vehicles operated between the hubs and the depots and between the depots and the satellites. Note that constraints (4) are redundant in this formulation.

[Gendron and Semet, 2009] compare different relaxations of the arc-based and path-based models on data obtained from the company. The real-life network has 93 depots, 320 satellites, and 701 customers. Based on this data, 32 instances were generated, varying the size of the network, the costs at the depots and satellites, and the capacity of the vehicles used between them. As expected, the linear relaxation of the path-based formulation is stronger than that of the arc-based formulation. However, when we consider binary relaxations, i.e., only the integer requirements on the variables $T_i, U_j$ and $H_{ij}$ are relaxed, both models provide the same lower bound. Observing that there is no fixed cost associated with a connection between a depot and a satellite, the binary relaxation of the path-based model can be modeled as a UFLP. Computational experiments were performed to confirm the theoretical results. On the real-life instances, the smallest gaps between the lower bounds computed with a particular relaxation and the best-known upper bounds were obtained by solving the path-based formulation with CPLEX with an average value of 1.7% after 2 hours of CPU Time. Based on these results and motivated by the real application, a heuristic approach was designed, which is described in the next section.

## 3 Heuristic approaches for the TUFLP-S with modular costs

Two heuristics were developed to solve the TUFLP-S with modular costs in Gendron et al. [2015]. First, we propose a variant of the Variable Neighborhood Search (VNS) metaheuristic [Mladenović and Hansen, 1997; Hansen et al., 2019] called the Multilayer VNS (MLVNS). Second, we adapt the classical slope scaling heuristic.

Variable Neighborhood Search (VNS) is a metaheuristic based on the exploitation of neighborhoods. First, VNS orders the neighborhood structures and then iteratively con-
siders them in that order. At each iteration, a new solution is obtained from the current best known solution by exploring the current neighborhood. Then, if a better solution is obtained, the local search continues considering the first neighborhood, otherwise the process continues from the best known solution by considering the next neighborhood.

MLVNS is a recursive variant of VNS based on partitioning the neighborhood structures into $l_{max}$ layers. Each layer $l$ has $k_{max}^l$ neighborhood structures associated with it. As in the classical VNS, the neighborhoods grow in complexity with the layer index. For each layer $l$, a VNS is invoked. First, a solution $x'$ is identified on a neighborhood $N_k(x)$ of the current solution $x$ until a first improving solution is found. To define the subset of $N_k(x)$, we use $m_{max}^k$ sub-neighborhood structures, resulting from a “natural” decomposition of $N_k$: $N_k^1, \dots, N_k^{m_{max}^k}$. At each step $m$, a solution $x'$ is generated by performing one move in $N_k^m(x)$. If $l=1$, we assume that evaluating this candidate solution is easy. Otherwise, when $l > 1$, MLVNS is called recursively up to layer $l-1$ to produce a final candidate solution $x'$. The generation of candidate solutions continues until $x'$ improves upon $x$ or $m > m_{max}^k$, i.e., all sub-neighborhood structures have been explored. Algorithm 1 summarizes the procedure.

In our implementation of MLVNS, we first generate an initial solution by a greedy constructive procedure. Then, we design a MLVNS algorithm that includes three layers. In layer 1, four simple neighborhood structures are invoked. They are based on exchange or closing operations. More precisely, they aim at: i) reassigning a customer from one satellite to another satellite; ii) reassigning a satellite from one depot to another depot; iii) closing a satellite; iv) closing a depot. Note that the exploration of the last two neighborhoods relies on calls to the first two. In layer 2, two neighborhood structures are used, which consist in opening a depot or a satellite presently closed. In each case, the subneighborhood structures are obtained by scanning the set of closed facilities, i.e., depots or satellites. The satellites/customers are then assigned to the closest facilities while keeping the assignment of at least one satellite/customer at each depot/satellite. The move is completed after applying MLVNS at layer 1 and the final configuration is determined. Layer 3 includes neighborhood structures designed to diversify the search
<caption>**Algorithm 1 MLVNS($x$: current solution, $l_{max}$: number of layers)**</caption>

1: **for** $l \leftarrow 1$ to $l_{max}$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;Order the $k_{max}^l$ neighborhood structures: $N_1, \dots, N_{k_{max}^l}$
3: &nbsp;&nbsp;&nbsp;&nbsp;**repeat**
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow 1$
5: &nbsp;&nbsp;&nbsp;&nbsp;**repeat**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Order (randomly) the $m_{max}^k$ sub-neighborhood structures: $N_k^1, \dots, N_k^{m_{max}^k}$
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m \leftarrow 1$
8: &nbsp;&nbsp;&nbsp;&nbsp;**repeat**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Generate $x' \in N_k^m(x)$
10: &nbsp;&nbsp;&nbsp;&nbsp;**if** $l > 1$ **then**
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;MLVNS($x', l-1$)
12: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
13: &nbsp;&nbsp;&nbsp;&nbsp;$m \leftarrow m + 1$
14: &nbsp;&nbsp;&nbsp;&nbsp;**until** ($f(x') < f(x)$) OR ($m > m_{max}^k$)
15: &nbsp;&nbsp;&nbsp;&nbsp;**if** $f(x') < f(x)$ **then**
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x \leftarrow x'$
17: &nbsp;&nbsp;&nbsp;&nbsp;**else**
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow 1$
19: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
20: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
21: &nbsp;&nbsp;&nbsp;&nbsp;**until** $k > k_{max}^l$
22: &nbsp;&nbsp;&nbsp;&nbsp;**until** some termination condition is satisfied
23: **end for**
through some cost perturbations. Four types of perturbations are considered: i) increas- ing the location cost for a selected subset of depots; ii) increasing the location cost for a selected subset of satellites; iii) increasing the transport costs from a depot belonging to a selected subset; iv) increase the transport costs from a satellite belonging to a selected subset. Once the costs have been modified, MLVNS at layer 2 is applied, then the costs are restored and MLVNS at layer 2 is applied again.

To assess the efficiency of the proposed MLVNS approach, we implemented an adapted variant of the slope scaling heuristic. As mentioned above (see Section 2), if we relax the integer requirements on the variables $T_i$, $U_j$ and $H_{ij}$ in the formulation of the TUFLP-S with modular costs, the model reduces to a UFLP model. The slope scaling heuristic consists in first determining an optimal solution to the UFLP. Given the values of vari- ables $Y_i$ and $X_{ijk}$, a solution of the TUFLP-S with modular costs is derived, and then improved by calling layer 1 of the MLVNS heuristic. Then, the costs in the UFLP model are updated so that the objective function value is equal to the value obtained with the formulation of the TUFLP-S with modular costs if the solution of the UFLP model does not change. For example, let $(\tilde{Y}, \tilde{X})$ be the solution obtained after performing layer 1. Let $i \in I$ be a depot such that : $\sum_{j \in J_i} \sum_{k \in K_j} v_k \tilde{X}_{ijk} > 0$. When the integer requirement on $V_i$ is relaxed, the corresponding value for $V_i$ is :

$$ \tilde{V}_i = \frac{1}{P} \sum_{j \in J_i} \sum_{k \in K_j} v_k \tilde{X}_{ijk} $$

Then the cost of $V_i$ in the binary relaxation of the model for the TUFLP-S with modular costs is updated as : $d_i[\tilde{V}_i]$. The derived UFLP formulation is solved, and the process is repeated until the objective values of the UFLP model and the solution with layer 1 of the MLVNS are equal, or until layer 1 provides the same solution for two consecutive calls, or until a time limit is reached.

[Gendron et al., 2015] compare the MLVNS with the slope scaling heuristic and with the solution of the path-based model for the TUFLP-S with modular costs with CPLEX on the testbed described in [Gendron and Semet, 2009] under a time limit of 24 minutes. Average gaps are computed relative to the lower bounds obtained by executing CPLEX
for up to two hours. First, we analyzed the contribution of each layer of the MLVNS to its performance. When layer 1 only is invoked, the average gap varies from 8.9% to 14.2 %. The gap ranges between 3.2% and 6.5% when layers 1 and 2 are applied. Last, when the MLVNS with the 3 layers is applied, it varies from 0.3 % for the smallest instances to 4.9% for the full network with the original costs of the company. The MLVNS performs better than both alternative approaches. CPLEX identifies poor solutions with an average gap that reaches 114.7%. The slope scaling heuristic is also outperformed, although the difference is weaker on the full network with an average gap of 5.1%.

In the path-based model, constraints (7),(8) and (9) are introduced to capture the complex modular cost function. The constraints associated with the design of the network defined by constraints (2),(3),(4), (5) and (6) correspond to those of the TUFLP-S, which is considered in what follows.

# 4 Formulations for the TUFLP-S

The TUFLP is defined on a graph where the set of nodes is partitioned into three subsets: *I*, *J* and *K*. *I* is associated with the potential depot locations, *J* with the potential satellite locations, and *K* with customer locations. The costs are divided into fixed costs and transportation costs. Fixed costs are $f_i$ for operating depot $i \in I$, $g_j$ for operating satellite $j \in J$ and $h_{ij}$ if satellite $j \in J$ is connected with depot $i \in I$. Finally, $c_{ijk}$ is the transportation cost from $i \in I$ to customer $k \in K$ via satellite $j \in J$.

In Barros and Labbé [1994], the authors introduce a model for the TUFLP, which four types of decision variables. First, $Y_i$ is a 0-1 variable which takes value 1 if depot $i \in I$ is open; $Z_j$ is a 0-1 variable which takes value 1 if satellite $j \in J$ is open; $T_{ij}$ is a 0-1 variable which takes value 1 if the connection between depot $i \in I$ to satellite $j \in J$ exists; $X_{ijk}$ is a 0-1 variable which takes value 1 if customer $k \in K$ is served from depot $i \in I$ via satellite $j \in J$. Then, the TUFLP can be stated as:
$$ (G) \min \sum_{i \in I} f_i Y_i + \sum_{j \in J} g_j Z_j + \sum_{i \in I} \sum_{j \in J} h_{ij} T_{ij} + \sum_{i \in I} \sum_{j \in J} \sum_{k \in K} c_{ijk} X_{ijk} \quad (16) $$

$$ \sum_{j \in J} \sum_{k \in K} X_{ijk} = 1 \qquad \forall k \in K \qquad (17) $$

$$ X_{ijk} \le T_{ij} \qquad \forall (i, j, k) \in I \times J \times K \qquad (18) $$

$$ \sum_{j \in J} X_{ijk} \le Y_i \qquad \forall (i, k) \in I \times K \qquad (19) $$

$$ \sum_{i \in I} X_{ijk} \le Z_j \qquad \forall (j, k) \in J \times K \qquad (20) $$

$$ Y_i \in \{0, 1\} \qquad \forall i \in I \qquad (21) $$

$$ Z_i \in \{0, 1\} \qquad \forall i \in I \qquad (22) $$

$$ T_{ij} \in \{0, 1\} \qquad \forall (i, j) \in I \times J \qquad (23) $$

$$ 0 \le X_{ijk} \le 1 \qquad \forall (i, j, k) \in I \times J \times K \qquad (24) $$

In this formulation, constraints (17) ensure that every customer of K is served once. Constraints (18,19,20) are the linking constraints between the decision variables to compute the total fixed costs associated with paths. The binary requirements on the $X_{ijk}$ variables are relaxed (24) since it is always possible to find an integer solution with the same cost when a solution with fractional values for $X_{ijk}$ variables is obtained.

In the TUFLP-S, an additional constraint forces each satellite to be connected to at most one depot. A weak model [Gendron et al., 2017] is obtained by adding the following constraints to (G):

$$ \sum_{i \in I} T_{ij} \le 1 \quad \forall j \in J \qquad (25) $$

Note that constraints (25) are the same constraints as (3) in (TCLDP). They can be viewed as generalized upper bound (GUB) constraints. This model can be improved by observing that $Z_j$ is equal to 1 if and only if it exists $i' \in I$ such that $T_{i'j}$ is equal to 1. This leads to the constraints:

$$ \sum_{i \in I} T_{ij} = Z_j \quad \forall j \in J \qquad (26) $$

which are stronger than (25) and lead to the so-called strong formulation (S) of the
TUFLP-S:

$$ (S) \min \sum_{i \in I} f_i Y_i + \sum_{j \in J} g_j Z_j + \sum_{(i,j) \in I \times J} h_{ij} T_{ij} + \sum_{(i,j,k) \in I \times J \times K} c_{ijk} X_{ijk} $$

$$ \sum_{(i,j) \in I \times J} X_{ijk} = 1 \qquad \forall k \in K $$

$$ \sum_{i \in I} T_{ij} = Z_j \qquad \forall j \in J $$

$$ X_{ijk} \le T_{ij} \qquad \forall (i,j,k) \in I \times J \times K $$

$$ \sum_{j \in J} X_{ijk} \le Y_i \qquad \forall (i,k) \in I \times K $$

$$ Y_i \in \{0, 1\} \qquad \forall i \in I $$

$$ Z_i \in \{0, 1\} \qquad \forall i \in I $$

$$ T_{ij} \in \{0, 1\} \qquad \forall (i,j) \in I \times J $$

$$ 0 \le X_{ijk} \le 1 \qquad \forall (i,j,k) \in I \times J \times K $$

Note that in (S), constraints (20) are removed since they are redundant. Then, we have the following result:

**Proposition 1 (Gendron et al. [2017])** The linear relaxation of (S) provides a lower bound which is greater than or equal to that of the weak model (G) with (25). This inequality can be strict.

Other formulations of the TUFLP-S can be easily obtained from (S). A model ($S_p$) with fewer decision variables can be obtained by projecting out variables $Z_j$ using constraints (26) and adding the constraints (25). Then, the integrality requirements on variable $Y_i$ can be relaxed in models (S) and ($S_p$) once we add the following constraints:

$$ T_{ij} \le Y_i \quad \forall (i,j) \in I \times J \qquad (27) $$

Thus, two alternative formulations ($S^C$) and ($S_p^C$) are derived. Note that the linear relaxation of the four models (S), ($S_p$), ($S^C$) and ($S_p^C$) provided the same bound [Gendron et al., 2017].
The four models plus the weak model were compared using CPLEX on a large testbed, including industrial instances derived from the data provided by the company for the TCLDP (see Section 2) and artificial instances obtained using the generator proposed by Kochetov and Ivanenko [2005] for the UFLP. These two sets have different characteristics with respect to the objective function. For industrial instances, the fixed costs associated with the satellites and with the connections between depots and satellites are set to zero, while the fixed costs associated with satellites have high values for artificial instances. Thus, the best results were obtained for industrial instances with ($S_p$) and for artificial instances with ($S^C$). The formulation ($S$) is a good compromise overall when a priority rule for branching is defined on variables $Y_i$ for industrial instances and on variables $Z_j$ for artificial instances.

# 5 Lagrangian-based methods for the TUFLP-S

We now present the exact and approximate solution methods proposed to the TUFLP-S [Gendron et al., 2016]. The main idea in developing the Lagrangian-based method was to recover a property already exploited in the relaxations of (TCLDP). As mentioned in Section 2, the binary relaxation of (TCLDP) can be modeled as a UFLP. Here, we aim to define a Lagrangian relaxation for a TUFLP-S formulation, which leads to define UFLPs as subproblems. To do so, we considered the model ($S_p^C$). ($S_p^C$) is the formulation ($S$) (see Section 4), where variables $Z_j$ are projected out, constraints (25) are added and the binary requirements on $Y_i$ are relaxed. When constraints (18) are relaxed, the Lagrangian subproblem $S_p^C(\lambda)$ is as follows, where ($\lambda$) are the Lagrange multipliers:
$$ (S_p^C(\lambda)) \min \sum_{i \in I} f_i Y_i + \sum_{(i,j) \in I \times J} (g_j + h_{ij} - \sum_{k \in K} \lambda_{ijk}) T_{ij} + \sum_{(i,j,k) \in I \times J \times K} (c_{ijk} + \lambda_{ijk}) X_{ijk} $$

$$ \sum_{(i,j) \in I \times J} X_{ijk} = 1 \qquad \forall k \in K $$

$$ \sum_{i \in I} T_{ij} \le 1 \qquad \forall j \in J $$

$$ \sum_{j \in J} X_{ijk} \le Y_i \qquad \forall (i,k) \in I \times K $$

$$ T_{ij} \le Y_i \qquad \forall (i,j) \in I \times J $$

$$ 0 \le Y_i \le 1 \qquad \forall i \in I $$

$$ T_{ij} \in \{0, 1\} \qquad \forall (i,j) \in I \times J $$

$$ 0 \le X_{ijk} \le 1 \qquad \forall (i,j,k) \in I \times J \times K $$

In the TUFLP-S, all variables are binary. Therefore, we can impose that variables $Y_i$ take 0-1 values. Then, we can relax the binary requirements on variable $T_{ij}$, since an integer solution can be easily derived from any optimal fractional solution. Finally, we can observe that all but one variable $X_{ijk}$ can be eliminated for a given pair $(i, k)$. Indeed, any optimal solution can be transformed into an equivalent solution where the only nonzero variable $X_{ijk}$ for $(i, k)$ is $X_{ij\ast k}$ with $c_{ij\ast k} + \lambda_{ij\ast k} = \min_{j \in J}(c_{ijk} + \lambda_{ijk})$. By defining new binary decision variables: $W_{ik} = X_{ij\ast k}$ and parameters : $d_{ik} = c_{ij\ast k} + \lambda_{ij\ast k}$ for all $(i, k) \in I \times K$, we get the compact formulation of the Lagrangian subproblem $S_p^C(\lambda)$:

$$ \min \sum_{i \in I} f_i Y_i + \sum_{(i,j) \in I \times J} (g_j + h_{ij} - \sum_{k \in K} \lambda_{ijk}) T_{ij} + \sum_{(i,k) \in I \times K} d_{ik} W_{ik} $$

$$ \sum_{i \in I} W_{ik} = 1 \qquad \forall k \in K $$

$$ \sum_{i \in I} T_{ij} \le 1 \qquad \forall j \in J \qquad (25) $$

$$ W_{ik} \le Y_i \qquad \forall i \in I $$

$$ T_{ij} \le Y_i \qquad \forall (i,j) \in I \times J $$

$$ Y_i \in \{0, 1\} \qquad \forall i \in I $$
$$0 \leq T_{ij} \leq 1$$

$$0 \leq W_{ik} \leq 1$$

$$\forall (i, j) \in I \times J$$

$$\forall (i, k) \in I \times K$$

This compact model is a formulation of an UFLP since GUB constraints (25) can be transformed into an equality by adding an artificial depot with zero assignment cost for each satellite.

Different approaches, such as the subgradient and bundle methods, were considered to solve the Lagrangian dual. However, they obtained tight lower bounds only when several iterations were performed, which was costly from a computational time point of view. Therefore, we used a simpler heuristic to obtain lower bounds. This approach consists in first solving the linear relaxation of ($S_p^C$), then the Lagrange multipliers are set to the values of the dual variables of constraints 18), and finally the Lagrangian subproblem $S_p^C(\lambda)$ is solved. This heuristic is fast and provides lower bounds, at least as good as the linear relaxation, since the Lagrangian subproblem does not have the integrality property.

A primal heuristic was devised from the solution of the Lagrangian subproblem obtained to obtain a feasible solution for the TUFLP-S. This heuristic alternatively solves two UFLPs, providing feasible solutions for the TUFLP-S. First, based on the solution of the subproblem, the set of open depots is fixed, and a UFLP, in which a connection between a depot and a satellite is viewed as a facility, is solved. Next, the connections between the satellites and the customers are fixed, and a UFLP, where the cost of assigning a satellite to a depot includes the total cost of serving the customers via the satellite from the depot, is solved. These two problems are solved to optimality repeatedly until no improvement of the best-known solution is observed.

These lower and upper bounds were embedded in an LP-based branch-and-bound algorithm. Two specific branching schemes were implemented. The first one is based on the GUB constraints (25). For satellite j, two branches are defined either $\sum_{i \in I} T_{ij} = 1$ or $\sum_{i \in I} T_{ij} = 0$. The second branching scheme implements a polytomic branching strategy based on the same constraint. As many child nodes as variables $T_{ij}$ in the sum are generated with $T_{ij} = 1$ and an additional node where the sum is set to zero. The most
effective branching strategy tested is to first branch with the GUB scheme and use the polytomic scheme if necessary.

Reliability branching [Achterberg et al., 2005] was implemented to perform the branch-
ing choice among the candidate branches. The main adaptation was related to the evalua-
tion of unreliable branches, i.e., branches for which the pseudocosts cannot be evaluated.
The unreliable branches are sorted according to the sum of absolute complementary
slackness violations in the Lagrangian subproblem. Satellites corresponding to the larger
values are considered first.

The efficiencies of the LP-based branch-and-bound algorithm and the primal heuristic
were assessed on the testbed described in Section 4. On the industrial instances, all
models presented in Section 4 provide the same gaps to optimal values at the root node.
The Lagrangian relaxation approach reduced these gaps significantly (by around 50%),
and the Lagrangian-based heuristic provides solutions with an optimality gap of 0.12% on
average. Considering the branch-and-bound algorithms, the best results were achieved
using only the polytomic branching. However, the approach was not competitive with the
solution of $S_p^C$ with CPLEX, even though the number of explored nodes was significantly
reduced (more than 14 times on average). On artificial instances, large gaps are obtained
at the root node with all models described in Section 4. Unfortunately, the Lagrangian
bounding procedure did not improve the values obtained with the linear relaxation of
$S_p^C$. However, the branch-and-bound algorithm with the branching strategy combining
the GUB scheme and the polytomic branching was able to solve all instances in less
than 1'600 seconds on average while CPLEX on formulation $S_p^C$ required up to 4'800
seconds on average. This difference in computation time is related to the number of
nodes explored in the branch-and-bound tree, which is more than 50 times larger for
the CPLEX solver. The implementation of reliability branching based on information
derived from solving the Lagrangian subproblem explains this reduction in the number
of nodes.
# 6 Conclusions

This article summarizes the research work of B. Gendron with his co-authors over a period of 15 years on the TUFLP-S. An industrial application in freight transportation drove the interest in this problem. This led to address the TUFLP-S with modular costs. The contribution was twofold. First, formulations of the classical TUFLP-S, on the one hand, and of the TUFLP-S with modular costs, on the other hand, were compared, both theoretically and experimentally, on a testbed including artificial and industrial instances. Second, heuristic and exact methods have been designed, and the computational results show their efficiencies. These results could be improved by developing an adaptive branching strategy choice for solving the models as well as for the specialized branch-and-bound algorithm. The latter would benefit to a great extent from a specialized algorithm for solving the UFLP. Regarding the industrial application, the solutions obtained with the MLVNS are within 5% of the best lower bounds obtained by solving the path-based model. A promising research project is the development of exact algorithms based on decomposition methods applied to the path-based model. Note that one of the main challenges lies in the modular cost structure of the objective function. We have shown that industrial instances can be solved easily when considered with the classical linear cost function.

The planned research avenues include the study of extensions to the multi-level facility location problems and to the multiperiod version. An interesting extension of the TUFLP-S with modular costs relies on the observation that the distribution network results from contracts between the company and carriers for longhaul and last-mile delivery services. Designing the network to meet time-varying demand while considering contracting with carriers is a promising research topic addressed recently in Clavijo López et al. [2024].
References

Aardal, K., Labbé, M., Leung, J., and Queyranne, M. (1996). On the two-level uncapacitated facility location problem. *INFORMS Journal on Computing*, 8(3):289–301.

Achterberg, T., Koch, T., and Martin, A. (2005). Branching rules revisited. *Operations Research Letters*, 33(1):42–54.

Amrani, H., Martel, A., Zufferey, N., and Makeeva, P. (2011). A variable neighborhood search heuristic for the design of multicommodity production-distribution networks with alternative facility configurations. *OR spectrum*, 33:989–1007.

Baïou, M. and Barahona, F. (2014). A polyhedral study of a two level facility location model. *RAIRO-Operations Research*, 48(2):153–165.

Barros, A. I., Dekker, R., and Scholten, V. (1998). A two-level network for recycling sand: a case study. *European Journal of Operational Research*, 110(2):199–214.

Barros, I. (1998). *Discrete and Fractional Programming Techniques for Location Models*. Kluwer Academic Publishers, Dordrecht, Boston, London.

Barros, I. and Labbé, M. (1994). A general model for the uncapacitated facility and depot location problem. *Location Science*, 2:173–191.

Bloemhof-Ruwaard, J. M., Salomon, M., and Van Wassenhove, L. N. (1994). On the co-ordination of product and by-product flows in two-level distribution networks: Model formulations and solution procedures. *European Journal of Operational Research*, 79(2):325–339.

Bloemhof-Ruwaard, J. M., Salomon, M., and Van Wassenhove, L. N. (1996). The capacitated distribution and waste disposal problem. *European Journal of Operational Research*, 88(3):490–503.

Chardaire, P., Lutton, J.-L., and Sutter, A. (1999). Upper and lower bounds for the two-level simple plant location problem. *Annals of Operations Research*, 86(0):117–140.
Clavijo López, C., Crama, Y., Pironet, T., and Semet, F. (2024). Multi-period distribution networks with purchase commitment contracts. *European Journal of Operational Research*, 312(2):556-572.

de Oliveira, P. B., Contreras, I., de Camargo, R. S., and de Miranda Júnior, G. (2020). A comparison of separation routines for benders optimality cuts for two-level facility location problems. *Expert Systems with Applications*, 141:112928.

Gao, L.-L. and Robinson Jr, E. P. (1992). A dual-based optimization procedure for the two-echelon uncapacitated facility location problem. *Naval Research Logistics (NRL)*, 39(2):191-212.

Gendron, B., Khuong, P.-V., and Semet, F. (2015). Multilayer variable neighborhood search for two-level uncapacitated facility location problems with single assignment. *Networks*, 66(3):214-234.

Gendron, B., Khuong, P.-V., and Semet, F. (2016). A lagrangian-based branch-and-bound algorithm for the two-level uncapacitated facility location problem with single-assignment constraints. *Transportation Science*, 50(4):1286-1299.

Gendron, B., Khuong, P.-V., and Semet, F. (2017). Comparison of formulations for the two-level uncapacitated facility location problem with single assignment constraints. *Computers & Operations Research*, 86:86-93.

Gendron, B. and Semet, F. (2009). Formulations and relaxations for a multi-echelon capacitated location-distribution problem. *Computers & Operations Research*, 36(5):1335-1355.

Gendron, B., Semet, F., and Strozyk, C. (2002). Adaptive Distribution Systems. In *TRANSTECH: Transport Technology Product & Process Innovation Management*, pages 36-42. Presses Universitaires de Valenciennes, Valenciennes.

Hansen, P., Mladenović, N., Brimberg, J., and Pérez, J. A. M. (2019). *Variable neighborhood search*. Springer.
Huang, Y., Chen, C.-W., and Fan, Y. (2010). Multistage optimization of the supply chains of biofuels. *Transportation Research Part E: Logistics and Transportation Review*, 46(6):820–830.

Ignacio, A. A. V., Ferreira Filho, V. J. M., and Galvão, R. D. (2008). Lower and upper bounds for a two-level hierarchical location problem in computer networks. *Computers & Operations Research*, 35(6):1982–1998.

Kaufman, L., Eede, M. V., and Hansen, P. (1977). A plant and warehouse location problem. *Journal of the Operational Research Society*, 28(3):547–554.

Klose, A. and Drexel, A. (2005). Facility location models for distribution system design. *European Journal of Operational Research*, 162(1):4–29.

Kochetov, Y. and Ivanenko, D. (2005). Computationally difficult instances for the uncapacitated facility location problem. *Metaheuristics: Progress as real problem solvers*, 32:351–367.

Krarup, J. and Pruzan, P. M. (1983). The simple plant location problem: Survey and synthesis. *European Journal of Operational Research*, 12(36-81):41.

Kratica, J., Dugošija, D., and Savić, A. (2014). A new mixed integer linear programming model for the multi level uncapacitated facility location problem. *Applied Mathematical Modelling*, 38(7):2118–2129.

Landete, M. and Marín, A. (2009). New facets for the two-stage uncapacitated facility location polytope. *Computational Optimization and Applications*, 44(3):487.

Marín, A. (2007). Lower bounds for the two-stage uncapacitated facility location problem. *European Journal of Operational Research*, 179(3):1126–1142.

Marín, A. and Pelegrín, B. (1999). Applying lagrangian relaxation to the resolution of two-stage location problems. *Annals of Operations Research*, 86(0):179–198.
Melo, M. T., Nickel, S., and Saldanha-Da-Gama, F. (2009). Facility location and supply chain management—a review. *European Journal of Operational Research*, 196(2):401–412.

Mitropoulos, P., Giannikos, I., and Mitropoulos, I. (2009). Exact and heuristic approaches for the locational planning of an integrated solid waste management system. *Operational Research*, 9:329–347.

Mladenović, N. and Hansen, P. (1997). Variable neighborhood search. *Computers & Operations Research*, 24(11):1097–1100.

Narula, S. C. and Ogbu, U. I. (1979). An hierarchal location—allocation problem. *Omega*, 7(2):137–143.

Ortiz-Astorquiza, C., Contreras, I., and Laporte, G. (2018). Multi-level facility location problems. *European Journal of Operational Research*, 267(3):791–805.

Pirkul, H. and Jayaraman, V. (1996). Production, transportation, and distribution planning in a multi-commodity tri-echelon system. *Transportation Science*, 30(4):291–302.

Pirkul, H. and Jayaraman, V. (1998). A multi-commodity, multi-plant, capacitated facility location problem: formulation and efficient heuristic solution. *Computers & Operations Research*, 25(10):869–878.

Ro, H.-b. and Tcha, D.-w. (1984). A branch and bound algorithm for the two-level uncapacitated facility location problem with some side constraints. *European Journal of Operational Research*, 18(3):349–358.

Şahin, G. and Süral, H. (2007). A review of hierarchical facility location models. *Computers & Operations Research*, 34(8):2310–2331.

Thanh, P. N., Péton, O., and Bostel, N. (2010). A linear relaxation-based heuristic approach for logistics network design. *Computers & Industrial Engineering*, 59(4):964–975.
Zhang, J. (2006). Approximating the two-level facility location problem via a quasi-greedy approach. *Mathematical Programming*, 108:159–176.