# A Catalog of Formulations for the Network Pricing Problem

Quang Minh Bui  CIRRELT and Département d’informatique et de recherche
opérationnelle, Université de Montréal  Bernard Gendron  CIRRELT and
Département d’informatique et de recherche opérationnelle, Université de
Montréal  Margarida Carvalho  CIRRELT and Département d’informatique et de
recherche opérationnelle, Université de Montréal

###### Abstract

We study the network pricing problem where the leader maximizes their revenue
by determining the optimal amounts of tolls to charge on a set of arcs, under
the assumption that the followers will react rationally and choose the
shortest paths to travel. Many distinct single-level reformulations to this
bilevel optimization program have been proposed, however, their relationship
has not been established. In this paper, we aim to build a connection between
those reformulations and explore the combination of the path representation
with various modeling options, allowing us to generate 12 different
reformulations of the problem. Moreover, we propose a new path enumeration
scheme, path-based preprocessing, and hybrid framework to further improve
performance and robustness when solving the final model. We provide numerical
results, comparing all the derived reformulations and confirming the
efficiency of the novel dimensionality reduction procedures.

##  1 Introduction

The network pricing problem (NPP) is a bilevel optimization program involving
two parties: a leader and multiple followers. First, the leader sets the
prices of several arcs (called tolled arcs) in a network. Afterwards, the
followers choose the optimal paths in the network subject to the prices set by
the leader. The leader’s objective is to maximize their profit, while the goal
of the followers is to minimize their own costs of transit across the network.
High prices do not always mean more profit for the leader if only a few
followers can afford them. Thus, the leader must seek a balance between the
prices and the demands of the followers reacting to those prices. In practice,
the leader could be a highway authority and the followers could be groups of
the population residing in different neighborhoods or cities.

###  1.1 Related Literature

The NPP was first introduced in Labbé et al. [13] as a single-commodity
problem (one follower), and later extended to the multi-commodity version
(multiple followers) by Brotcorne et al. [5]. The problem is proven to be NP-
hard even when there is only one follower [15]. The first mixed integer linear
program (MILP) formulation was also introduced in [5]. Since then, many
improvements to the MILP of [5] were proposed such as the shortest-path graph
model (SPGM) [17], preprocessing techniques [1], path-based models [1, 8],
valid inequalities, and tight bounds for big-M parameters [7]. Other methods
were also explored, including multipath enumeration [3] and tabu search [2].
Variants of the network pricing problem include the joint network design and
pricing problem [6], the complete toll NPP (clique pricing problem) [12], and
the logit NPP [10].

The NPP has a natural bilevel optimization formulation. There are general-
purpose methods to solve bilevel optimization problems [9, 16]. However, due
to their complexity, they are usually considered per case to exploit specific
problem structures. A common approach to solve a bilevel optimization problem
is to convert it to a single-level reformulation. For the NPP, the follower
problems are linear. Thus, there are three methods generally used for this
conversion based on optimality conditions: strong duality, complementary
slackness, and value function. In the first two methods, the primal and the
dual constraints of the follower problems are included in the single-level
reformulation and then accompanied by either strong duality or complementary
slackness constraints. In the third method, only the primal constraints are
needed and the optimality of the follower problems is ensured by the value
function constraints.

Strong duality was used in the original formulation proposed in [5]. The
discovery of the path representation led to the use of complementary slackness
in the path-based formulations developed in [8]. Value function formulations
were mentioned in [11] and [1]. However, these models were presented as
separate formulations with little connection in between.

###  1.2 Contribution and Paper Organization

In the first part of this work, our aim is to explore the links between all of
the previously mentioned single-level reformulations of the NPP. First, we
show that in strong duality and complementary slackness reformulations, the
path representation can be applied not only to the primal follower problems,
but also to the dual problems. This allows us to show that the value function
method is a special case of strong duality when the dual problems are written
in the path representation. Second, we provide a systematic way to formalize a
reformulation by breaking it down into components with different options for
each component. This general method can also be used to categorize all the
models mentioned above and explain the connection between them.

The second part of this work is dedicated to the path enumeration process.
This is important because a well-performed path-based reformulation requires
an efficient path enumeration algorithm. Existing methods must enumerate
either all the paths [1] or a set of relevant paths but require multiple calls
to a linear solver [8]. Based on previous discoveries on the properties of
redundant paths in [1, 8], we develop a new path enumeration process that can
enumerate relevant paths without employing a linear solver. Besides that,
based on this set of relevant paths, we also derive a new preprocessing method
that can be applied to arc-based reformulations. This preprocessing method is
crucial if for a fair comparison between the reformulations, since path
enumeration is also counted as a preprocessing method for the path-based
reformulations.

Regardless of the efficiency of the path enumeration process, the potential
size of the set of all relevant paths is exponential, which may require an
enormous amount of time to enumerate them. In the last part of this work, we
introduce a method to circumvent this problem by taking advantage of the
multi-commodity nature of the NPP. By deciding which commodities are worth
applying path enumeration and mixing different reformulations into one, we can
save time by spending less time enumerating very large sets of paths and more
time solving the actual optimization model. In this paper, we refer to this
solution as the hybrid framework for multi-commodity problems.

This paper is structured as follows. In Section 2, we describe the NPP,
namely, its bilevel optimization program, and explain the general method to
formalize a reformulation. Section 3 shows how path enumeration is performed
in the novel preprocessing method. The hybrid framework is introduced in
Section 4. Section 5 presents computational results, including a comparison of
all formulations and several experiments to prove the efficiency of the new
preprocessing method and the hybrid framework.

##  2 The Network Pricing Problem

In this section, we provide the bilevel programming formulation for the
network pricing problem in Section 2.1 and a systematic way to generate
single-level reformulations through the combination of different components in
Section 2.2. The conversion from the bilevel formulation to single-level
reformulations produces bilinear terms. The linearization of these bilinear
terms is the topic of discussion in Section 2.3.

###  2.1 Problem Formulation

Let us consider a graph $G=(V,A)$ where $A$ is partitioned into a set of toll
arcs $A_{1}$ and a set of toll-free arcs $A_{2}$. Each arc $a\in A$ has an
initial cost $c_{a}>0$. Let $K$ be the set of commodities (O-D pairs) and
$o^{k},d^{k},\eta^{k}$ be the origin, the destination, and the demand of
commodity $k\in K$, respectively. We define the following variables:

  * •

$T_{a},a\in A_{1}$: the toll of arc $a$;

  * •

$x_{a}^{k},a\in A_{1},k\in K$: the flow of commodity $k$ on toll arc $a$;

  * •

$y_{a}^{k},a\in A_{2},k\in K$: the flow of commodity $k$ on toll-free arc $a$.

The network pricing problem is then formulated as a bilevel program:

$\displaystyle(\text{NPP})\quad\max_{T\geq 0,x,y}$ $\displaystyle\sum_{k\in
K}\sum_{a\in A_{1}}\eta^{k}T_{a}x_{a}^{k}$ $\displaystyle\text{s.t. }\forall
k\in K$ $\displaystyle\left\\{\begin{aligned}
(x^{k},y^{k})\in\arg\min_{\hat{x},\hat{y}}&\sum_{a\in
A_{1}}(c_{a}+T_{a})\hat{x}_{a}+\sum_{a\in A_{2}}c_{a}\hat{y}_{a},\\\
\text{s.t. }&\sum_{a\in A^{+}_{1}(i)}\hat{x}_{a}+\sum_{a\in
A^{+}_{2}(i)}\hat{y}_{a}-\sum_{a\in A^{-}_{1}(i)}\hat{x}_{a}-\sum_{a\in
A^{-}_{2}(i)}\hat{y}_{a}=b_{i}^{k}&&i\in V,\\\ &\hat{x}_{a}\in\\{0,1\\}&&a\in
A_{1},\\\ &\hat{y}_{a}\in\\{0,1\\}&&a\in A_{2},\end{aligned}\right.$

where $b_{i}^{k}=1$ if $i=o^{k}$, $-1$ if $i=d^{k}$, and $0$ otherwise. For a
node $i$, the set of its outgoing tolled (toll-free) arcs is $A^{+}_{1}(i)$
($A^{+}_{2}(i)$) and the set of its incoming tolled (toll-free) arcs is
$A^{-}_{1}(i)$ ($A^{-}_{2}(i)$). The leader controls the toll prices $T_{a}$,
while each follower $k\in K$ decides the taken paths by setting $x_{a}^{k}$
and $y_{a}^{k}$ equal to 1 if the arc belongs to the path, and 0 otherwise. We
consider the optimistic version of this bilevel problem, _i.e._ , if a
follower has multiple optimal solutions with respect to the prices set by the
leader, then they will choose the solution that benefits the leader the most.
For each commodity, we assume that there exists at least a toll-free path (a
path without tolled arcs). Otherwise, the leader could gain infinite revenue
by exploiting that particular O-D pair.

###  2.2 Generalized Single-level Reformulations

A generalized single-level reformulation of the network pricing problem
consists of three components:

  * •

Primal representation (arc or path);

  * •

Dual representation (arc or path);

  * •

Optimality condition (strong duality or complementary slackness);

First, the primal and the dual representations of the follower problems need
to be chosen. There are two options for each: the arc representation and the
path representation. Consider the follower problem for commodity $k$:

$\displaystyle(\text{PA}^{k})$ $\displaystyle\min_{x^{k},y^{k}}$
$\displaystyle\sum_{a\in A_{1}}(c_{a}+T_{a})x_{a}^{k}+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}$ s.t. $\displaystyle\sum_{a\in
A^{+}_{1}(i)}x_{a}^{k}+\sum_{a\in A^{+}_{2}(i)}y_{a}^{k}-\sum_{a\in
A^{-}_{1}(i)}x_{a}^{k}-\sum_{a\in A^{-}_{2}(i)}y_{a}^{k}=b_{i}^{k}$
$\displaystyle i\in V,$ (1) $\displaystyle x_{a}^{k}\in\\{0,1\\}$
$\displaystyle a\in A_{1},$ $\displaystyle y_{a}^{k}\in\\{0,1\\}$
$\displaystyle a\in A_{2}.$

This is called the _primal-arc representation_ of the follower problem.
Because the set of constraints of the follower problem forms a totally
unimodular matrix, integrality conditions for $x_{a}^{k}$ and $y_{a}^{k}$ can
be dropped. This enables us to write the follower problem in the _dual-arc
representation_ :

$\displaystyle(\text{DA}^{k})$ $\displaystyle\max_{\lambda^{k}}\ $
$\displaystyle\lambda^{k}_{o^{k}}-\lambda^{k}_{d^{k}}$ s.t.
$\displaystyle\lambda^{k}_{i}-\lambda^{k}_{j}\leq c_{a}+T_{a}$ $\displaystyle
a\equiv(i,j)\in A_{1},$ (2) $\displaystyle\lambda^{k}_{i}-\lambda^{k}_{j}\leq
c_{a}$ $\displaystyle a\equiv(i,j)\in A_{2}.$ (3)

In [1, 8], instead of writing the follower problem using arc-flow
$x_{a}^{k},y_{a}^{k}$, the authors replaced these variables with $z^{k}_{p}$
representing the path-flow or the selection of paths $p$ in the set of all
elementary paths $P^{k}$ which is finite. We will refer this as the _primal-
path representation_ :

$\displaystyle(\text{PP}^{k})$ $\displaystyle\min_{z^{k}}\ $
$\displaystyle\sum_{p\in P^{k}}\left(\sum_{a\in
A}\delta_{a}^{p}c_{a}+\sum_{a\in A_{1}}\delta_{a}^{p}T_{a}\right)z^{k}_{p}$
s.t. $\displaystyle\sum_{p\in P^{k}}z^{k}_{p}=1,$ (4) $\displaystyle
z^{k}_{p}\in\\{0,1\\}$ $\displaystyle p\in P^{k}.$

In the above formulation, $\delta_{a}^{p}=1$ if the path $p$ includes the arc
$a$ and $\delta_{a}^{p}=0$ if not. Once again, the constraints form a totally
unimodular matrix, hence $z^{k}_{p}$ is not required to be binary and thus the
follower problem also has a _dual-path representation_ :

$\displaystyle(\text{DP}^{k})$ $\displaystyle\max_{L^{k}}\ $ $\displaystyle
L^{k}$ s.t. $\displaystyle L^{k}\leq\sum_{a\in
A}\delta_{a}^{p}c_{a}+\sum_{a\in A_{1}}\delta_{a}^{p}T_{a}$ $\displaystyle
p\in P^{k}.$ (5)

In total, there are 4 combinations of primal-dual representations (arc-arc,
arc-path, path-arc, and path-path). Note that the primal and the dual
representations are not required to match each other, so primal-arc could also
be paired with dual-path.

After deciding the representations, we need to connect the primal and the dual
representations by either strong duality or complementary slackness
constraints. Strong duality simply connects the primal and the dual objective
functions together. Below are the strong duality constraints for the arc-arc,
arc-path, path-arc, and path-path combinations, respectively:

$\displaystyle\sum_{a\in A_{1}}(c_{a}+T_{a})x_{a}^{k}+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}$ $\displaystyle=\lambda^{k}_{o^{k}}-\lambda^{k}_{d^{k}},$
(6) $\displaystyle\sum_{a\in A_{1}}(c_{a}+T_{a})x_{a}^{k}+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}$ $\displaystyle=L^{k},$ (7) $\displaystyle\sum_{p\in
P^{k}}\left(\sum_{a\in A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}\right)z^{k}_{p}$
$\displaystyle=\lambda^{k}_{o^{k}}-\lambda^{k}_{d^{k}},$ (8)
$\displaystyle\sum_{p\in P^{k}}\left(\sum_{a\in
A}\delta_{a}^{p}c_{a}+\sum_{a\in A_{1}}\delta_{a}^{p}T_{a}\right)z^{k}_{p}$
$\displaystyle=L^{k}.$ (9)

Complementary slackness matches bounded variables in the primal with
inequality constraints in the dual. The pairing is simple if the primal and
the dual have the same representation (either arc-arc or path-path):

$\displaystyle(c_{a}+T_{a}-\lambda^{k}_{i}+\lambda^{k}_{j})x_{a}^{k}$
$\displaystyle=0$ $\displaystyle a\equiv(i,j)\in A_{1},$ (10)
$\displaystyle(c_{a}-\lambda^{k}_{i}+\lambda^{k}_{j})y_{a}^{k}$
$\displaystyle=0$ $\displaystyle a\equiv(i,j)\in A_{2},$ (11)
$\displaystyle\left(\sum_{a\in A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}-L^{k}\right)z^{k}_{p}$ $\displaystyle=0$
$\displaystyle p\in P^{k}.$ (12)

However, if the primal and the dual have different representations, variables
in the primal representation must be converted to their dual counterparts. For
the path-arc combination, conversion from $z$ to $(x,y)$ is done by using the
identities $x_{a}^{k},y_{a}^{k}=\sum_{p\in P^{k}}\delta_{a}^{p}z^{k}_{p}$:

$\displaystyle(c_{a}+T_{a}-\lambda^{k}_{i}+\lambda^{k}_{j})\left(\sum_{p\in
P^{k}}\delta_{a}^{p}z^{k}_{p}\right)$ $\displaystyle=0$ $\displaystyle
a\equiv(i,j)\in A_{1},$ (13)
$\displaystyle(c_{a}-\lambda^{k}_{i}+\lambda^{k}_{j})\left(\sum_{p\in
P^{k}}\delta_{a}^{p}z^{k}_{p}\right)$ $\displaystyle=0$ $\displaystyle
a\equiv(i,j)\in A_{2}.$ (14)

For the arc-path combination, the conversion from $(x,y)$ to $z$ is more
complicated. We will use the definition that $z^{k}_{p}=1$ if and only if
$x_{a}^{k}=1$ and $y_{a}^{k}=1$ for all arcs $a$ belonging to the path $p$
(given that $(x,y)$ produce an elementary path). Mathematically,
$z^{k}_{p}=\prod_{a\in A_{1}\mid\delta_{a}^{p}=1}x_{a}^{k}\prod_{a\in
A_{2}\mid\delta_{a}^{p}=1}y_{a}^{k}$. The complementary slackness constraint
for the arc-path combination is:

$\displaystyle\left(\sum_{a\in A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}-L^{k}\right)\prod_{\begin{subarray}{c}a\in A_{1}\\\
\delta_{a}^{p}=1\end{subarray}}x_{a}^{k}\prod_{\begin{subarray}{c}a\in
A_{2}\\\ \delta_{a}^{p}=1\end{subarray}}y_{a}^{k}$ $\displaystyle=0$
$\displaystyle p\in P^{k}.$ (15)

An important rule of the conversion is that it must cover all possible
solutions in the primal representation. Otherwise, the optimality condition
can be dodged by choosing the primal solution that is not covered. If we leave
the graph and $P^{k}$ as inputted, this rule is satisfied. However, in Section
3, we will introduce methods that reduce the size of the graph and $P^{k}$,
thus requiring additional attention. We will discuss this later in Section
3.2.

For completeness, accordingly with the choice of the primal representation,
here are the objective function of the single-level formulations:

$\displaystyle\sum_{k\in K}\sum_{a\in A_{1}}\eta^{k}T_{a}x_{a}^{k},$ (16)
$\displaystyle\sum_{k\in K}\sum_{a\in A_{1}}\sum_{p\in
P^{k}}\eta^{k}\delta_{a}^{p}T_{a}z^{k}_{p}.$ (17)

Based on the primal, the dual, and the optimality condition, we can categorize
all the single-level formulations as in Table 1. The variables, constraints
and objective functions of these formulations are listed in Table 2. The
variables are bounded implicitly ($x,y,z,T$ are non-negative; $\lambda,L$ are
unbounded).

Dual Primal | Arc | Path
---|---|---
Arc | Standard (STD) |  Path-Arc Standard (PASTD)
Path | Value Function (VF) |  Path Value Function (PVF)
(a) Strong duality

Dual Primal | Arc | Path
---|---|---
Arc |  Complementary Slackness (CS) |  Path-Arc Complementary Slackness (PACS)
Path |  Value Function Complementary Slackness (VFCS) |  Path Complementary Slackness (PCS)
(b) Complementary slackness

Table 1: Complete map of all single-level reformulations. |  |  | Constraints
---|---|---|---
Label | Variables | Objective | Primal | Dual | Opt. Cond.
Strong duality |  |
STD | $T,x,y,\lambda$ | (16) | (1) | (2), (3) | (6)
VF | $T,x,y,L$ | (16) | (1) | (5) | (7)
PASTD | $T,z,\lambda$ | (17) | (4) | (2), (3) | (8)
PVF | $T,z,L$ | (17) | (4) | (5) | (9)
Complementary slackness |  |
CS | $T,x,y,\lambda$ | (16) | (1) | (2), (3) | (10), (11)
VFCS | $T,x,y,L$ | (16) | (1) | (5) | (15)
PACS | $T,z,\lambda$ | (17) | (4) | (2), (3) | (13), (14)
PCS | $T,z,L$ | (17) | (4) | (5) | (12)
Table 2: List of variables, constraints, and objective functions of all
single-level reformulations.

There are some notable reformulations in the list. The most canonical way to
convert a bilevel linear program to its single-level version is the standard
formulation (STD):

$\displaystyle(\text{STD})$ $\displaystyle\max\ $ $\displaystyle\sum_{k\in
K}\sum_{a\in A_{1}}\eta^{k}T_{a}x_{a}^{k}$ s.t. $\displaystyle\sum_{a\in
A^{+}_{1}(i)}x_{a}^{k}+\sum_{a\in A^{+}_{2}(i)}y_{a}^{k}-\sum_{a\in
A^{-}_{1}(i)}x_{a}^{k}-\sum_{a\in A^{-}_{2}(i)}y_{a}^{k}=b_{i}^{k}$
$\displaystyle k\in K,i\in V,$
$\displaystyle\lambda^{k}_{i}-\lambda^{k}_{j}\leq c_{a}+T_{a}$ $\displaystyle
k\in K,a\equiv(i,j)\in A_{1},$
$\displaystyle\lambda^{k}_{i}-\lambda^{k}_{j}\leq c_{a}$ $\displaystyle k\in
K,a\equiv(i,j)\in A_{2},$ $\displaystyle\sum_{a\in
A_{1}}(c_{a}+T_{a})x_{a}^{k}+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}=\lambda^{k}_{o^{k}}-\lambda^{k}_{d^{k}}$ $\displaystyle
k\in K,$ $\displaystyle x_{a}^{k}\geq 0$ $\displaystyle k\in K,a\in A_{1},$
$\displaystyle y_{a}^{k}\geq 0$ $\displaystyle k\in K,a\in A_{2},$
$\displaystyle T_{a}\geq 0$ $\displaystyle a\in A_{1}.$

This formulation has been referred to in most papers in the literature [11, 8,
1, 7], including the first paper on the NPP [5]. It only uses arc
representations, so the need for path enumeration is eliminated. Its
linearized version (Section 2.3) only requires $x$ to be integer. Furthermore,
the number of arcs is fixed, hence, it is reliable compared to the path-based
formulations whose sizes depend on the number of paths of each commodity which
is unknown. Overall, the standard formulation is the most straightforward
method to solve the NPP.

The second reformulation in which we are interested is the value function
formulation (VF):

$\displaystyle(\text{VF})$ $\displaystyle\max\ $ $\displaystyle\sum_{k\in
K}\sum_{a\in A_{1}}\eta^{k}T_{a}x_{a}^{k}$ s.t. $\displaystyle\sum_{a\in
A^{+}_{1}(i)}x_{a}^{k}+\sum_{a\in A^{+}_{2}(i)}y_{a}^{k}-\sum_{a\in
A^{-}_{1}(i)}x_{a}^{k}-\sum_{a\in A^{-}_{2}(i)}y_{a}^{k}=b_{i}^{k}$
$\displaystyle k\in K,i\in V,$ $\displaystyle L^{k}\leq\sum_{a\in
A}\delta_{a}^{p}c_{a}+\sum_{a\in A_{1}}\delta_{a}^{p}T_{a}$ $\displaystyle
k\in K,p\in P^{k},$ $\displaystyle(c_{a}+T_{a})x_{a}^{k}+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}=L^{k}$ $\displaystyle k\in K,$ $\displaystyle
x_{a}^{k}\geq 0$ $\displaystyle k\in K,a\in A_{1},$ $\displaystyle
y_{a}^{k}\geq 0$ $\displaystyle k\in K,a\in A_{2},$ $\displaystyle T_{a}\geq
0$ $\displaystyle a\in A_{1}.$

If we combine the second and the third constraints, a new constraint emerges:

$\sum_{a\in A_{1}}(c_{a}+T_{a})x_{a}^{k}+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}\leq\sum_{a\in A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}\hskip 28.45274ptp\in P^{k}.$ (18)

The left hand side of Eq. (18) is the cost of the current path, while the
right hand side is the cost of path $p$. Eq. (18) means that the cost of the
current path must not surpass the cost of any path in $P^{k}$, which restricts
$(x,y)$ to choose the path with least cost. This is called the value function
constraint and it is also a popular method to generate the single-level
formulation of a bilevel problem. Here, we have shown that the value function
method is just a special case of strong duality when we write the dual
representation of the follower problems in the solution space (dual-path in
the case of the NPP). The value function formulation for the NPP was
previously mentioned in [11].

The last formulation in which we focus is the path complementary slackness
formulation (PCS):

$\displaystyle(\text{PCS})$ $\displaystyle\max\ $ $\displaystyle\sum_{k\in
K}\sum_{a\in A_{1}}\sum_{p\in P^{k}}\eta^{k}\delta_{a}^{p}T_{a}z^{k}_{p}$ s.t.
$\displaystyle\sum_{p\in P^{k}}z^{k}_{p}=1$ $\displaystyle k\in K,$
$\displaystyle L^{k}\leq\sum_{a\in A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}$ $\displaystyle k\in K,p\in P^{k},$
$\displaystyle\left(\sum_{a\in A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}-L^{k}\right)z^{k}_{p}=0$ $\displaystyle k\in K,p\in
P^{k},$ $\displaystyle z^{k}_{p}\geq 0$ $\displaystyle k\in K,p\in P^{k},$
$\displaystyle T_{a}\geq 0$ $\displaystyle a\in A_{1}.$

This formulation is the path-based model of Didi-Biha et al. [8]. It is the
complete opposite of the standard formulation. It only uses path
representations, hence its performance heavily relies on the number of paths
in $P^{k}$. Methods for reducing the size of $P^{k}$ will be discussed in
Section 3. In this formulation, only the selection of $p$ is important and any
information on the structure of the graph is discarded. The advantage of this
formulation is its simplicity, especially when the size of $P^{k}$ is small.

###  2.3 Linearization

In all eight single-level formulations, there are bilinear terms (or in the
case of VFCS, multilinear terms). In order to solve these formulations using
MILP solvers, these terms need to be linearized. We will consider formulations
using strong duality and those using complementary slackness separately.

####  2.3.1 Linearization of Strong Duality Formulations

In the strong duality formulations, the bilinear terms arise from the leader’s
revenue, which are either $\sum_{a\in A_{1}}T_{a}x_{a}^{k}$ or $\sum_{a\in
A_{1}}\sum_{p\in P^{k}}\delta_{a}^{p}T_{a}z^{k}_{p}$. Because $x_{a}^{k}$ and
$z^{k}_{p}$ can only take 0 and 1 as their values, we can force them to be
binary and add new variables $t_{a}^{k}=T_{a}x_{a}^{k}=T_{a}\sum_{p\in
P^{k}}\delta_{a}^{p}z^{k}_{p}$ with the following constraints:

$\displaystyle 0\leq t_{a}^{k}\leq M_{a}^{k}x_{a}^{k}$ $\displaystyle a\in
A_{1},$ (19) $\displaystyle 0\leq T_{a}-t_{a}^{k}\leq N_{a}(1-x_{a}^{k})$
$\displaystyle a\in A_{1}$ (20)

for the primal-arc representation and

$\displaystyle 0\leq t_{a}^{k}\leq M_{a}^{k}\sum_{p\in
P^{k}}\delta_{a}^{p}z^{k}_{p}$ $\displaystyle a\in A_{1},$ (21) $\displaystyle
0\leq T_{a}-t_{a}^{k}\leq N_{a}\left(1-\sum_{p\in
P^{k}}\delta_{a}^{p}z^{k}_{p}\right)$ $\displaystyle a\in A_{1}$ (22)

for the primal-path representation. $M_{a}^{k}$ and $N_{a}$ are big-M
parameters. We refer this method of linearization as _direct linearization_.
Tight bounds for $M_{a}^{k}$ and $N_{a}$ can be found in Dewez et al. [7]. The
leader revenue for each commodity becomes $\sum_{a\in
A_{1}}\eta^{k}t_{a}^{k}$, which we can use to replace the bilinear terms. The
strong duality constraints of the four linearized formulations are as follows:

$\displaystyle\sum_{a\in A_{1}}(c_{a}x_{a}^{k}+t_{a}^{k})+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}$ $\displaystyle=\lambda^{k}_{o^{k}}-\lambda^{k}_{d^{k}},$
(23) $\displaystyle\sum_{a\in A_{1}}(c_{a}x_{a}^{k}+t_{a}^{k})+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}$ $\displaystyle=L^{k},$ (24) $\displaystyle\sum_{p\in
P^{k}}\sum_{a\in A}\delta_{a}^{p}c_{a}z^{k}_{p}+\sum_{a\in A_{1}}t_{a}^{k}$
$\displaystyle=\lambda^{k}_{o^{k}}-\lambda^{k}_{d^{k}},$ (25)
$\displaystyle\sum_{p\in P^{k}}\sum_{a\in
A}\delta_{a}^{p}c_{a}z^{k}_{p}+\sum_{a\in A_{1}}t_{a}^{k}$
$\displaystyle=L^{k}.$ (26)

The objective function of the linearized formulation is independent from the
primal representation:

$\displaystyle\sum_{k\in K}\sum_{a\in A_{1}}\eta^{k}t_{a}^{k}.$ (27)

A summary of all linearized strong duality formulations is provided in Table
3. The variables with stars are required to be binary.

|  |  | Constraints
---|---|---|---
Label | Variables | Obj. | Primal | Dual | Opt. Cond. | Linearization
STD | $T,t,x^{*},y,\lambda$ | (27) | (1) | (2), (3) | (23) | (19), (20)
VF | $T,t,x^{*},y,L$ | (27) | (1) | (5) | (24) | (19), (20)
PASTD | $T,t,z^{*},\lambda$ | (27) | (4) | (2), (3) | (25) | (21), (22)
PVF | $T,t,z^{*},L$ | (27) | (4) | (5) | (26) | (21), (22)
Table 3: List of variables, constraints, and objective functions of all
linearized strong duality formulations.

####  2.3.2 Linearization of Complementary Slackness Formulations

In complementary slackness formulations, bilinear terms appear in all
complementary slackness constraints. In the arc-arc (CS) and path-path (PCS)
combinations, each complementary slackness constraint is gated by a primal
variable ($x,y$, or $z$), which can only take binary values. Thus, we could
use them as branching condition for the constraints. Here are the linearized
constraints of the arc-arc combination (CS):

$\displaystyle\lambda^{k}_{i}-\lambda^{k}_{j}$ $\displaystyle\geq
c_{a}+T_{a}-R_{a}^{k}(1-x_{a}^{k})$ $\displaystyle a\equiv(i,j)\in A_{1},$
(28) $\displaystyle\lambda^{k}_{i}-\lambda^{k}_{j}$ $\displaystyle\geq
c_{a}-R_{a}^{k}(1-y_{a}^{k})$ $\displaystyle a\equiv(i,j)\in A_{2},$ (29)

and of the path-path combination (PCS):

$\displaystyle L^{k}$ $\displaystyle\geq\sum_{a\in
A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}-S^{k}_{p}(1-z^{k}_{p})$ $\displaystyle p\in P^{k}.$
(30)

$R_{a}^{k}$ and $S^{k}_{p}$ are big-M parameters depending on the dual
representation ($R_{a}^{k}$ for dual-arc and $S^{k}_{p}$ for dual-path). When
$x,y$ or $z$ are equal to 1, the big-M terms are dropped. Combining with the
dual constraints, these complementary slackness constraints force them to be
active. If the primal variables are equal to 0, the big-M terms make the
constraints redundant, essentially removing them from the formulation.

To calculate the values of these big-M parameters, recalling that $N_{a}$ is
the upper bound of $T_{a}$ defined in (20), denote:

  * •

$\underline{\lambda}^{k}_{i}$: the minimum cost from $i$ to $d^{k}$ when
$T_{a}=0,\forall a\in A_{1}$;

  * •

$\overline{\lambda}^{k}_{j}$: the minimum cost from $j$ to $d^{k}$ when
$T_{a}=N_{a},\forall a\in A_{1}$;

  * •

$\underline{L}^{k}$: the minimum cost from ${o^{k}}$ to ${d^{k}}$ when
$T_{a}=0,\forall a\in A_{1}$.

Then, the upper bounds for the big-M parameters are:

$\displaystyle R_{a}^{k}$
$\displaystyle=c_{a}+N_{a}-\underline{\lambda}^{k}_{i}+\overline{\lambda}^{k}_{j}$
$\displaystyle a\equiv(i,j)\in A_{1},$ $\displaystyle R_{a}^{k}$
$\displaystyle=c_{a}-\underline{\lambda}^{k}_{i}+\overline{\lambda}^{k}_{j}$
$\displaystyle a\equiv(i,j)\in A_{2},$ $\displaystyle S^{k}_{p}$
$\displaystyle=\sum_{a\in A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}N_{a}-\underline{L}^{k}$ $\displaystyle p\in P^{k}.$

In the calculation of $\overline{\lambda}^{k}_{j}$, excluding all the tolled
arcs (set $T=\infty$) is a valid option, nevertheless, it may disconnect $j$
and $d^{k}$ which may lead to infinite $R^{K}_{a}$. Thus, since we aim for the
smallest $R^{k}_{a}$, we do not remove them. A similar technique can be used
for the path-arc combination (PACS):

$\displaystyle\lambda^{k}_{i}-\lambda^{k}_{j}$ $\displaystyle\geq
c_{a}+T_{a}-R_{a}^{k}\left(1-\sum_{p\in P^{k}}\delta_{a}^{p}z^{k}_{p}\right)$
$\displaystyle a\equiv(i,j)\in A_{1},$ (31)
$\displaystyle\lambda^{k}_{i}-\lambda^{k}_{j}$ $\displaystyle\geq
c_{a}-R_{a}^{k}\left(1-\sum_{p\in P^{k}}\delta_{a}^{p}z^{k}_{p}\right)$
$\displaystyle a\equiv(i,j)\in A_{2}.$ (32)

For the arc-path combination (VFCS), the constraints have multilinear terms.
We will use the expression

$\sum_{a\in A}\delta_{a}^{p}-\sum_{a\in
A_{1}}\delta_{a}^{p}x_{a}^{k}-\sum_{a\in A_{2}}\delta_{a}^{p}y_{a}^{k}$

as the branching condition. The number of arcs in path $p$ is $\sum_{a\in
A}\delta_{a}^{p}$. This expression is only equal to 0 if all the arcs along
the path $p$ are active, _i.e._ $p$ is chosen. The linearized constraints of
the arc-path combination (VFCS) are:

$\displaystyle L^{k}$ $\displaystyle\geq\sum_{a\in
A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}-S^{k}_{p}\left(\sum_{a\in
A}\delta_{a}^{p}-\sum_{a\in A_{1}}\delta_{a}^{p}x_{a}^{k}-\sum_{a\in
A_{2}}\delta_{a}^{p}y_{a}^{k}\right)$ $\displaystyle p\in P^{k}.$ (33)

After linearizing all the bilinear terms in the complementary slackness
constraints, there are still bilinear terms in the objective function (for the
leader’s revenue). We could use the direct linearization method as in the
strong duality case which adds more big-M constraints to the formulations.
However, we can take advantage of the special structure of the NPP: the
leader’s revenue appears twice, once in the objective function, and once in
the strong duality constraint. Furthermore, since strong duality constraints
are not utilized in complementary slackness formulations, the leader’s revenue
can be extracted from these constraints and then be substituted in the
objective function. We will refer to this method as _linearization by
substitution_. Let $\tau^{k}=\sum_{a\in A_{1}}T_{a}x_{a}^{k}=\sum_{a\in
A_{1}}\sum_{p\in P^{k}}\delta_{a}^{p}T_{a}z^{k}_{p}$ be the leader’s revenue
for commodity $k$ per unit of demand. The strong duality constraints are used
to extract $\tau^{k}$:

$\displaystyle\sum_{a\in A_{1}}c_{a}x_{a}^{k}+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}+\tau^{k}$
$\displaystyle=\lambda^{k}_{o^{k}}-\lambda^{k}_{d^{k}},$ (34)
$\displaystyle\sum_{a\in A_{1}}c_{a}x_{a}^{k}+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}+\tau^{k}$ $\displaystyle=L^{k},$ (35)
$\displaystyle\sum_{p\in P^{k}}\sum_{a\in
A}\delta_{a}^{p}c_{a}z^{k}_{p}+\tau^{k}$
$\displaystyle=\lambda^{k}_{o^{k}}-\lambda^{k}_{d^{k}},$ (36)
$\displaystyle\sum_{p\in P^{k}}\sum_{a\in
A}\delta_{a}^{p}c_{a}z^{k}_{p}+\tau^{k}$ $\displaystyle=L^{k}.$ (37)

The objective function becomes:

$\displaystyle\sum_{k\in K}\eta^{k}\tau^{k}.$ (38)

Table 4 is the summary of all linearized complementary slackness formulations.
The variables with stars are required to be binary. Since there are two
linearization methods for the bilinear terms in the objective function, a
suffix is added to the label, where 1 means direct linearization and 2 means
linearization by substitution. Overall, the linearization of the complementary
slackness formulations is more complicated than that of the strong duality
formulations. Complementary slackness formulations usually have more binary
variables and more constraints, which may lead to worse performance.

|  |  | Constraints
---|---|---|---
Label | Variables | Obj. | Primal | Dual | Opt. Cond. | Linearization
Direct linearization |  |  |
CS1 | $T,t,x^{*},y^{*},\lambda$ | (27) | (1) | (2), (3) | (28), (29) | (19), (20)
VFCS1 | $T,t,x^{*},y^{*},L$ | (27) | (1) | (5) | (33) | (19), (20)
PACS1 | $T,t,z^{*},\lambda$ | (27) | (4) | (2), (3) | (31), (32) | (21), (22)
PCS1 | $T,t,z^{*},L$ | (27) | (4) | (5) | (30) | (21), (22)
Linearization by substitution |  |  |
CS2 | $T,\tau,x^{*},y^{*},\lambda$ | (38) | (1) | (2), (3) | (28), (29) | (34)
VFCS2 | $T,\tau,x^{*},y^{*},L$ | (38) | (1) | (5) | (33) | (35)
PACS2 | $T,\tau,z^{*},\lambda$ | (38) | (4) | (2), (3) | (31), (32) | (36)
PCS2 | $T,\tau,z^{*},L$ | (38) | (4) | (5) | (30) | (37)
Table 4: List of variables, constraints, and objective functions of all
linearized complementary slackness formulations.

Combining with the four strong duality formulations, in total, there are 12
different MILP formulations for the NPP. Some of them are completely new:
(PASTD), (CS), and (VFCS), while others are already mentioned in other works.
Table 5 shows how these reformulations fit in the big picture. In the left and
the middle columns, the names and the labels used in the original works are
listed. The corresponding labels used in this paper are shown in the right
column. We remark that in Bouhtou et al. [1], the path-based formulation
(PMIP) is a path value function formulation (PVF) with a minor modification:
the bilinear terms are linearized separately per path by replacing
$t_{a}^{k}=T_{a}x_{a}^{k}$ with $r^{k}_{pa}=T_{a}x_{a}^{k}z^{k}_{p}$.

Name in original work | Original label | General label
---|---|---
Brotcorne et al. [5] |  |
Mixed integer formulation | CPLEX | STD
Heilporn et al. [11] |  |
Mixed integer formulation | TP2 | STD
New formulation | TP3 | VF
Didi-Biha et al. [8] |  |
Arc-based formulation | MIP I | STD
Arc-path formulation | MIP II | PACS1
Path-based formulation | MIP III | PCS2
Bouhtou et al. [1] |  |
Arc-based formulation | AMIP | STD
Path-based formulation | PMIP | PVF*
Dewez et al. [7] |  |
Arc-based formulation | TOP-ARCS | STD
Path-based formulation | PATH | PCS2
Table 5: The mapping of previous formulations to the general framework.

##  3 Path Enumeration and Preprocessing

This section describes a new path enumeration process in details, combining
all the ingredients previously explored in the literature. Then, we discuss
the use of this process as a preprocessing method for arc-based formulations.
The latter will enable us to fairly compare arc and path-based reformulations,
since we make preprocessing available for both.

###  3.1 Path Enumeration

Formulations using path representations require the set $P^{k}$, and thus, the
explicit enumeration of all paths from ${o^{k}}$ to ${d^{k}}$. The size of
$P^{k}$ can be exponential. However, not all paths are relevant. We can
eliminate many of them by using the dominance rule in Bouhtou et al. [1].

######  Definition 1.

Given a commodity $k$, a path is bilevel feasible if it is optimal in the
follower problem for some value of $T$.

Mathematically, if $p$ is bilevel feasible, then there exists $T$ such that
for all $q\in P^{k}$:

$\sum_{a\in A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}\leq\sum_{a\in A}\delta_{a}^{q}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{q}T_{a}.$

######  Lemma 2 (Bouhtou et al. [1]).

Consider any commodity $k\in K$. Let $p$ and $q$ be two different paths in
$P^{k}$. If for all $a\in A_{1}$, $\delta_{a}^{p}\leq\delta_{a}^{q}$ and

$\sum_{a\in A}\delta_{a}^{p}c_{a}<\sum_{a\in A}\delta_{a}^{q}c_{a},$

then $q$ is not bilevel feasible, _i.e._ , $q$ cannot be the optimal path for
any value of $T$.

###### Proof.

Suppose that $q$ is the optimal path for some fixed value of $T$. The
follower’s cost of $q$ must be minimal:

$\sum_{a\in A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}\geq\sum_{a\in A}\delta_{a}^{q}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{q}T_{a}.$

However, because $\delta_{a}^{p}\leq\delta_{a}^{q}$, this means $\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}\leq\sum_{a\in A_{1}}\delta_{a}^{q}T_{a}$; and
because $\sum_{a\in A}\delta_{a}^{p}c_{a}<\sum_{a\in A}\delta_{a}^{q}c_{a}$,
this leads to a contradiction. ∎

The dominance rule in Lemma 2 is the only necessary rule to eliminate all non-
bilevel-feasible paths. Any remaining path is bilevel feasible.

######  Theorem 3.

Any path which is not eliminated by the dominance rule in Lemma 2 is bilevel
feasible.

###### Proof.

We will prove that if path $p$ is not eliminated, then it is optimal for the
following value of $T$:

$T_{a}=\begin{cases}0&\text{if }\delta_{a}^{p}=1,\\\
\infty&\text{otherwise}.\end{cases}$

Suppose that $p$ is not optimal for the above value of $T$ and let $q$ be the
optimal path. Consequently, $q$ can only use arcs with $\delta_{a}^{p}=1$
because all other tolled arcs are disabled. In other words,
$\delta_{a}^{q}\leq\delta_{a}^{p}$. Also because $T_{a}=0$ for
$\delta_{a}^{p}=1$, the costs of $p$ and $q$ in this case are exactly their
initial costs, hence $\sum_{a\in A}\delta_{a}^{q}c_{a}<\sum_{a\in
A}\delta_{a}^{p}c_{a}$. By Lemma 2, $p$ should have been eliminated which is a
contradiction. Therefore, $p$ must be optimal for the given value of $T$. ∎

Definition 1 alone does not ensure that any bilevel feasible path is relevant
to solve the NPP. Suppose we have two different paths with the same initial
cost and the same set of tolled arcs. Then they can be both bilevel feasible.
In such case, however, only one path is required. Thus, we employ Assumption 4
to make the term “bilevel feasible” coincide with the term “relevant”.

######  Assumption 4.

Two different bilevel feasible paths must have different initial costs.

This assumption also implies that two different bilevel feasible paths will
have different sets of tolled arcs. This can be fulfilled if we add a
minuscule random perturbation to each arc.

Yen’s algorithm [18] is usually used to enumerate the set of paths of a
commodity. The algorithm outputs the shortest path, the second shortest path,
etc, up to the $K$-th shortest path between two nodes in the graph. In this
case, the shortest path is the path with the minimum initial cost (the cost
when set $T=0$). The algorithm should be stopped when it has found the first
toll-free path as in the following corollary:

######  Corollary 5.

Given a commodity $k\in K$, a path $p$ cannot be bilevel feasible if its
initial cost is greater than that of the shortest toll-free path $\pi^{k}$.

In Bouhtou et al. [1], the authors reduced the size of the graph by
transforming the original graph to the shortest path graph model (SPGM). In
SPGM, all nodes that are not incident to any tolled arc will be removed and
the toll-free arcs connecting to it are replaced by the outer product of the
list of incoming arcs and the list of outgoing arcs. Then, some elimination
rules (which are only applicable on the SPGM) are applied to reduce the number
of arcs. Finally, the set of paths is enumerated using the Yen’s algorithm on
the SPGM and the dominance rule is applied afterward. A drawback of the SPGM
is the inflation of the set of arcs which is caused by the replacement of
nodes by arcs. Besides that, the elimination rules in SPGM still leave many
redundant nodes and arcs. Thus, there is a sizable number of irrelevant paths
generated by the Yen’s algorithm.

In Didi-Biha et al. [8], the authors used a modified version of the Lawler’s
procedure [14], a general form of the Yen’s algorithm, to enumerate the set of
bilevel feasible paths directly without using the SPGM. SPGM is not needed
here because it is only a method to reduce the size of the graph before Yen’s
algorithm is applied which is not used by the authors. Due to Assumption 4,
only paths with different combinations of tolled arcs are explored, ignoring
most redundant paths. However, the algorithm in [8] requires solving a
shortest path problem with some arcs fixed to 1 which cannot be accomplished
by the Dijsktra’s algorithm and requires a linear program solver. We propose
an improved version of this algorithm which just needs a shortest path
algorithm. Algorithm 1 shows this enumeration process in details. In the
algorithm, $C$ is the set of candidate paths, $N$ is the maximum number of
paths we want to enumerate, $p^{(j)}$ is the $j$-th shortest path, $s(q^{i})$
is the spur node of path $q^{i}$, and $R(q^{i})$ is the set of excluded tolled
arcs of path $q^{i}$.

Algorithm 1 Path enumeration

Input: The graph $G=(V,A)$, the O-D pair $o^{k},d^{k}$, the initial costs
$c_{a}$, the maximum number of paths to be enumerated $N$.

Output: A list of $N$ paths, most of which are bilevel feasible.

Find the first shortest path $p^{(1)}$

$C\leftarrow\\{p^{(1)}\\}$

$s(p^{(1)})\leftarrow{o^{k}}$

$R(p^{(1)})\leftarrow\varnothing$

$j\leftarrow 1$

while $j\leq N$ do

Output the path with minimum cost in $C$, call it $p^{(j)}$

Remove $p^{(j)}$ from $C$

if $p^{(j)}$ has no tolled arcs then

Stop the algorithm

end if

Let $a_{1},a_{2},\ldots,a_{m}$ be the ordered tolled arcs of $p^{(j)}$ from
$s(p^{(j)})$ to ${d^{k}}$

$\hat{s}\leftarrow s(p^{(j)})$

for $i$ from 1 to $m$ do

$q^{i}\leftarrow$ subpath of $p^{(j)}$ from ${o^{k}}$ to $\hat{s}$

$s(q^{i})\leftarrow\hat{s}$

$R(q^{i})\leftarrow R(p^{(j)})\cup\\{a_{i}\\}$

Remove nodes in $q^{i}$ (for this loop only)

Remove arcs in $R(q^{i})$ (for this loop only)

if $\hat{s}$ and ${d^{k}}$ are connected then

Append to $q^{i}$ the shortest path from $\hat{s}$ to ${d^{k}}$

Add $q^{i}$ to $C$

end if

$\hat{s}\leftarrow$ target of $a_{i}$

end for

$j\leftarrow j+1$

end while

Algorithm 1 is a version of the Lawler’s algorithm [14] which is an
enumeration scheme for binary problems such as the shortest path problem. Let
$S(p^{(j)})$ be the set of tolled arcs of $p^{(j)}$ from ${o^{k}}$ to
$s(p^{(j)})$. We can see that $p^{(j)}$ is a solution (not necessarily
optimal) of the shortest path problem when all arcs in $R(p^{(j)})$ are fixed
to 0 and all arcs in $S(p^{(j)})$ are fixed to 1. We will call $p^{(j)}$ a
solution of the problem $(R(p^{(j)}),S(p^{(j)}))$. The Lawler’s algorithm
tells us that we need to spawn $n$ subproblems where $n$ is the number of non-
fixed variables. We can re-order these non-fixed variables in any order.
Hence, the variables with value 1 ($x_{a_{1}}$ to $x_{a_{m}}$) come first and
the variables with value 0 ($x_{a_{m+1}}$ to $x_{a_{n}}$) follow. The former
$m$ variables are the tolled arcs of $p^{(j)}$ from $s(p^{(j)})$ to ${d^{k}}$
(all the tolled arcs preceding $s(p^{(j)})$ are fixed by assumption). These
$m$ variables must also be in order of appearance as mentioned in Algorithm 1.
According to this order, the $n$ subproblems of the Lawler’s algorithm are
generated by fixing all the previously fixed variables with these additional
constraints (one constraint for each subproblem):

$\displaystyle(1)$ $\displaystyle x_{a_{1}}=0$ $\displaystyle(2)$
$\displaystyle x_{a_{1}}=1,x_{a_{2}}=0$
$\displaystyle\vdots\leavevmode\nobreak\ $ $\displaystyle\hskip
28.45274pt\vdots$ $\displaystyle(m)$ $\displaystyle
x_{a_{1}}=x_{a_{2}}=\ldots=x_{a_{m-1}}=1,x_{a_{m}}=0$ $\displaystyle(m+1)$
$\displaystyle x_{a_{1}}=\ldots=x_{a_{m}}=1;x_{a_{m+1}}=1$
$\displaystyle(m+2)$ $\displaystyle
x_{a_{1}}=\ldots=x_{a_{m}}=1;x_{a_{m+1}}=0,x_{a_{m+2}}=1$
$\displaystyle\vdots\leavevmode\nobreak\ $ $\displaystyle\hskip
28.45274pt\vdots$ $\displaystyle(n)$ $\displaystyle
x_{a_{1}}=\ldots=x_{a_{m}}=1;x_{a_{m+1}}=\ldots=x_{a_{n-1}}=0,x_{a_{n}}=1$

We can verify that the feasible regions of all subproblems are mutually
exclusive and their union is equivalent to the feasible region of the problem
$(R(p^{(j)}),S(p^{(j)}))$ except $p^{(j)}$. The subproblems in the second
group have a common constraint $x_{a_{1}}=\ldots=x_{a_{m}}=1$. However,
because of Lemma 2, $x_{a_{1}}=\ldots=x_{a_{m}}=1$ implies that the resulting
paths cannot be bilevel feasible (dominated by $p^{(j)}$). Therefore, we only
need to consider the first $m$ subproblems.

Consider the subproblem $i$. Let $q^{i}$ be the path constructed as in
Algorithm 1. We call $q^{i}$ a child of $p^{(j)}$. There are two cases: we can
construct $q^{i}$ and we cannot (there is no path from $\hat{s}$ to
${d^{k}}$). If we can and $q^{i}$ exists, then similar to $p^{(j)}$, $q^{i}$
is a solution of the problem $(R(q^{i}),S(q^{i}))$ where
$R(q^{i})=R(p^{(j)})\cup\\{a_{i}\\}$ and
$S(q^{i})=S(p^{(j)})\cup\\{a_{1},\ldots,a_{i-1}\\}$. Although we do not assume
that $p^{(j)}$ and $q^{i}$ are optimal in their own subproblems, they satisfy
two special properties. First, the costs of $p^{(j)}$ and its child $q^{i}$
follow the correct enumeration order:

######  Lemma 6.

The cost of $q^{i}$ is at least the cost of $p^{(j)}$.

Hereafter, we will use the term _part_ to refer to a subpath of $p^{(j)}$ or
$q^{i}$ of which one end is either the origin ${o^{k}}$ or the destination
${d^{k}}$. The term _segment_ will refer to any other subpath.

###### Proof.

Because $p^{(j)}$ is constructed by Algorithm 1, its part from $s(p^{(j)})$ to
${d^{k}}$ is the shortest path excluding arcs in $R(p^{(j)})$ and all nodes
preceding $s(p^{(j)})$. As a result, its part from $s(q^{i})$ to ${d^{k}}$ is
also the shortest path with respect to the same set of constraints. Now,
consider the path $q^{i}$. Its part from ${o^{k}}$ to $s(q^{i})$ is identical
to that of $p^{(j)}$. The other part from $s(q^{i})$ to ${d^{k}}$ is the
shortest path, while excluding arcs in $R(q^{i})$ and all nodes preceding
$s(q^{i})$. Given that $R(p^{(j)})\subset R(q^{i})$ and $s(p^{(j)})$ comes
before $s(q^{i})$, the new shortest path problem is a restriction. Thus, the
cost of this part of $q^{i}$ must be at least the cost of the same part of
$p^{(j)}$. ∎

Lemma 6 is fundamental for using Corollary 5 as the stopping condition. If the
enumeration order is not maintained, then we may stop prematurely after
encountering the first toll-free path. The other property is expressed in the
lemma below and it will be used for proving the correction of another pruning
condition.

######  Lemma 7.

Any toll-free segment of $q^{i}$ (or any path generated by the algorithm) is
the shortest segment while excluding all arcs in $R(q^{i})$ and all the
preceding nodes of that segment.

###### Proof.

Assume that $q^{i}$ is a child of $p^{(j)}$ and the lemma is true for
$p^{(j)}$. The part of $p^{(j)}$ from ${o^{k}}$ to $s(q^{i})$ is reused for
$q^{i}$ and does not contain any arc in $R(q^{i})$, so the lemma is true for
this part. Consider the other part of $q^{i}$ from $s(q^{i})$ to ${d^{k}}$. By
construction in Algorithm 1, this part is the shortest path, while excluding
arcs in $R(q^{i})$ and nodes preceding $s(q^{i})$. As a result, any subpath of
this part is also the shortest subpath with respect to the same set of
constraints. Later, for any toll-free segment in this part, we can extent the
set of excluded nodes from the set of nodes that precedes $s(q^{i})$ to the
set of nodes that precedes that segment. This means that the lemma is also
true for this part of $q^{i}$. Recall that $s(q^{i})$ is the target of some
tolled arc, thus there are no toll-free segments containing $s(q^{i})$ in the
middle and the lemma is true for $q^{i}$. To finish the induction proof, we
consider the base case of the first shortest path for which the lemma is also
true. ∎

Now, we consider the other case where we cannot construct $q^{i}$. This means
$s(q^{i})$ and ${d^{k}}$ are disconnected when we remove arcs in $R(q^{i})$
and nodes preceding $s(q^{i})$. We will prove that there is no bilevel
feasible path in $(R(q^{i}),S(q^{i}))$:

######  Lemma 8.

If $q^{i}$ does not exists, any path in $(R(q^{i}),S(q^{i}))$ is not bilevel
feasible.

###### Proof.

If the problem $(R(q^{i}),S(q^{i}))$ is infeasible, then the above statement
is trivially true. Suppose it has a feasible solution $\hat{q}^{i}$ which is a
simple path (path with no loops). We will construct another path (not
necessarily a solution of $(R(q^{i}),S(q^{i}))$) which dominates $\hat{q}^{i}$
by replacing a segment in $\hat{q}^{i}$ by a shorter segment in $p^{(j)}$.
First, we need to sort the arcs in $S(q^{i})$ by the order of appearance in
$p^{(j)}$:

$\hat{a}_{1},\hat{a}_{2},\hat{a}_{3},\ldots,\hat{a}_{n-1},\hat{a}_{n}$

where $n=|S(q^{i})|$. We also define $\hat{a}_{0}$ as a virtual arc whose
target is ${o^{k}}$ (its source is not relevant). Next, for $0\leq i<n$, we
will try to replace the segment of $\hat{q}^{i}$ between $\hat{a}_{i}$ and
$\hat{a}_{i+1}$ with the same segment of $p^{(j)}$ if $\hat{a}_{i+1}$ comes
after $\hat{a}_{i}$ in $\hat{q}^{i}$ and if the latter segment has smaller
cost. If we are able to make such a replacement, due to the definition of
$S(q^{i})$, the segment connecting $\hat{a}_{i}$ and $\hat{a}_{i+1}$ in
$p^{(j)}$ is toll-free; hence, we only remove tolled arcs while building a
better path, and thus the resulting path dominates $\hat{q}^{i}$. If we cannot
replace any segment, this means that either all segments in $\hat{q}^{i}$ are
identical to those in $p^{(j)}$ or there is a segment in $\hat{q}^{i}$ with
smaller cost.

Consider the first case: all segments in $\hat{q}^{i}$ are identical to those
in $p^{(j)}$. In this case, the parts of $\hat{q}^{i}$ and $p^{(j)}$ from
${o^{k}}$ to $s(q^{i})$ are identical. The other part of $\hat{q}^{i}$ must be
a path from $s(q^{i})$ to ${d^{k}}$. Because we assume that $\hat{q}^{i}$ is a
simple path, we must exclude all nodes preceding $s(q^{i})$. However, since we
also assume that $q^{i}$ does not exist, this implies $s(q^{i})$ and ${d^{k}}$
to be disconnected if we exclude those nodes. This is a contradiction.

Consider the second case: there is a segment in $\hat{q}^{i}$ with smaller
cost. Let $\hat{a}_{m}$ be the first arc such that the segment of
$\hat{q}^{i}$ from $\hat{a}_{m}$ to $\hat{a}_{m+1}$ is cheaper than its
counterpart of $p^{(j)}$. All segments preceding $\hat{a}_{m}$ in both
$\hat{q}^{i}$ and $p^{(j)}$ must be identical since we assume that we cannot
replace any segment. Consequently, the part from ${o^{k}}$ to $\hat{a}_{m}$ of
both paths are the same. By Lemma 7, the segment of $p^{(j)}$ from
$\hat{a}_{m}$ to $\hat{a}_{m+1}$ is the optimal path while excluding all
preceding nodes. However, the same segment of $\hat{q}^{i}$ is better, which
means it must violate that constraint and must repeat some node preceding
$\hat{a}_{m}$. This is again a contradiction as we assumed that $\hat{q}^{i}$
is a simple path.

Since both cases result in a contradiction, we can always replace some segment
in $\hat{q}^{i}$ with a better segment in $p^{(j)}$. Therefore, $\hat{q}^{i}$
cannot be bilevel feasible. ∎

######  Theorem 9.

Given a sufficiently large value of $N$, the output of Algorithm 1 includes
all bilevel feasible paths.

###### Proof.

At any step of the algorithm, the subproblems generated by the Lawler’s
procedure do not remove any part of the feasible region of the original
problem. The subproblems are mutually exclusive and there is a finite number
of paths, thus the algorithm is guaranteed to stop in finite time. A
subproblem $i$ is pruned if and only if it is in the second group of
subproblems or if $q^{i}$ does not exists. Both cases imply that all solutions
in the feasible region of the subproblem are not bilevel feasible (Lemma 8).
Besides that, the stopping condition in Corollary 5 is ensured by Lemma 6.
Therefore, no bilevel feasible solution is ruled out. ∎

Although the set of paths enumerated in Algorithm 1 is close to the final set
of bilevel feasible paths, there are still redundant paths. To obtain the
final set of bilevel feasible paths, the dominance rule in Lemma 2 still needs
to be applied once more.

######  Example 10.

Consider the graph in Figure 1. The number of each arc represents the initial
cost $c_{a}$. The first shortest path is ${o^{k}}-u-v-{d^{k}}$ with the cost
of 3. By Algorithm 1, three subproblems are generated, each producing a
candidate path. They are the second shortest path ${o^{k}}-u-{d^{k}}$ (cost
4), the third shortest path ${o^{k}}-u-v-w-{d^{k}}$ (cost 6), and the toll-
free path ${o^{k}}-{d^{k}}$ (cost 10). All 4 paths will be returned, however,
the path ${o^{k}}-u-{d^{k}}$ dominates ${o^{k}}-u-v-w-{d^{k}}$. The final set
of bilevel feasible paths only has 3 paths: ${o^{k}}-u-v-{d^{k}}$,
${o^{k}}-u-{d^{k}}$, and ${o^{k}}-{d^{k}}$.

${o^{k}}$$u$$v$$w$${d^{k}}$11122103 Figure 1: Graph with redundant paths
(dashed arcs are tolled arcs).

###  3.2 Path-based Preprocessing

In [8, 1], the authors compare their proposed path-based formulations to an
unprocessed or SPGM-processed standard formulation. In these comparisons,
path-based formulations are more advantageous compared to their counterparts.
Although SPGM can reduce the complexity of a graph to some degree, path-based
formulations essentially bypass all the redundant arcs and nodes in the graph,
hence they have more preprocessing power. However, we could harness this
preprocessing power of path enumeration and apply it to arc-based
formulations. We propose a new preprocessing method which uses the set of
enumerated paths to remove all redundant arcs and nodes completely. The idea
is simple: given the set of bilevel feasible paths, any arc or node which does
not appear in any path will be removed. Then, to further reduce the size of
the graph, chains of toll-free arcs will be replaced with a single toll-free
arc with cost equal to the sum of the costs of the associated arcs. This
preprocessing method will let the arc-based formulations use the information
generated from the path enumeration process and enable us to compare all the
formulations in a fair manner.

In addition, SPGM can still be applied on top of the processed graph. This
does not remove any additional tolled arcs because the set of tolled arcs
obtained from path-based preprocessing is already minimal. However, SPGM
transforms the graph by replacing nodes with toll-free arcs, which can
sometimes reduce both the number of nodes and the number of toll-free arcs in
sparse graphs.

As mentioned in Section 2.3, the conversions used in the (PACS) and the (VFCS)
formulations must cover all solutions in the primal representation. When we
apply the path-based preprocessing to the arc representation and use the set
of bilevel feasible paths for the path representation, the (PACS) formulation
still satisfies this rule. This is because the set of primal solutions of
(PACS) is the set of bilevel feasible paths, all of them can be converted to
arcs in the processed graph. On the other hand, the (VFCS) formulation does
not satisfy this rule anymore. In the processed graph, there are many paths
which are not bilevel feasible. If these paths are not covered in the (VFCS)
formulation, the optimality condition can be bypassed by selecting these paths
as the primal solutions. However, enumerating all paths (including non-
bilevel-feasible) in the processed graph would be too costly. Therefore, a
special treatment is neccesary and we propose a cutting-plane method just for
(VFCS): the processed graph and the set of bilevel feasible paths will still
be used to generate the initial set of constraints. Then, when we encounter a
solution, we will examine if the path in that solution is a bilevel feasible
path. If it is not, we will use it to generate an additional complementary
slackness constraint to cover the new path and continue solving.

##  4 Hybrid Model for Multi-Commodity Problems

Although path-based formulations and preprocessing can give the MIP solver a
boost by removing redundant variables, they require time to enumerate the
paths. Since the number of paths can be exponential, sometimes, it is faster
to just solve the problem without relying on path enumeration (unprocessed or
SPGM-processed STD or CS). However, it is unknown which option is better until
the enumeration process is finished. Instead of a full enumeration, we could
enumerate until a predetermined number of paths and then decide if it is worth
to continue enumerating or to cancel the enumeration process and fallback to a
formulation which does not require path enumeration. Since we have a multi-
commodity problem, this probing method can be applied for each commodity
separately. The end result is a hybrid model, assigning different formulations
and preprocessing for different commodities depending on their numbers of
paths. This process is shown in Algorithm 2. In the algorithm, $N$ is called
the breakpoint, which is a threshold to decide whether the enumeration process
should be continued. The commodities with only 1 path are redundant, because
the only path must be a toll-free path which does not bring any profit to the
leader. Commodities with less than $N$ paths are assigned to a path-based
formulation or an arc-based formulation with path preprocessing, while
commodities with more than $N$ paths are assigned to a fallback formulation
which does not need path enumeration. All variables are managed independently
according to their corresponding commodities, except for $T_{a}$ which is
shared.

Algorithm 2 Hybrid model

Input: The graph $G=(V,A)$ with costs $c_{a}$, all commodities
$(\eta^{k},{o^{k}},{d^{k}})$.

Output: A hybrid model $(obj,constraints)$.

$obj\leftarrow 0$

$constraints\leftarrow\varnothing$

for $k\in K$ do

Enumerate the first $N+1$ paths of commodity $k$

Let $\hat{P}^{k}$ be the set of enumerated paths

if $|\hat{P}^{k}|>1$ then

if $|\hat{P}^{k}|\leq N$ then

Assign a formulation $F^{k}$ which is either path-based or arc-based with
path-based preprocessing

else

Assign an arc-based formulation $F^{k}$ without path-based preprocessing

end if

$obj\leftarrow obj+obj(F^{k})$

$constraints\leftarrow constraints\cup constraints(F^{k})$

end if

end for

Use a solver to solve $(obj,constraints)$

As an example, suppose we divide the set of commodities $K$ into 3 parts based
on the number of paths: $K_{1}$ for commodities having only 1 path, $K_{2}$
for commodities with less than $N$ paths, and $K_{3}$ for commodities with
more than $N$ paths. If we assign (PCS2) to $K_{2}$ and (STD) to $K_{3}$, then
we get the following hybrid model:

$\displaystyle\max\ $ $\displaystyle\sum_{k\in
K_{2}}\eta^{k}\tau^{k}+\sum_{k\in K_{3}}\sum_{a\in A_{1}}\eta^{k}t_{a}^{k}$
s.t. $\displaystyle\eqref{pp},\eqref{dp}$ $\displaystyle k\in K_{2},$
$\displaystyle L^{k}\geq\sum_{a\in A}\delta_{a}^{p}c_{a}+\sum_{a\in
A_{1}}\delta_{a}^{p}T_{a}-S^{k}_{p}(1-z^{k}_{p})$ $\displaystyle k\in
K_{2},p\in P^{k},$ $\displaystyle\sum_{p\in P^{k}}\sum_{a\in
A}\delta_{a}^{p}c_{a}z^{k}_{p}+\tau^{k}=L^{k}$ $\displaystyle k\in K_{2},$
$\displaystyle z^{k}_{p}\in\\{0,1\\}$ $\displaystyle k\in K_{2},p\in P^{k},$
$\displaystyle\eqref{pa},\eqref{da1},\eqref{da2},\eqref{directa1},\eqref{directa2}$
$\displaystyle k\in K_{3},$ $\displaystyle\sum_{a\in
A_{1}}(c_{a}x_{a}^{k}+t_{a}^{k})+\sum_{a\in
A_{2}}c_{a}y_{a}^{k}=\lambda^{k}_{o^{k}}-\lambda^{k}_{d^{k}}$ $\displaystyle
k\in K_{3},$ $\displaystyle x_{a}^{k}\in\\{0,1\\}$ $\displaystyle k\in
K_{3},a\in A_{1}.$

The hybrid model can be extended to have multiple breakpoints and different
assignment schemes. For example, path-based formulations are less effective
for commodities which have many paths because the more paths a commodity has,
the more variables and constraints are added. In contrast, arc-based
formulations with path preprocessing only remove arcs and nodes which
guarantee complexity reduction. We could use the hybrid model to balance
between these two kinds of formulations, by assigning path-based formulations
to commodities with few paths and arc-based formulations for those with many
paths.

##  5 Experiments

We conducted computational experiments to evaluate in practice the
reformulations and preprocessing presented in the previous sections. In
Section 5.1, the experimental methodology is described, including the
implementation details, the generation and properties of the instances used.
Section 5.2 provides a comparison between the 12 reformulations listed in
Section 2. Section 5.3 validates the efficiency of the new path-based
preprocessing. In Section 5.4, we justify the use of the hybrid model
introduced in Section 4.

###  5.1 Methodology

The instances used in the tests is randomly generated. We employ a generation
method similar to the one described in [4] to make the problems challenging.
For a given graph and a set of O-D pairs, first, the shortest path of each
commodity is found. Next, we count the number of paths passing through each
arc and sort them in the descending order of that count. Following that order,
each arc is converted into a tolled arc until 2/3 of the desired number of
tolled arcs is reached. The last 1/3 is selected randomly among all remaining
arcs. An arc is converted only if it does not remove the last toll-free path
for all commodities. To make the generated data more realistic, we also
enforce the properties of the arc to be symmetrical, which means that any arc
and its reversed arc will have the same cost, and both must be either tolled
or toll-free. The cost of 80% of all arcs will be distributed uniformly from 5
to 35, while the remaining 20% will have the maximum cost of 35. The cost of
tolled arcs are halved after the conversion. The proportion of tolled arcs is
20%.

There are 200 generated problem instances (or in short, problems) divided into
four sets which are summarized in Table 6. Each set has 50 problems with 5
different numbers of commodities: 30, 35, 40, 45, and 50 (10 problems
generated for each value). The topologies are inspired from [6] and are
illustrated in Figure 2. The set G provides a dataset similar to the one used
in [4] and [5], while the three other sets have more paths and are more
challenging to solve. The cumulative distribution of the number of bilevel
feasible paths of each set is shown in Figure 3. The horizontal axis
represents the breakpoint $N$ while the vertical axis shows the proportion of
commodities with no more than $N$ bilevel feasible paths.

Label | Topology | Dimensions | $|V|$ | Avg. $|A|$
---|---|---|---|---
G | Grid |  $5\times 12$ nodes | 60 | 206
H | Grid |  $12\times 12$ nodes | 144 | 528
D | Delaunay | - | 144 | 832
V | Voronoi | - | 144 | 410
Table 6: Properties of generated data. GridDelaunayVoronoi Figure 2:
Illustration of the three types of topology.
$\displaystyle{10^{0}}$$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{10^{4}}$Breakpoint
(N)$\displaystyle{0.0}$$\displaystyle{0.2}$$\displaystyle{0.4}$$\displaystyle{0.6}$$\displaystyle{0.8}$$\displaystyle{1.0}$Proportion
of commoditiesGHDV Figure 3: Cumulative distribution of the number of bilevel
feasible paths.

Our code is implemented in C++. All test runs are executed single-threaded on
the Béluga cluster (Intel Xeon 2.4 GHz) under Linux. All tested models are
instantiations of the hybrid model described in Algorithm 2. Given a problem,
the path enumeration will be applied first, then a formulation will be
assigned to each commodity and the single-level reformulation is synthesized.
Finally, the single-level reformulation is solved directly by CPLEX 12.9. If a
model does not rely on path enumeration, the breakpoint $N$ is set to 1 and
path enumeration is still applied. This will remove the commodities with only
1 path and provide basic preprocessing for all models. The main formulation
will be assigned for the commodities with less than $N$ bilevel feasible
paths. The default fallback formulation for commodities with more than $N$
paths is the (STD) with no preprocessing. By increasing $N$, we essentially
replace the fallback formulation with the main formulation. If the main
formulation is better, then we should see an improvement in performance when
increasing $N$ and vice versa. Some formulations are only good up to a certain
breakpoint, which can also be observed if there is a change in the direction
of their performance over $N$.

We imposed a time limit of 1 hour for each problem. The time for path
enumeration is included in the total time of a test run. If the run exceeds 1
hour, it is stopped and the optimality gap is recorded. The problems are
divided into two groups: easy and hard. The easy group consists of 108
problems which can be solved by at least one run. All the remaining 92
problems form the hard group, for which an optimal solution is not available.
When we compare the performance of any two models, we are interested in three
criteria: the number of problems solved by each model, the average time it
takes to solve the easy group, and the average optimality gap of the the hard
group.

###  5.2 Comparison of all Formulations

In this section, we compare the performance of all 12 formulations. These will
be assigned as the main formulations of the hybrid models. If the main
formulation is arc-based, then path-based preprocessing is applied. If the
main formulation is a mixed arc-path or path-arc formulation, then we use
path-based preprocessing on the arc representation. If the main formulation is
a path-based formulation, no further preprocessing is applied. In all cases,
path enumeration is required for the main formulations.

Table 7 shows the performance of all models over $N$. At first look, we can
see that formulations with suffix 1 worsen over $N$. This means they are even
worse than the fallback formulation. These are the formulations with
complementary slackness as optimality condition and direct linearization. The
reason may be because both complementary slackness and direct linearization
use big-M constraints. All other formulations perform roughly the same (except
for (VFCS2)), with the standard formulation taking the lead in all three
criteria.

Figure 4 plots the performance over $N$ of four notable formulations: (STD),
(VF), (CS2), and (PCS2). (STD), (VF), and (CS2) represent three basic
paradigms to convert a bilevel linear problem to a single-level reformulation.
(PCS2) is a path-based formulation which is mentioned in [8] and [7]. In these
papers, their experiments suggest that (PCS2) outperforms (STD), but this
comparison is done when the (STD) formulation is preprocessed by the SPGM
method. Here, we observe the opposite result: the standard formulation
outperformed (PCS2). The key is the new path-based preprocessing, which
theoretically provides arc-based formulations with the same preprocessing
power of the path enumeration, thus it levels the playing field. Figure 4 also
shows a drawback of path-based formulations: the larger $N$ is, the more
complex they become. (PCS2) must add an extra binary variable for every
enumerated path, while (STD) and (CS2) use the set of paths to eliminate arcs
and nodes. Thus, the complexity of (PCS2) increases over $N$ while the
complexity of the other two formulations is capped by its unprocessed version.
In the graph, we can see that the performance of (PCS2) starts dropping after
$N=500$. The other formulations using primal-path representation ((PASTD),
(PVF), (PACS)) also suffer from this problem. (VF) uses dual-path, so it only
adds more constraints instead of binary variables which allows it to continue
improving over $N$.

$N$ | 10 | 20 | 50 | 100 | 200 | 500 | 1000
---|---|---|---|---|---|---|---
std | 61 | 64 | 69 | 72 | 78 | 92 | 94
vf | 59 | 64 | 68 | 68 | 75 | 90 | 93
pastd | 56 | 64 | 68 | 71 | 76 | 86 | 85
pvf | 59 | 64 | 67 | 69 | 79 | 80 | 81
cs1 | 47 | 32 | 8 | 1 | 0 | 0 | 0
cs2 | 58 | 62 | 67 | 67 | 70 | 71 | 73
pacs1 | 49 | 33 | 8 | 1 | 0 | 0 | 0
pacs2 | 58 | 62 | 63 | 68 | 70 | 76 | 69
vfcs1 | 43 | 10 | 0 | 0 | 0 | 0 | 0
vfcs2 | 56 | 55 | 44 | 36 | 29 | 23 | 21
pcs1 | 48 | 33 | 9 | 2 | 0 | 0 | 0
pcs2 | 61 | 61 | 72 | 71 | 77 | 80 | 79
(a) Number of problems solved

$N$ | 10 | 20 | 50 | 100 | 200 | 500 | 1000
---|---|---|---|---|---|---|---
std | 1999 | 1864 | 1739 | 1604 | 1430 | 1127 | 1029
vf | 2022 | 1870 | 1730 | 1644 | 1520 | 1231 | 1155
pastd | 2056 | 1863 | 1724 | 1662 | 1498 | 1315 | 1220
pvf | 1980 | 1867 | 1770 | 1699 | 1486 | 1409 | 1429
cs1 | 2396 | 2915 | 3495 | 3585 | 3601 | 3601 | 3601
cs2 | 2069 | 1973 | 1787 | 1822 | 1697 | 1599 | 1518
pacs1 | 2566 | 3368 | 3601 | 3601 | 3601 | 3601 | 3601
pacs2 | 2135 | 2189 | 2438 | 2687 | 2890 | 2986 | 3050
vfcs1 | 2359 | 2818 | 3491 | 3587 | 3601 | 3601 | 3601
vfcs2 | 1992 | 1916 | 1852 | 1712 | 1673 | 1572 | 1676
pcs1 | 2351 | 2861 | 3448 | 3576 | 3601 | 3601 | 3601
pcs2 | 2012 | 1903 | 1696 | 1599 | 1544 | 1529 | 1607
(b) Time of the easy group

$N$ | 10 | 20 | 50 | 100 | 200 | 500 | 1000
---|---|---|---|---|---|---|---
std | 12.8 | 12.7 | 11.9 | 11.3 | 10.5 | 9.3 | 8.5
vf | 13.3 | 12.6 | 12.1 | 11.4 | 11.2 | 10.0 | 9.8
pastd | 13.2 | 12.8 | 11.7 | 11.7 | 10.8 | 9.7 | 9.2
pvf | 13.6 | 12.7 | 12.1 | 11.7 | 11.1 | 10.1 | 10.6
cs1 | 16.3 | 21.8 | 42.1 | 70.0 | 94.2 | 99.6 | 100.0
cs2 | 13.3 | 12.9 | 12.8 | 12.2 | 11.9 | 10.9 | 10.9
pacs1 | 18.9 | 34.4 | 70.3 | 92.6 | 99.1 | 100.0 | 100.0
pacs2 | 13.7 | 13.4 | 14.6 | 15.1 | 15.5 | 17.1 | 19.2
vfcs1 | 16.3 | 20.3 | 36.3 | 58.2 | 86.2 | 97.7 | 99.5
vfcs2 | 13.6 | 13.5 | 13.5 | 13.3 | 12.8 | 12.1 | 12.1
pcs1 | 15.6 | 20.3 | 39.5 | 63.2 | 88.4 | 98.7 | 99.9
pcs2 | 13.2 | 13.3 | 12.4 | 11.8 | 10.8 | 10.6 | 11.1
(c) Gap of the hard group

Table 7: Performance of all formulations over $N$.
$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{60}$$\displaystyle{70}$$\displaystyle{80}$$\displaystyle{90}$N.o.
problems
solvedstdvfcs2pcs2$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{1000}$$\displaystyle{1200}$$\displaystyle{1400}$$\displaystyle{1600}$$\displaystyle{1800}$$\displaystyle{2000}$Time
of Easy group
(s)stdvfcs2pcs2$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$Breakpoint
(N)$\displaystyle{9}$$\displaystyle{10}$$\displaystyle{11}$$\displaystyle{12}$$\displaystyle{13}$Gap
of Hard group (%)stdvfcs2pcs2 Figure 4: Performance of notable formulations
over $N$.

###  5.3 Path-based Preprocessing

In this section, we compare the new path-based preprocessing with the state-
of-the-art preprocessing method (SPGM). Figure 5 shows the reduction ratio
between those two preprocessing methods. Given a breakpoint $N$ (horizontal
axis), the reduction ratio is the quotient of the total number of
nodes/arcs/tolled arcs of all commodities with no more than $N$ bilevel
feasible paths after being processed by the path-based preprocessing over the
same sum in the original graph (absolute ratio) or in the SPGM-processed graph
(relative ratio). Each criterion (node, arc, tolled arc) has a different
impact on the final single-level reformulation, which is summarized in Table
8. The tolled arc criterion is the most important because it is linked with
the number of binary variables in the final formulation. In Figure 5, we can
see that the path-based preprocessing excels in all three criteria. On
average, graphs processed by the path-based preprocessing have 10% less nodes,
66% less arcs, and 49% less tolled arcs compare to the SPGM method. Compared
to the original graph, the path-based preprocessing removes 75% of all tolled
arcs, which reduces the size of the problems significantly. The efficiency
varies between different instance sets, with the path-based preprocessing
being more favorable in more connected graphs such as in dataset D and H.

$\displaystyle{10^{0}}$$\displaystyle{10^{2}}$$\displaystyle{10^{4}}$$\displaystyle{0.0}$$\displaystyle{0.1}$$\displaystyle{0.2}$$\displaystyle{0.3}$NodeAbsoluteGHDVAvg$\displaystyle{10^{0}}$$\displaystyle{10^{2}}$$\displaystyle{10^{4}}$$\displaystyle{0.6}$$\displaystyle{0.8}$$\displaystyle{1.0}$$\displaystyle{1.2}$Relative
to
SPGMGHDVAvg$\displaystyle{10^{0}}$$\displaystyle{10^{2}}$$\displaystyle{10^{4}}$$\displaystyle{0.00}$$\displaystyle{0.05}$$\displaystyle{0.10}$$\displaystyle{0.15}$ArcGHDVAvg$\displaystyle{10^{0}}$$\displaystyle{10^{2}}$$\displaystyle{10^{4}}$$\displaystyle{0.4}$$\displaystyle{0.6}$$\displaystyle{0.8}$$\displaystyle{1.0}$GHDVAvg$\displaystyle{10^{0}}$$\displaystyle{10^{2}}$$\displaystyle{10^{4}}$$\displaystyle{0.0}$$\displaystyle{0.1}$$\displaystyle{0.2}$$\displaystyle{0.3}$Tolled
ArcGHDVAvg$\displaystyle{10^{0}}$$\displaystyle{10^{2}}$$\displaystyle{10^{4}}$$\displaystyle{0.4}$$\displaystyle{0.5}$$\displaystyle{0.6}$GHDVAvg

Figure 5: Reduction ratio of Path-based Preprocessing. Criteria | Impact
---|---
Node |  • Num. of constraints in primal-arc
|  • Num. of variables in dual-arc
Arc |  • Num. of variables in primal-arc (strong duality)
|  • Num. of binary variables in primal-arc
| (complementary slackness)
|  • Num. of constraints in dual-arc
|  • Num. of constraints in complementary slackness
Tolled arc |  • Num. of binary variables in primal-arc
|  • Num. of variables and constraints in direct linearization
Table 8: Impact of the reduction ratio criteria.

Figure 6 shows the real performance of the two preprocessing methods over $N$.
The main formulation in both models is the (STD), but one is matched with
path-based preprocessing, while the other is matched with SPGM preprocessing.
Once again, we can observe that the path-based preprocessing excels in all
three criteria. It can solve 30 more problems (over 200 in total), in half the
time and half the gap compared to the SPGM model. Besides that, the SPGM model
hardly improves when $N>1000$ while path-based preprocessing keeps thriving
beyond this threshold. In our experiments, applying SPGM on top of the path-
based preprocessing does not produce significantly better results, and
sometimes can even be detrimental. This concludes that path-based
preprocessing is the superior preprocessing method.

$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{10^{4}}$$\displaystyle{60}$$\displaystyle{70}$$\displaystyle{80}$$\displaystyle{90}$$\displaystyle{100}$N.o.
problems solvedPath-
basedSPGM$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{10^{4}}$$\displaystyle{0}$$\displaystyle{500}$$\displaystyle{1000}$$\displaystyle{1500}$$\displaystyle{2000}$Time
of Easy group (s)Path-
basedSPGM$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{10^{4}}$Breakpoint
(N)$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$Gap
of Hard group (%)Path-basedSPGM Figure 6: Performance comparison between Path-
based preprocessing and SPGM preprocessing.

###  5.4 Hybrid Model

In this work, we show numerical results to justify the idea of hybrid model
introduced in Section 4 as a way to compromise between the time spent for path
enumeration and the time spent to solve the reformulation. We will use the
best reformulation which is the (STD) for all commodities in this test. Figure
7 shows its performance with $N$ varying from 10 to 100000 (horizontal axis).
If a commodity has no more than $N$ bilevel feasible paths, then path-based
preprocessing is applied. Otherwise, the unprocessed graph is used instead.
The larger $N$ is, the longer it takes to enumerate, but more commodities will
be processed. If $N$ is very large, all commodities will be processed and all
paths are enumerated. If $N$ is set to 0, then it is identical to an
unprocessed model. From Figure 7, we observe that larger $N$ generally leads
to better performance until a certain point around 10000 paths. After that,
the time spent for path enumeration does not have a positive impact on the
performance anymore. The average time of the easy group is not affected by
large $N$ because they do not have many paths to enumerate in the first place.
At $N=100000$, there are two unsolved problems in class D which require more
than 1 hour for path enumeration and this contributes to the rise at the end
in the average optimality gap of the hard group. In conclusion, for the given
instances, stopping the path enumeration at $N=10000$ provides a good balance.
Stopping early also improves the robustness of the program, because although
commodities with more than 10000 paths are not common, encountering one could
make the program stuck for a long time.

$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{10^{4}}$$\displaystyle{10^{5}}$$\displaystyle{60}$$\displaystyle{70}$$\displaystyle{80}$$\displaystyle{90}$$\displaystyle{100}$N.o.
problems
solved$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{10^{4}}$$\displaystyle{10^{5}}$$\displaystyle{750}$$\displaystyle{1000}$$\displaystyle{1250}$$\displaystyle{1500}$$\displaystyle{1750}$$\displaystyle{2000}$Time
of Easy group
(s)$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{10^{4}}$$\displaystyle{10^{5}}$Breakpoint
(N)$\displaystyle{8}$$\displaystyle{10}$$\displaystyle{12}$Gap of Hard group
(%) Figure 7: Performance of the standard model over $N$.

##  6 Conclusion

In this work, we related many modeling concepts which enrich and unify the
general toolbox to solve bilevel problems. If the follower problem can be
formulated in multiple ways, then we can mix the primal and the dual (if it
exists) of any two formulations to write the single-level formulation. In the
case of a (mixed-integer) linear follower problem, this alternative
formulation always exists: the follower problem can be written as a linear
combination of the extreme points of the convex hull of the feasible set.
Bilevel feasibility, and its enumeration, can be used as a preprocessing
method. The hybrid framework is applicable in any multi-commodity problem with
multiple reformulations. The network pricing problem is a good example to
demonstrate all these modeling techniques together.

The NPP makes some theoretical assumptions whose relaxation enables a more
close mirroring of reality. From a practical point of view, future work must
consider additional real-life features, including arc construction cost,
capacity on arcs and nodes, congestion, competition, and uncertainty, many of
which are highly non-linear and make the problem significantly harder.

## Acknowledgements

The authors are grateful for the support of Institut de valorisation des
données (IVADO) and Fonds de recherche du Québec (FRQ) through the FRQ-IVADO
Research Chair, and of the Natural Sciences and Engineering Council of Canada
(NSERC) through its Discovery Grant program.

This research was enabled in part by support provided by Calcul Québec
(www.calculquebec.ca) and Compute Canada (www.computecanada.ca).

## References

  * [1] Mustapha Bouhtou, Stan Hoesel, Anton Kraaij, and Jean-Luc Lutton.  Tariff optimization in networks.  INFORMS Journal on Computing, 19, 08 2007.
  * [2] Luce Brotcorne, F Cirinei, Patrice Marcotte, and Gilles Savard.  A tabu search algorithm for the network pricing problem.  Computers & Operations Research, 39(11):2603–2611, 2012.
  * [3] Luce Brotcorne, Fabien Cirinei, Patrice Marcotte, and Gilles Savard.  An exact algorithm for the network pricing problem.  Discrete Optimization, 8(2):246–258, 2011.
  * [4] Luce Brotcorne, Martine Labbé, Patrice Marcotte, and Gilles Savard.  A bilevel model and solution algorithm for a freight tariff-setting problem.  Transportation Science, 34(3):289–302, 2000.
  * [5] Luce Brotcorne, Martine Labbé, Patrice Marcotte, and Gilles Savard.  A bilevel model for toll optimization on a multicommodity transportation network.  Transportation science, 35(4):345–358, 2001.
  * [6] Luce Brotcorne, Martine Labbé, Patrice Marcotte, and Gilles Savard.  Joint design and pricing on a network.  Operations research, 56(5):1104–1115, 2008.
  * [7] Sophie Dewez, Martine Labbé, Patrice Marcotte, and Gilles Savard.  New formulations and valid inequalities for a bilevel pricing problem.  Operations Research Letters, 36(2):141 – 149, 2008.
  * [8] Mohamed Didi-Biha, Patrice Marcotte, and Gilles Savard.  Path-based formulations of a bilevel toll setting problem, pages 29–50.  Springer US, Boston, MA, 2006.
  * [9] Matteo Fischetti, Ivana Ljubić, Michele Monaci, and Markus Sinnl.  A new general-purpose algorithm for mixed-integer bilevel linear programs.  Operations Research, 65(6):1615–1637, 2017.
  * [10] François Gilbert, Patrice Marcotte, and Gilles Savard.  A numerical study of the logit network pricing problem.  Transportation Science, 49(3):706–719, 2015.
  * [11] G Heilporn, M Labbé, P Marcotte, and G Savard.  New formulations and valid inequalities for the toll setting problem.  IFAC Proceedings Volumes, 39(3):431–436, 2006.
  * [12] Géraldine Heilporn, Martine Labbé, Patrice Marcotte, and Gilles Savard.  A polyhedral study of the network pricing problem with connected toll arcs.  Networks: An International Journal, 55(3):234–246, 2010.
  * [13] Martine Labbé, Patrice Marcotte, and Gilles Savard.  A bilevel model of taxation and its application to optimal highway pricing.  Management science, 44(12-part-1):1608–1622, 1998.
  * [14] Eugene L Lawler.  A procedure for computing the k best solutions to discrete optimization problems and its application to the shortest path problem.  Management science, 18(7):401–405, 1972.
  * [15] Sébastien Roch, Gilles Savard, and Patrice Marcotte.  An approximation algorithm for stackelberg network pricing.  Networks: An International Journal, 46(1):57–67, 2005.
  * [16] Sahar Tahernejad, Ted K Ralphs, and Scott T DeNegre.  A branch-and-cut algorithm for mixed integer bilevel linear optimization problems and its implementation.  Mathematical Programming Computation, pages 1–40, 2020.
  * [17] CPM van Hoesel, Anton F van der Kraaij, Carlo Mannino, Gianpaolo Oriolo, and Mustapha Bouhtou.  Polynomial cases of the tarification problem.  METEOR, Maastricht University School of Business and Economics, 2003.
  * [18] Jin Y Yen.  Finding the k shortest loopless paths in a network.  Management Science, 17(11):712–716, 1971.