Ons

Roua

Amel

Aymen, Jalal Isaias George

Keywords: Combinatorial optimization, Column generation, Branch-and-Price, Valid inequalities, Benders decomposition, Branch-and-Benders-cut, Heuristic, Virtual Network Functions, Software-Dened Networking, Service Functions Chaining. vi Optimisation combinatoire, Génération de colonnes, Branch-and-Price, inégalités valides, décomposition de Benders, Branch-and-Benders-Cut, Heuristique, Fonctions virtualisées . . . . . . .

I would like to thank my husband Youcef for being a source of love, strength, inspiration and permanent support during the last three years,

Résumé

La Virtualisation des fonctions réseau et les réseaux dénis par logiciel sont deux nouvelles technologies prometteuses qui émergent dans la nouvelle génération des réseaux de télécommunication. Leur utilisation permet la minimisation du temps de traitement des services et un gain en énergie et en coûts. Dans cette thèse, nous étudions un problème auquel sont confrontés les fournisseurs des services réseau, intitulé le problème de placement des fonctions virtuelles et de routage des services dans les réseaux dénis par logiciel.

Étant donnés un ensemble de fonctions de réseau virtuelles (VNF) et un ensemble de paires de n÷uds origine-destination (commodités), le problème étudié dans cette thèse consiste à trouver les installations optimales des fonctions virtuelles sur les noeuds du réseau de télécommunication, et le chemin associé satisfaisant les contraintes de latence et passant par les fonctions virtuelles installées. Ce problème est un problème d'optimisation combinatoire dicile à résoudre. Évidemment, le problème devient encore plus dicile, si en plus, les ux de données doivent être acheminés en utilisant le concept de Chaînage de fonctions de service (SFC) pour lequel les VNFs doivent être traitées le long du chemin de routage dans un ordre prédéni. D'autres contraintes techniques peuvent être considérées dans ce problème. L'objectif est de minimiser les coûts d'installation des fonctions virtuelles sur les noeuds du réseau et les coûts d'activation des noeuds.

Dans ce manuscrit, nous commençons par prouver que le problème est NP-Dicile au sens fort, même en considérant une seule commodité, et en relâchant les contraintes de capacité sur les n÷uds et sur les fonctions virtuelles, et aussi les contraintes de latence et de précédence. Nous proposons ensuite une formulation PLNE (Programmation linéaire en nombres entiers) compacte pour modéliser le problème. Cette formulation ne semble pas être assez forte pour trouver des solutions réalisables en un temps raisonnable à l'aide d'un solveur standard. An de remédier à cela, nous fournissons, une heuristique basée sur une formulation PLNE.

Nous proposons également, deux formulations PLNE étendues pour modéliser le problème et nous dénissons l'algorithme de Branch-and-Price associé. Pour chacune de ces deux formulations, nous présentons, la procédure de génération de

Introduction

Typically, telecommunication networks are composed of dierent devices such as switches, routers, and middleboxes, which are used in order to provide network functions such as Proxies, Firewalls (FW), Intrusion Detection Systems (IDS), Load Balancers (LB), etc. The hardware middleboxes are costly, energy-intensive and have short lifecycles [START_REF] Juliver | Resource allocation in NFV: A comprehensive survey[END_REF]. Constantly, network service providers have to deal with a large number of trac requests in order to satisfy customer demands. These trac requests represent services such as video streaming, virtual private network (VPN), online gaming, etc.

Moreover, each service requires a specic Service Function Chain (SFC), representing a sequence of hardware devices that must be traversed by the data ows in a pre-dened order. Routing packets through these xed devices is very time consuming and reduces the QoS (e.g., end-to-end latency). Also, when a new service is required, new hardware devices should be used, placed, and connected to the network elements [START_REF] Li | Software-dened network function virtualization: A survey[END_REF].

Furthermore, in traditional telecommunication networks, no device has visibility of the whole network. Each switch has its own data plane, to forward the trac, and its own control plane, to decide where to send the data. Having a considerable number of devices implies managing a considerable number of control and data planes. This requires a manual conguration, and makes packet forwarding very complicated within a network [START_REF] Feamster | The road to SDN: an intellectual history of programmable networks[END_REF][START_REF] Pfa | Extending networking into the virtualization layer[END_REF]. Therefore, network management becomes very challenging [START_REF] Astuto | A survey of software-dened networking: Past, present, and future of programmable networks[END_REF],

which is very expensive in terms of CAPEX and OPEX.

During the last decade, the legacy telecommunication networks have started their transformation by introducing two promising technologies, namely, the Network Function Virtualization and Software-Dened Networking. The Network Function Virtualization allows the network functions to be executed as software on commercial o-theshelf equipment, thus being instantiated on-demand without installing new equipment.

Software-dened networking is a new paradigm that simplies network management and makes the deployment of new services easier by separating the data plane from the control plane. This is done by creating one or multiple central controllers having a global view of the whole network state.

Introduction

NFV and SDN provide various benets: specically, the VNFs could be installed ondemand on network nodes, and the routing paths associated with the Service Function

Chains could be computed dynamically to encounter the installed VNFs [START_REF] Allybokus | Virtual function placement for service chaining with partial orders and anti-anity rules[END_REF]. These two technologies bring new exibility in network management. Nevertheless, their introduction generates hard decision problems combining the choice of locations to deploy the VNFs and the computation of ecient routing paths meeting the installed VNFs in the right order. Both problems are widely studied separately in the literature, but their combination, specied by the VNF chaining problem, introduces a new challenge [6].

In this thesis, we study the Virtual Network Functions Placement and Routing Problem faced by Network Service Providers (NSPs). Given a telecommunication network, a set of VNFs and a set of trac requests (or commodities), the VNFPR problem consists in determining the optimal VNFs placement at network nodes, and the associated latency-constrained routing paths for each commodity, satisfying the SFC constraints.

The VNF-installation and the node activation costs are to be minimized. Further technical constraints, allowing a better Quality of the Service (QoS) and Quality of Experience (QoE), can be considered along with nodes capacities, VNFs capacities, or conict between VNFs of the same service.

In the rst part of the thesis, we investigate the problem's basic properties, we show that the problem is strongly NP-hard even for its simplest version, and we propose an MILP compact formulation to model it. This formulation does not seem to be strong enough to nd a solution using an o-the-shelf solver in a reasonable time. To tackle the problem from a computational perspective, we propose a path-based heuristic that provides optimal solutions for realistic instances from the literature.

Afterwards, we propose two extended formulations based on Dantzig-Wolfe decomposition to model the problem, namely: path formulation (PF) and Dantzig-Wolfe (DW) formulation. In order to strengthen their LP-bounds, we propose several families of valid inequalities and demonstrate their benets. We show, theoretically, that the Dantzig-Wolfe formulation is stronger than the path formulation in terms of linear relaxation. We present a branching scheme for each one and develop a Branch-and-Price algorithm. We compare both algorithms with the MILP compact formulation and the automatic Benders of Cplex.

In the last part of the dissertation, we study a variant of the problem in which we relax the VNF-capacity and conict constraints. We provide theoretical results that allow us to reformulate the problem using Benders decomposition and three families of valid inequalities to strengthen the LP-bounds. We combine all these ingredients in a Branch-and-Benders-Cut framework, and we test it on a set of realistic benchmark instances.

This manuscript is organized as follows. In Chapter 1, we briey present the basic notions of combinatorial optimization, its most important approaches, and dene some notions of the telecommunication network. This chapter also includes a review of the state-of-the-art methods for the VNFPRP. In Chapter 2, we present the VNFPRP, its proprieties, and an MILP compact formulation to model it. Chapter 3 is dedicated to the MILP-based heuristic. In Chapter 4, we present two extended formulations and their Branch-and-Price algorithms. In Chapter 5, we discuss a variant of the problem, and we provide theoretical results that allow us to derive Benders reformulation and to develop a Branch-and-Benders-Cut algorithm.

Chapter 1

Preliminaries and State-of-the-Art

This chapter introduces and presents the main notions and denitions necessary to ease the manuscript's understanding. First, we summarize the most important methods in combinatorial optimizaion, such as cutting plane, Branch-and-Cut, Benders reformulation, column generation, and the Branch-and-Price algorithms. Furthermore, the principal graph theory's notions and the key denitions in the telecommunication network are reported. The last section of the chapter is related to the review of the literature on the Virtual Network Functions Placement and Routing problem. Outline of the chapter. In Section 1.1, we present the most important reformulation and resolution approaches used in combinatorial optimization. Some important notions of graph theory are provided in Section 1.2. Section 1.3 is dedicated to the telecommunication network notions and denitions. The last section summarizes the literature review related to the VNFPRP.

1.1 Polyhedra and Integer Linear Programming methods Some denitions in this section have been collected from the book chapter of A.R Mahjoub [START_REF] Ridha | Polyhedral approaches[END_REF].

Elements of polyhedral theory

Let x ∈ R n be a vector, where n is a positive integer. We say that x is a linear combination of x 1 , x 2 , . . ., x m ∈ R n if there exist m scalars λ 1 , λ 2 , . . ., λ m such that x = m i=1 λ i x i . If m i=1 λ i = 1, then x is said to be a ane combination of x 1 , x 2 , . . ., x m . Moreover, if λ i ≥ 0, for all i ∈ {1, . . . , m}, we say that x is a convex combination of x 1 , x 2 , . . ., x m . Given a set S = {x 1 , . . . , x m } ∈ R n×m , the convex hull of S is the set of points x ∈ R n which are convex combination of x 1 , . . ., x m (see Figure 1.1), that is conv(S) = {x ∈ R n |x is a convex combination of x 1 , . . . , x m }.

The points x 1 , . . ., x m ∈ R n are linearly independent if the unique solution of the system m i=1 λ i x i = 0, is λ i = 0, for all i ∈ {1, . . . , m}. Furthermore, they are anely independent if the unique solution of the system m i=1 λ i x i = 0, m i=1 λ i = 1, is λ i = 0, i = 1, . . ., m.

A polyhedron P is the set of solutions of a linear system Ax ≤ b, that is P = {x ∈ R n |Ax ≤ b}, where A is a m-row n-columns matrix and b is a vector, b ∈ R m . A polytope is a bounded polyhedron. A point x of P will be also called a solution of P .

An inequality ax ≤ α is valid for a polyhedron P ⊆ R n if for every solution x ∈ P , ax ≤ α. This inequality is said to be tight for a solution x ∈ P , if ax = α. The inequality ax ≤ α is violated by x ∈ P , if ax > α. Let ax ≤ α be a valid inequality for the polyhedron P . F = {x ∈ P |ax = α} is called a face of P .

An inequality ax ≤ α is redundant, if the system {A x ≤ b } obtained by removing this inequality from Ax ≤ b denes the same polyhedron P . This is the case when ax ≤ α can be written as a linear combination of inequalities of the system A x ≤ b .

A solution is an extreme point of the polyhedron P , if and only if it cannot be written as the convex combination of two dierent solutions of P . The polyhedron P can also be described by its extreme points. In fact, every solution of P can be written as a convex combination of some extreme points of P . Figure 1.2 illustrates the main denitions given is this section.

Cutting plane method

Usually, the characterization of the convex hull of a combinatorial optimization problem is done by dening a large number of inequalities, which is exponential in most of the cases. Dene all these linear inequalities is not always obvious. In order to overcome this diculty, the cutting plane method can be used. The cutting plane method is based on the so-called separation problem. This consists, given a polyhedron P of R n and a point x * ∈ R n , in verifying if x * belongs to P or not, and if not, to identify an inequality a T x ≥ b, valid for P and violated by x * . In the latter case, we say that the hyper-plane a T x = b separates P and x * (see Figure 1.3). The cutting plane method consists in solving successive linear programs, with possibly a large number of inequalities, as follows. Let LP = min{cx | Ax ≥ b} be a Preliminaries and State-of-the-Art linear program and LP a linear program obtained by considering a small number of inequalities among Ax ≥ b. Let x * be the optimal solution of the latter system. We solve the separation problem associated with Ax ≥ b and x * . This phase is called the separation phase. If every inequality of Ax ≥ b is satised by x * , then x * is also optimal for LP . If not, let ax ≥ α be an inequality violated by x * . Then we add ax ≥ α to LP and repeat this process until an optimal solution is found. Algorithm 1 summarizes the dierent cutting plane steps.

Algorithm 1: A cutting plane algorithm [START_REF] Magnouche | The multi-terminal vertex separator problem : Complexity, Polyhedra and Algorithms[END_REF] Data: A linear program LP and its system of inequalities Ax ≥ b Result: Optimal solution x * of LP 1 : Consider a linear program LP with a small number of inequalities of LP ; 2 : Solve LP and let x * be an optimal solution; 3 : Solve the separation problem associated with Ax ≥ b and x * ; if an inequality ax ≥ α of LP is violated by x * then Add ax ≥ α to LP ; Go to step 2 ;

end else

x * is optimal for LP ; return x * ; end

The cutting-plane algorithm aims to solve the LP-relaxation of a given problem. The integrality of the provided solution has to be checked, if an integer solution is required and the Branch-and-Bound algorithm must be applied if the solution is fractional. Some inequalities have to be added to the problem during the Branch-and-Bound algorithm to achieve the optimality. This means that we have to combine both algorithms in order to get an optimal integer solution; the resulting algorithm is called Branch-and-Cut algorithm.

Branch-and-Cut algorithm

Consider a combinatorial optimization problem P which can be stated as: min{cx|Ax ≥ b, x ∈ {0, 1} n }, where Ax ≥ b has a large number of inequalities. A Branch-and-Cut algorithm starts by creating a Branch-and-Bound tree whose root node corresponds to a linear program LP 0 = min{cx|A 0 x ≥ b 0 , 0 ≤ x ≤ 1}, where A 0 x ≥ b 0 is a subsystem of Ax ≥ b having a small number of inequalities. Then, we solve the linear relaxation of the problem that is LP = {cx|Ax ≥ b, 0 ≤ x ≤ 1} using a cutting plane algorithm starting from LP 0 . Let x * 0 denote its optimal solution and A 0 x ≥ b 0 denote the set of inequalities added to LP 0 at the end of the cutting plane phase. If x * 0 is integral, then it is optimal. Otherwise, if x * 0 is fractional, we perform a branching phase. This step consists of choosing a variable, say x 1 , with a fractional value and adding two nodes P 1 and P 2 in the Branch-and-Cut tree. The node P 1 corresponds to the linear program LP 1 = min{cx|A 0 x ≥ b 0 , A 0 x ≥ b 0 , x 1 = 0} and LP 2 = min{cx|A 0 x ≥ b 0 , A 0 x ≥ b 0 , x 1 = 1}. We then solve the linear program LP 1 = min{cx|Ax ≥ b, x 1 = 0} (resp., LP 2 = min{cx|Ax ≥ b, x 1 = 1}) by a cutting plane method, starting from LP 1 (resp. LP 2 ). If the optimal solution of LP 1 (resp. LP 2 ) is integral then, it is feasible for the problem. Its value is then an upper bound of the optimal solution, and the node P 1 (resp. P 2 ) becomes a leaf of the Branch-and-Cut tree. If the solution is fractional, then we select a variable with a fractional value and add two children to the node P 1 (resp. P 2 ), and so on.

Similarly, as in the Branch-and-bound algorithm, infeasible nodes can be generated by adding constraints of type x i = 0 and x i = 1, with x i fractional variable. Also, this can generate nodes with a worst objective value than the best known lower bound.

Both types of nodes are pruned. The Branch-and-cut algorithm terminates when all nodes have been explored.

A good upper bound can be used in the Branch-and-Cut algorithm in order to prune nodes that do not allow its improvement. This allows reducing the size of the Branchand-Cut tree, and consequently, reduce the time used by the algorithm. Furthermore, this upper bound may be improved by applying a primal heuristic which aims to produce a feasible solution from the solution obtained at a given node of the Branchand-Cut tree when this later solution is fractional (and hence infeasible). When the solution computed is better than the best known upper bound, it may signicantly reduce the number of generated nodes, as well as the CPU time. Moreover, this guarantees to have an approximation of the optimal solution before visiting all the nodes of Branch-and-Cut tree, for example, when a CPU time limit has been reached.

Benders decomposition

Benders decomposition is one of the most famous decomposition tools for Mathematical Programming, proposed by JF Benders in [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF][START_REF] Benders | Partitioning in mathematical programming[END_REF], in order to reformulate and solve specic large-scale optimization problems admitting a set of complicating variables.

Complicating variables can be dened as those variables that, when temporarily xed, make the remaining optimization problem considerably more tractable. This means that xing their values reduces the given problem to an ordinary linear program [START_REF] Arthur M Georion | Generalized benders decomposition[END_REF].

Benders method is based on the cutting-plane method and utilizes the duality theory in order to derive families of cuts. Hence, instead of considering the problem with all its decision variables and constraints of a large-scale problem, Benders decomposition partitions the problem into multiple smaller problems called Benders subproblems and one master problem called Benders master problem. The computational diculty of optimization problems increases signicantly with a large number of variables and constraints. Solving these smaller problems iteratively can be more ecient than solving a single large problem. The Benders method is generalized for MINLP by Georion in [START_REF] Arthur M Georion | Generalized benders decomposition[END_REF].

Let (P ) be an MILP admitting two set of variables x and y. (P ) is dened as follows:

(P ) : min cx + dy

s.t Ax + By ≥ b x ∈ R p + y ∈ Y,
where Y is a complicated polyhedron, A and B are matrices, b, c and d are vectors with appropriate dimensions. Once complicated variables y are xed, the problem (P ) becomes signicantly easier to solve. The obtained problem represents a linear programming problem in x variables, which (depending on the structure of A) can be decomposed into smaller subproblems for subsets of x. Benders decomposition separates problem (P ) into two problems: (i) a master problem that contains the y variables, and (ii) subproblems that contain the x variables. We rst note that problem (P ) can be written with xed y variables as follows:

(P ) : min For the feasible region of (D) we have:

F = {u : u t A ≤ c t , u ≥ 0}
where :

F represents a polyhedron independent of y.

F is composed of extreme points u p , p = 1, . . . , P and extreme rays r q , q = 1, . . . , Q 1) In the case where (D) is unbounded this means that the problem (SP ) is infeasible. Therefore, there exists an unboundedness direction (extreme ray) r such that:

(b -B ȳ) t r > 0 ⇒ φ D (ȳ) → ∞.
In order to eliminate ȳ ∈ Y which causes the Benders subproblem infeasibility, a Benders feasibility cut is added to (P ):

(b -B ȳ) t r ≤ 0.

All infeasible y ∈ Y can be eliminated by adding constraints (b -By) t r q ≤ 0, q = 1, . . . , Q.

2) (D) is feasible for ȳ ∈ Y , by LP-duality theory we have:

φ(ȳ) = φ D (ȳ) = max u≥0 {(b -B ȳ) t u : A t u ≤ c}.
As each problem can be expressed by its extreme points and the optimality is attained at one of them, we have:

φ(ȳ) = φ D (ȳ) = max p∈{1,...,P } {(b -B ȳ) t u p }.
from which we can derive the Benders optimality cuts (see below). Then, (P ) can be reformulated by the following program:

(P ) : min y∈Y {dy + max p∈{1,...,P } (b -By) t u p } s.t (b -By) t r q ≤ 0, q = 1, . . . , Q

Preliminaries and State-of-the-Art

Using an auxiliry variables w to bound the inner maximization problem, the problem (P ) is equivalent to the following MIP: Since there is, typically, an exponential number of extreme points and extreme rays associated with the dual formulation, the Benders formulation admits an exponential number of constraints. Generating all constraints of type (OptCut) and (FeasCut) is very time-consuming and implies adding non-ecient constraints to the model. Instead, using the Benders algorithm, Benders cuts are separated at each iteration. We start the algorithm by solving the relaxed master problem with a subset of constraints. Let y * be the optimal solution of this problem. Then we solve the associated dual problem based on y * in order to generate Benders feasibility and optimality cuts. This procedure is repeated until no further notated Benders cuts are found. This method converges to an optimal solution in a nite number of iterations [START_REF] Arthur M Georion | Generalized benders decomposition[END_REF].

There exist two ways to apply the Benders method, old and modern one. In the old version, the relaxed master problem is solved as an MILP after adding cuts. Similarly to the Branch-and-Cut algorithm, the modern version (also called Branch-and-Benderscut ) consists in solving, at each node of the branching tree, the relaxed master problem using Benders method. Some successful applications using Benders method are: large-scale stochastic optimization [START_REF] Slyke | L-shaped linear programs with applications to optimal control and stochastic programming[END_REF], network design, xed-charge network design problems [START_REF] Alysson | A survey on Benders decomposition applied to xed-charge network design problems[END_REF], linear and quadratic facility location problems [START_REF] Fischetti | Benders decomposition without separability: A computational study for capacitated facility location problems[END_REF][START_REF] Fischetti | Redesigning Benders decomposition for large-scale facility location[END_REF] and covering location problems [START_REF] Cordeau | Benders decomposition for very large scale partial set covering and maximal covering location problems[END_REF].

Dantzig-Wolfe decomposition

The Dantzig-Wolfe decomposition denition have been collected from the paper of M.E. Lübbecke and J. Desrosiers [START_REF] Marco | Selected topics in column generation[END_REF]. Let P = {x ∈ R n + | Dx ≥ d} = ∅ be a polyhedron, {p q } q∈Q its extreme points and {p r } r∈R its extreme rays, with Q and R nite. Each x ∈ P can be written as convex combination of extreme points plus non-negative combination of extreme rays as follows:

In this subsection, we briey introduce the

x = q∈Q p q λ q + r∈R p r λ r , q∈Q λ q = 1, λ ∈ R |Q|+|R| +
By applying the linear transformation c j = c T p j and a j = Ap j , j ∈ Q ∪ R we obtain the following formulation:

(D -W ) : min q∈Q c q λ q + r∈R c r λ r (1.1) s.t. q∈Q a q λ q + r∈R a r λ r ≥ b (u) (1.2) q∈Q λ q = 1 (v) (1.3) λ ≥ 0 (1.4)
The (D-W) model admits a huge number of variables, but it possibly admits less constraints than (C). Constraints (1.3) are called convexity constraints. The compact formulation (C) and the extensive formulation (D-W) are equivalent and they give the same optimal objective function value. Nevertheless, they respective polyhera are not combinatorially equivalent [8,[START_REF] Marco | Selected topics in column generation[END_REF].

As the (D-W) formulation has a large number of variables, it should be solved using a column generation procedure (see Subsection 1.1.6). First, the model is initialized by a subset of variables, (it is called the restricted master problem, RMP), then the variables Preliminaries and State-of-the-Art needed to solve its linear relaxation are added at each iteration of the column generation procedure. Given a dual optimal solution (ū, v) of the RMP, the subproblem in Dantzig-Wolfe decomposition determines the column j which is optimal for min j∈Q∪R {c j -ūT a j -v}.

This objective function representing the minimum reduced cost c * , which is, using the linear transformation, equivalent to:

c * := min{(c T -ūT A)x -v|Dx ≥ d, x ≥ 0}
If c * ≥ 0, this means that no column with negative reduced cost exists and the column generation procedure terminates. If c * < 0, this means that the optimal solution of the pricing problem is either an extreme point or an extreme ray, and the respective column is added to the RMP.

Column generation procedure

Some denitions in this section, Section 1.1.7 and Section 1.2 have been token from the thesis of Y.Magnouche [START_REF] Magnouche | The multi-terminal vertex separator problem : Complexity, Polyhedra and Algorithms[END_REF].

Generally, compact MILP formulations performance is fairly limited because they often provide a weak linear relaxation. According to Sadykov and Vanderbeck [START_REF] Sadykov | Column generation for extended formulations[END_REF],

working in an extended variable space allows one to develop tight reformulations for mixed-integer programs. Those reformulations can admit an exponential number of variables which cannot be considered explicitly in the model. In this section, we introduce a method that allows for solving the LP-relaxation of extended formulations. This method is called Column generation method.

The column generation method is utilized to solve linear programs admitting an exponential number of variables. Only a small number of variables is considered at the beginning of the procedure. This method was pioneered by Dantzig and Wolfe in 1960 [START_REF] George | Decomposition principle for linear programs[END_REF] in order to solve problems that could not be managed eciently (in terms of CPU time and memory consumption) because of their size. Column generation is usually applied either for problems whose structure is suitable for a Dantzig-Wolfe decomposition (see Subsection 1.1.5), or for problems with a large number of variables.

Gilmore and Gomory [START_REF] Paul | A linear programming approach to the cutting-stock problem[END_REF][START_REF] Paul | A linear programming approach to the cutting stock problempart ii[END_REF] utilized this method to solve a cutting stock problem with a huge number of variables.

The column generation procedure aims to solve a sequence of linear programs with a restricted number of variables (also referred to as columns). The algorithm starts by solving a linear program having a small number of variables, and such that a feasible solution for the original problem may be identied using this basis. At each iteration of the algorithm, a so-called pricing problem is solved with the objective to identify the variables which must enter the current base. A negative reduced cost characterizes these variables. The reduced cost associated with a variable is computed using the dual variables associated with the problem's constraints. Then, the linear program that is obtained by adding the generated variables is solved; this procedure is repeated until no variable with negative reduced cost can be identied by the pricing problem. In this case, the solution obtained from the last restricted program is optimal for the original model. The main step of the column generation procedure is summarized in Algorithm 2.

Algorithm 2: A column generation algorithm [START_REF] Magnouche | The multi-terminal vertex separator problem : Complexity, Polyhedra and Algorithms[END_REF] Data : A linear program MP (Master Problem) with a huge number of variables Output : An optimal solution x * of MP 1: Consider a linear program RMP (Restricted Master Problem) including only a small subset of variables of the MP; 2: Solve RMP and let x * be its optimal solution; 3: Solve the pricing problem associated with the dual variables obtained by the resolution of the RMP; 4: If there exists a variable x with a negative reduced cost then; 5: add x to RMP.

6:

go to 2.

7: else 8:

x * is optimal for MP. 9:

return x * .
The column generation procedure can be seen as the dual of the cutting plane method since it adds columns (variables) instead of rows (inequalities) in the linear program.

Moreover, the pricing problem may be NP-hard. One can then use heuristic procedures to solve it. For more details on column generation algorithms, the reader is suggested to consult [START_REF] Desrosiers | A primer in column generation[END_REF][START_REF] Lübbecke | Selected Topics in Column Generation[END_REF][START_REF] Vanderbeck | Decomposition and column generation for integer programs[END_REF].

Branch-and-Price algorithm

The column generation procedure aims to solve the LP-relaxation of an optimization problem, having a huge number of variables. The integrality of the solution provided is not guaranteed. Hence to provide an optimal integer solution for an ILP formulation,

Preliminaries and State-of-the-Art the column generation method has to be merged with a Branch-and-bound algorithm; the resulting algorithm is called Branch-and-Price algorithm. This algorithm can be considered as the dual of the Branch-and-Cut algorithm.

The Branch-and-Price algorithm consists of solving each node of the Branch-and-Bound tree using the column generation procedure, all columns with negative reduced cost are added to the model in order to improve its LP-relaxation. When no column with negative reduced cost exists and the current solution is not integer, the algorithm starts the branching phase.

The Branch-and-Price algorithm has been used on various elds to solve large scale integer programming problems, and even real-life problems such as Cutting stock problem [START_REF] Alves | A stabilized branch-and-price-and-cut algorithm for the multiple length cutting stock problem[END_REF], Generalized Assignment Problem (GAP) [START_REF] Savelsbergh | A Branch-And-Price Algorithm for the Generalized Assignment Problem[END_REF], Airline Crew Scheduling [START_REF] Barnhart | Airline Crew Scheduling[END_REF],

Multi-commodity Flow Problems [START_REF] Barnhart | Using branch-and-price-and-cut to solve origin-destination integer multicommodity ow problems[END_REF], etc.

Graph theory

The current section is dedicated to introducing basic denitions in graph theory, which will be used throughout the chapters of this dissertation. For more details, we refer the reader to [START_REF] Schrijver | Combinatorial optimization: polyhedra and eciency[END_REF].

Undirected graphs

An undirected graph is denoted by G = (V, E) where V is the set of vertices or nodes and E is the set of edges. If e is an edge between two vertices u and v, then u and v are called the ends of e, and we write e = uv. If u is an end-point of e, then u (resp. e) is said to be incident to e (resp. u). Similarly, two vertices u and v forming an edge are said to be adjacent.

Let u and v be two vertices of V . A path p between u and v is an alternating sequence of vertices (v 0 , v 1 , v 2 , ..., v k-1 , v k ), where v 0 = u, v k = v; v i-1 v i represents an edge connecting v i-1 and v i in the path p for i = 1, . . . , k. p is called elementary if it does not visit more than once the same node in G.

Directed graphs

A directed graph is denoted by D = (N, A) where N is the set of nodes and A the set of arcs.

If a ∈ A is an arc connecting a node u to node v, then u will be called initial end and v is called nal end and we write a = (u, v). We say that a is an outgoing arc of node u and an incoming arc of node v. The vertices u and v are called ends of a. If v is an end (initial or nal) of a, then v (resp. a) is said to be incident to a (resp. v).

Shortest path algorithms

In this subsection, we present the well-known algorithms proposed for nding the shortest path between two nodes in a graph.

Dijkstra's algorithm

Given a graph G = (V, A) (resp. G = (V, E)), with positive arc (resp. edge) length l uv ∈ R + , for each (u, v) ∈ A (resp. uv ∈ E) the Dijkstra's algorithm aims to nd a shortest path to all nodes from a single source in a network [START_REF] Edsger | A note on two problems in connexion with graphs[END_REF], [START_REF] Donald | A note on dijkstra's shortest path algorithm[END_REF]. The original algorithm aims to nd the shortest paths from a dened source to all nodes in the graph (destinations). This algorithm can be adapted to provide the shortest path between a end end end return previous 1.2.3.2 Moore-Bellman-Ford's algorithm Given a graph G = (V, A) (resp. G = (V, E)), with arc (resp. edge) length l uv ∈ R, for each (u, v) ∈ A (resp. uv ∈ E) the Moore-Bellman-Ford's algorithm aims to nd a shortest path from a single source to all nodes in a network, without negative length

for i = 1 to |V | -1 do for each arc (v, w) ∈ A do if dist[w] > l uv + dist[v] then dist[w] := l uv + dist[v] predesessor[w] := v end end end return predesessor 1.2.

Yen's algorithm

Given a graph G = (V, A), Yen's algorithm aims to nd K-shortest loopless paths in the graph G between source and destination nodes [START_REF]Finding the k shortest loopless paths in a network[END_REF]. Yen's algorithm employs any shortest path algorithm (Dijkstra, for example) in order to nd the shortest path, then proceeds to nd K -1 deviations of this path, by deleting at each iteration one arc belonging to the shortest path. The complexity of Yen's algorithm based on Djikstra for nbP ath = 1 to K do for each node i in Sol[nbP ath -1] do This section is dedicated to introducing basic denitions and notations in the telecommunication networks, which will be used throughout the chapters of this dissertation.

First, we dene the network structure briey; we present some examples of network functions. Then, we describe the notion of Network Function Virtualization (NFV),

we dene the Software-Dened Networking (SDN) architecture and explain the Service Functions Chaining (SFC).

Network Structure

The telecommunication network is composed of three parts, as shown in Figure 1.4:

Access network is the part of a telecommunication network that gives the user access to services.

Core network (or backbone) is the part of a network that connects the dierent parts of the access network. The core network also provides the gateway to other networks.

Data center is a dedicated space within a building, or a group of buildings used to house computer systems and associated components, such as telecommunications and storage systems.

Preliminaries and State-of-the-Art

Network devices

Computer Networks consist of a large number of network devices such as routers, switches, and many types of middleboxes and links connecting them. In the following, we dene some network nodes.

Network interface is a software or a hardware interface connecting two parts of equipment. It has the ability to process low-level network information.

Repeater and hub is an electronic equipment that receives a network signal, cleans it of undesirable noise, and regenerates it. The signal is re-transmitted at a higher power level, or to the other side of obstruction so that the signal can traverse longer distances without degradation.

Bridge is a network device that connects and lters trac between two network segments at the data link layer 2 of the OSI [1] model to create a single network.

Switch: is a device that connects several segments (cables or bers) in a computer and telecommunication network, making it possible to create virtual circuits. Switching is one of two frame transport modes within computer and com- Router is a device that aims to determine routes and forwards packets between networks by processing the routing information included in the packets. They have a total view of the network.

Modem: (MOdulator-DEModulator) is a hardware device that transforms data into a suitable format for a transmission medium so that it can be transmitted from one computer to another.

Firewall is a network device for controlling network security and access rules. Typical examples include rewalls, video optimizers, load balancers, parental control, etc. Traditionally, SFs were implemented on physical middleboxes, which are intermediary devices (such as applications, functions, machines, etc.) used to treat data packages dierently from standard routers. In the following, we will give denitions of some network functions, published in memo [2] by Brian Edward Carpenter.

Preliminaries and State-of-the-Art

Network Address Translator (NAT) aims to dynamically assign a global and unique address to a host that does not have one.

Application-level gateway (ALG) needs to keep state for the sessions they are handling, and if the state is lost, the session will normally break irrevocably.

IP Firewall the simplest form of rewall is a router that screens and rejects packets based purely on elds in the IP and transport headers.

Proxies Fielding et al. [START_REF] Fielding | Hypertext transfer protocolhttp/1.1[END_REF], dene a Web proxy as an intermediary program which acts as both a server and a client for the purpose of making requests on behalf of other clients.

Routing schemes

We can distinguish several routing schemes in order to route ow packets through a network. They change in how they transfer messages:

Unicast transmits a message between a sender and a single destination.

Broadcast transfers a message from one sender to all nodes in the network.

Multicast transmits a message from one sender to a group of nodes that have expressed interest in receiving the message.

Anycast transfers a message from one sender to anyone out of a group of nodes.

Geocast delivers a message from a source to a group of nodes in a network based on their geographic location. The hardware middleboxes imply large capital expenditure (CAPEX), as well as an operational expenditure (OPEX) [START_REF] Li | The virtual network function placement problem[END_REF], caused by the considerable number of these devices in the network which is comparable with the number of switches and routers [START_REF] Justine | A survey of enterprise middlebox deployments[END_REF]. Furthermore, the trac has to be conducted in such a way that trac ows pass through these devices, which causes long installation delays and high maintenance costs. Moreover, they do not allow to add new features, are energy-intensive, and have short life-cycles [START_REF] Juliver | Resource allocation in NFV: A comprehensive survey[END_REF]. In average and for a large network, the middleboxes cost can reach up a million dollars over ve years [START_REF] Sherry | Making middleboxes someone else's problem: network processing as a cloud service[END_REF]. Besides, each middlebox provides Preliminaries and State-of-the-Art a dierent functionality i.e., there exists a heterogeneous set of these devices in the network, which necessitates a large management team with dierent expertises.

From 1990 to the year 2000, the network virtualization captivated the attention of network service providers and network operators. Network virtualization can be seen as a representation of one or more logical network topologies, which denes how the data should be transferred, in the same infrastructure. Using virtualization, we can use, for example, multiple logical routers on a single platform, and it also permits resource isolation in terms of CPU and memory.

In addition, network virtualization oers the possibility to transform the way operators plan, use, and maintain their network infrastructures. Its purpose is to deal with a signicant number of specic hardware devices deployed in the operator's networks and the very high generated costs.

In network function virtualization technology the network functions are executed as software on commercial o-the-shelf equipment, this allows their instantiation on demand without the installation of new equipment. For instance, an open-source software-based rewall can be run on an x86 platform in a virtual machine [START_REF] Han | Network function virtualization: Challenges and opportunities for innovations[END_REF].

Network Function Virtualization (NFV) makes it possible for NSPs to employ various Virtual Network Functions (VNFs) without installing new equipment [START_REF] Xu | Online joint placement and allocation of virtual network functions with heterogeneous servers[END_REF]. In [START_REF] Mosharaf | A survey of network virtualization[END_REF],

authors present the history and the state-of-the-art of network virtualization, they present several layers (application, link, network) and levels (link, node) of virtualization and give some examples of virtualized networks. Li and Chen [START_REF] Li | Software-dened network function virtualization: A survey[END_REF] give an overview of the history of NFV, presenting how middleboxes evolve to virtual network functions.

Chiosi et al. in [START_REF] Chiosi | Network functions virtualisation: An introduction, benets, enablers, challenges and call for action[END_REF] present a technical report in which they introduce the network function virtualization, its benets, enablers and challenges. They explain how software replaced hardware in the network architecture. They also give some use cases and raise the challenges of the virtualization. Some works in NFV resource allocation emerged in very recent years. Exhaustive surveys are already given in [START_REF] Juliver | Resource allocation in NFV: A comprehensive survey[END_REF] and [START_REF] Garay | Service description in the NFV revolution: Trends, challenges and a way forward[END_REF]. 1) Network Function Virtualization Infrastructure (NFVI) is the environment that provides the virtual resources required to support the execution of the Virtual Network Functions [START_REF] Chiosi | Network Functions Virtualisation (NFV) Network Operator Perspectives on Industry Progress[END_REF]. It can be physical (servers and switches) or virtual (virtual machines and virtual switches). This layer is divided into three parts:

Herrera and Botero

Hardware resources it includes computing resources (servers and RAM), storage (disks storage), and networking resources (switches, routers, and rewalls).

Virtualization layer it decouples virtual from physical resources; it is responsible for abstracting physical resources into virtual resources. It permits the software to grow separately from the hardware.

Software resources it includes virtual computing resources, virtual storage, and virtual networking resources.

2) Virtual Network Function (VNF) it represents the virtualized network elements such that vrouter, vbase station, vrewall, vIDS, etc. It can be connected

or combined together to oer service chains.

Preliminaries and State-of-the-Art

3) Management and Orchestration (MANO) it includes the orchestration and life-cycle management of physical and/or software resources that support the infrastructure virtualization, and the lifecycle management of VNFs [START_REF] Chiosi | Network Functions Virtualisation (NFV) Network Operator Perspectives on Industry Progress[END_REF]. It is composed of three sotware parts:

Virtualized Infrastructure Manager (VIM) it is responsible for controlling and managing the NFVI computing, network, and storage resources.

VNF Manager: it a software that aims to manage the life-cycle of VNF instances, is responsible for initializing, updating, querying, scaling, and terminating one or more VNF instances.

NFV Orchestrator: it manages the life-cycle of network services, including instantiation, policy management, performance management, and KPI monitoring.

4) Operation Support System/Business Support System (OSS/BSS): OSS handles network, fault, conguration, service, and element management. BSS deals with client, operations, order, billing, and revenue management. In the NFV architecture, the current BSS/OSS of an operator may be combined with the NFV Management and Orchestration using standard interfaces.

Software Dened Networking (SDN)

In traditional telecommunication networks, each switch (device) has its own data plane, to forward the trac, and its own control plane, to decide where to send the data. Also, it has information only about the state of its neighbors to interact with them within the network. No single device has visibility of the whole network. Each switch has to work separately in order to decide where and how to send data, and then all devices synchronize together at a given time slot to route packets through the network.

Supposing that the device uses a command-line interface, this necessitates a manual conguration, which is very costly and time-consuming. Furthermore, having a considerable number of devices implies managing a considerable number of control and data planes. Moreover, if each one has a dierent operating system and interface, it will be challenging to create one application that can be installed on each of them, allowing us to change its control and data plane. Therefore, network management is very challenging [START_REF] Astuto | A survey of software-dened networking: Past, present, and future of programmable networks[END_REF]. This leads to increasing costs in terms of CAPEX and OPEX, in order to minimize the QoS (e.g., end-to-end latency), and makes packet forwarding very complicated within a network [START_REF] Feamster | The road to SDN: an intellectual history of programmable networks[END_REF][START_REF] Pfa | Extending networking into the virtualization layer[END_REF]. Software-dened networking is a new paradigm that simplies network management, and makes the deployment of new services easier by separating the data plane from the control plane. This is done by introducing one controller having a global view of the network state. Several papers were proposed in order to dene, explain, and give examples of Software-Dened Networking (SDN) architecture. A software-based controller is responsible for managing the forwarding information of one or more switches [START_REF] Lara | Network innovation using openow: A survey[END_REF].

As mentioned by Rao in [START_REF] Sridhar | SDN and its use-cases-NV and NFV[END_REF], SDN is a technology that gives the network designer the freedom to refactor the control plane. The author gives the relationship between Network Virtualization, Network Function Virtualization, and Software-Dened Networking, showing how SND plays a signicant role in NV and NFV. A list of companies and their SDN products/solutions is presented in [START_REF] Sridhar | SDN and its use-cases-NV and NFV[END_REF]. Nunes et al. in [START_REF] Astuto | A survey of software-dened networking: Past, present, and future of programmable networks[END_REF], explain that in SDN architecture, the network intelligence is logically centralized in software-based controllers (the control plane), and network devices become simple packet forwarding devices (the data plane) that can be programmed via an open interface.

According to authors in [START_REF]Software-dened networking: The new norm for networks[END_REF], Software-Dened Networking provides a new dynamic Preliminaries and State-of-the-Art network architecture that transforms the traditional networks into service-delivery platforms. They also show the benets of OpenFlow-Based Software-Dened Networks.

OpenFlow is the rst standard protocol interface designed explicitly for SDN, providing high-performance and granular trac control. OpenFlow provides an open protocol to program the ow-table in dierent switches and routers [START_REF] Mckeown | Openow: enabling innovation in campus networks[END_REF], and allow direct access and manipulation of the forwarding plane of network devices. In [START_REF] Manzalini | Software-dened networks for future networks and services[END_REF] authors give a list of the most important technical challenges for the development and deployment of SDN. Guerzoni et al. in [START_REF] Guerzoni | Network functions virtualisation: an introduction, benets, enablers, challenges and call for action, introductory white paper[END_REF], present the relationship between NFV and SDN. Li

and Chen [START_REF] Li | Software-dened network function virtualization: A survey[END_REF] propose a survey that investigates the development of NFV under the software-dened NFV architecture.

Within an SDN architecture, the trac engineering mechanisms become much more eciently implemented with respect to legacy network approaches (e.g., IP, ATM, and MPLS) [START_REF] Tomassilli | Towards Next Generation Networks with SDN and NFV[END_REF]. Furthermore, all network information are easily retrieved; network elements are dynamically and proactively programmed without having to handle them individually [START_REF] Ian F Akyildiz | A roadmap for trac engineering in sdn-openow networks[END_REF]. Moreover, SDN provides various benets, namely, fast detection time [START_REF] Sharma | Openow: Meeting carrier-grade recovery requirements[END_REF], re-routing [START_REF] Berde | ONOS: towards an open, distributed SDN OS[END_REF] and shows its capacity to identify and grow from failure under the requirement below 50 ms [START_REF] Brungard | Requirements of an mpls transport prole[END_REF]. With SDN, the network becomes programmable, dynamic, and exible according to frequent network changes.

Examples of centralized controllers are: Beacon [START_REF] Erickson | The beacon openow controller[END_REF], Ryu [128], OpenDayLight [START_REF] Medved | Opendaylight: Towards a model-driven SDN controller architecture[END_REF],

and Maestro [START_REF] Cai | Maestro: A system for scalable openow control[END_REF]. Nevertheless, having a single controller in the network signies having a unique point of failure; if the node (controller) fails, there is no other point that makes decisions and manages the data plane. Accordingly, the notion of a distributed controller is proposed. It represents a set of nodes that can be either centralized or physically distributed. Some distributed controllers are Onix [START_REF] Koponen | Onix: A distributed control platform for large-scale production networks[END_REF], ONOS [START_REF] Berde | ONOS: towards an open, distributed SDN OS[END_REF], and DISCO [START_REF] Phemius | Disco: Distributed multidomain SDN controllers[END_REF]. Application plane allows network managers to quickly congure, manage, secure and optimize network resources via dynamic automated SDN programs.

Control plane exercises direct control over the data plane using an Application

Programming Interface (API), which denes the information exchange between the two planes.

Data plane is a set of network elements such as routers, switches, and middleboxes that oer ecient and programmable packet forwarding devices without any software to make autonomous decisions.

OpenFlow is a protocol represented by an open interface that manages to update the ow table of the device.

Service Functions Chaining (SFC)

Service Function Chain (SFC) is dened as a sequence of network functions that should be traversed by a given data ow in a predened order [START_REF] Quinn | Problem statement for service function chaining[END_REF]. For example, the network administrator may specify a policy that all HTTP trac should follow the policy chain:

rewall → IDS → proxy [START_REF] Li | The virtual network function placement problem[END_REF]; Also for instance, the Intrusion Detection System must inspect packets before compressing or encrypting them [START_REF] Savi | Impact of processing costs on service chain placement in network functions virtualization[END_REF].

For each service, a specic set of required functions needed to handle its packets, and their order is dened by the service function chain; this allows exible management and classication of services and permits to dene the requirements associated with each one [START_REF] Ahmed M Medhat | Service function chaining in next generation networks: State of the art and research challenges[END_REF].

Traditionally, the service function chain represents an ordered sequence of hardware devices that must be traversed by the data ows to support some service. Nevertheless, new hardware devices should be used, placed, and connected to the network elements if a new service is required. This generates additional costs and is time-consuming [START_REF] Li | Software-dened network function virtualization: A survey[END_REF].

Thanks to NFV and SDN, the VNFs could be installed on-demand on network nodes, and the routing paths associated with the Service Function Chains could be computed In the Virtual Network Embedding problem, the service chains can be seen as virtual graphs (paths) to embed into the original network. Each node in the virtual graph represents a VNF, and the links between nodes represent the order. Bhamare et al. in [START_REF] Bhamare | A survey on service function chaining[END_REF] dene Service Function Chaining as a mechanism that allows various service functions to be connected to form a service allowing operators to benet from the software-dened virtualized infrastructure. Medhat et al. in [START_REF] Ahmed M Medhat | Service function chaining in next generation networks: State of the art and research challenges[END_REF], explain that the optimization models are needed for SFC conguration and allocation to achieve optimal network performance, satisfy user demands, accommodate SFCs dynamically, and minimize network costs. Zhang et al. [START_REF] Zhang | Enabling ecient service function chaining by integrating NFV and SDN: architecture, challenges and opportunities[END_REF] proposes a typical SFC functional framework integrating SDN and NFV including service modeling and resource allocation.

Literature review 35

Many works were also proposed for the Service Functions Chaining; for example, a near-optimal solution is provided by a heuristic in [START_REF] Li | Network functions virtualization with soft real-time guarantees[END_REF] for the real-time NFV, aiming to maximize the total number of requests assigned to the cloud for each SFC, while the deadline constraints are satised. Anity-based approach to overcome with limited applicability of the model.

Hantouti et al. [START_REF] Hantouti | Trac steering for service function chaining[END_REF] present a comprehensive analysis and classication in three categories of the current SDN-based SFC approaches using eciency criteria such as deployment cost, scalability, and exibility.

Some studies have been proposed in the literature in order to dene how to place an ordered chain of functions in a network. Dierent objective functions were considered such as: minimizing the number of activated nodes, the end-to-end latency [START_REF] Mehraghdam | Specifying and placing chains of virtual network functions[END_REF], the energy consumption [START_REF] Huin | Energyecient service function chain provisioning[END_REF], or the bandwidth [START_REF] Huin | Optimal network service chain provisioning[END_REF].

Problem variants

Dierent objective functions variants In the context of NFV and SDN, the underlying problem consists of nding an optimal installation of VNFs, so that the given trac requests can be routed within the given network infrastructure while respecting the SFC constraints [5,[START_REF] Gupta | Joint virtual network function placement and routing of trac in operator networks[END_REF][START_REF] Mehraghdam | Specifying and placing chains of virtual network functions[END_REF][START_REF] Moens | VNF-P: A model for ecient placement of virtualized network functions[END_REF]. Many dierent variants of the problem have been studied in the recent literature. When it comes to the denition of the objective function, some authors consider the minimization of the number of activated nodes, the VNFs installation costs [START_REF] Allybokus | Virtual function placement for service chaining with partial orders and anti-anity rules[END_REF], the end-to-end latency [START_REF] Mehraghdam | Specifying and placing chains of virtual network functions[END_REF], the energy consumption [START_REF] Huin | Energyecient service function chain provisioning[END_REF],

or the bandwidth [START_REF] Huin | Optimal network service chain provisioning[END_REF].

Simplest variants Several straightforward variants of the problem were studied, for example, in [7], authors consider the VNF chaining problem with a single type of VNFs for all demands. Furthermore, they suppose that all trac requests represent exactly the same service. They propose an extension of their work that allows using multiple VNFs and dierent services. Also, Casado et al. [START_REF] Casado | Virtualizing the network forwarding plane[END_REF] consider the problem with a single type of VNF and present a heuristic algorithm towards solving the placement problem.

Basta et al. [START_REF] Basta | Applying NFV and SDN to LTE mobile core gateways, the functions placement problem[END_REF] [START_REF] Cohen | Near optimal placement of virtual network functions[END_REF], provide a near-optimal approximation algorithm to address the problem of placing VNFs on the physical network without order constraints. Their objective is to minimize installation and path costs. Moreover, in order to minimize the number of used network functions, Sang et al. [START_REF] Sang | Provably ecient algorithms for joint placement and allocation of virtual network functions[END_REF] x the routing paths and focus their attention on the placement of VNFs without considering any chaining constraints.

Closest variants Allybokus et al. [START_REF] Allybokus | Virtual function placement for service chaining with partial orders and anti-anity rules[END_REF] consider a generic version of the VNF placement and routing problem in which many features (partial order on the functions, bound on end-to-end latency, conicts between network functions, resource limitations and processing capacities) are already taken into account. They proposed several MILP formulations, either to minimize the total deployment cost or to minimize the number of rejected demands. A heuristic algorithm based on a continuous relaxation of one formulation is also proposed and appears to be very ecient on instances derived from the Geant topology. However, the MILP models do not seem to be strong enough for direct resolution with the Cplex solver. Similarly, Addis et al. [5] formulate the VNFPRP as an MILP, taking into account ow constraints, latency constraints, node and function capacity constraints, while minimizing the number of used links to route trac and the number of used nodes to install functions. The authors also propose a math-heuristic approach for the problem. Contrary to our work, they allow for compression/decompression at the VNF nodes.

Preliminaries and State-of-the-Art Furthermore, Mehraghdam et al. [START_REF] Mehraghdam | Specifying and placing chains of virtual network functions[END_REF] showed that the problem could be considered as the Location-Routing Problem [START_REF] Nagy | Location-routing: Issues, models and methods[END_REF], that aims to create several paths between dierent VNFs, and then connect them to get source-destination paths satisfying precedence constraints. We test the proposed model on a set of realistic instances derived from SNDlib.

Motivations

The Also, Huawei, in the last years, deployed 560 SDN/NFV commercial projects around the world [START_REF]Huawei Releases SDN/NFV Commercial and Technological Innovations[END_REF]. Furthermore, in 2013, Google declared the use of SDN to interconnect its data centers all over the planet [START_REF] Jain | B4: Experience with a globally-deployed software dened wan[END_REF]. Using this technology, Google was able to achieve several advantages, including ecient network control, more accessible and faster innovation cycles of networks and services, better network exploitation, and a decrease of both OPEX (Operating Expenditure) and CAPEX (Capital Expenditure).

Both technologies create new challenging problems and on eof them is the Virtual

Network Functions Placement and Routing problem. This problem can be decomposed into two NP-hard problems, namely, the multi-Commodity Flow (MCF) problem with unsplittable commodity paths and the Facility Location (FL) problem (as shown in Figure 2.1). Further constraints can be added, such as, the service chain and conict constraints. All these technical constraints make the problem very challenging. Specifically, the challenges are: (i) where to install the required VNFs on network nodes at minimum cost (ii) dening, for each demand, a routing path that traverses suitable nodes in the correct order to satisfy the requirements of a given service chain, and

(iii) taking into account network load and other dynamic characteristics when routing through existing VNFs [START_REF] Dwaraki | Adaptive service-chain routing for virtual network functions in software-dened networks[END_REF].

The methods proposed in the literature to deal with the problem mainly consist in heuristics, or two-stage methods, that solve the routing part and the installation part of the problem separately. Alternatively compact MILP formulations have been proposed that cannot solve large scale instances (see, [5,[START_REF] Allybokus | Virtual function placement for service chaining with partial orders and anti-anity rules[END_REF][START_REF] Gouareb | Delay sensitive virtual network function placement and routing[END_REF][START_REF] Hmaity | Latency-and capacity-aware placement of chained virtual network functions in fmc metro networks[END_REF][START_REF] Tashtarian | Distributed VNF scaling in large-scale datacenters: An ADMM-based approach[END_REF]). Generally, to overcome the fairly limited computational performance of compact formulations, and to improve the capabilities of nding exact solutions for larger instances, reformulation and decomposition methods are recommended in the literature (see [START_REF] Bilde | Sharp lower bounds and ecient algorithms for the simple plant location problem[END_REF], [START_REF] Chopra | An extended formulation of the convex recoloring problem on a tree[END_REF], [START_REF] Sadykov | Column generation for extended formulations[END_REF], [START_REF] Vanderbeck | Decomposition and column generation for integer programs[END_REF], [START_REF] Vanderbeck | Reformulation and decomposition of integer programs[END_REF]).

The Virtual Network Functions Placement and Routing Problem

Motivated by all of the above, the main objective of this thesis is to model and solve the VNFPRP using exact and heuristic methods. To this end, we use some decomposition and reformulation approaches from the literature in order to solve large scale instances. Outline of the chapter The chapter is organized as follows. In Section 2.2, we dene the problem and introduce the main notations. In Section 2.3, we present the problem's properties. An illustrative example is shown in Section 2.4. Section 2.5 is devoted to the problem's complexity analysis. The compact MILP formulation is introduced in Section 2.6. In Section 2.7, we present some computational results and we end the chapter by some conclusions and remarks in Section 2.8.

Problem denition 2.2.1 Notation

The telecommunication network is modeled as a bi-directed graph G = (N, A). The set of all physical locations equipped with hardware devices allowing VNFs installation is denoted by N and called the set of nodes. The set representing all connections between nodes is called the set of arcs and is denoted by A. At each node u ∈ N at most c u (c u ∈ N) VNFs can be installed, and an activation cost ψ u > 0 has to be paid. The arc latency l uv ≥ 0 is dened for each arc (u, v) ∈ A.

Let F denote the set of all virtual network functions. Each VNF f ∈ F has a limited (bandwidth) capacity m f ∈ N to manage the trac demand. For installing each single copy of f at a node u ∈ N , an installation cost ψ f u > 0 has to be paid. For the set of all trac requests/commodities, denoted by C, each commodity k ∈ C is characterized by: a source node s k ; a destination node d k (s k = d k ); a bandwidth b k > 0; a maximum latency value l k > 0; a subset of VNFs, F k ⊆ F , and a set of incompatibility constraints A k between VNFs. Finally, recall that each commodity requires a specic Service Function Chain to handle its data packages in a specic order. The fact that the VNF f has to be executed before VNF g is expressed by f ≺ k g in the path associated to the commodity k. End-to-end latency constraints the sum of arc latencies belonging to the routing path of each commodity k ∈ C should not exceed the given upper-bound

Parameters

m f : Capacity of VNF f , f ∈ F . c u : Capacity of node u, u ∈ N . l uv : Latency of links (u, v), (u, v) ∈ A. ψ f u : Installation cost of VNF f at node u, f ∈ F , u ∈ N . ψ u : Activation cost of node u, u ∈ N . s k : Source node for commodity k, k ∈ C. d k : Destination node for commodity k, k ∈ C. l k : Maximum latency delay of commodity k, k ∈ C. b k : Bandwidth of commodity k, k ∈ C.
l k ∈ R + . Installation constraints each VNF f ∈ F k required for commodity k ∈ C
should be installed at one of the nodes of the routing path, but not at the source node.

Precedence constraints for each commodity k ∈ C, the VNFs composing its Service Function Chain should be traversed in the right order by the routing path.

Conict constraints (also called anti-anity/incompatibility constraints): For dierent reasons (e.g. resiliency, privacy) [START_REF] Tomassilli | Energyecient service chains with network function virtualization[END_REF] VNFs pairs that are in conict should not be installed at the same node.

The VNFPR problem aims to minimize the sum of VNF installation costs ψ f u plus the sum of node activation costs ψ u . A solution satisfying the constraints described above is said to be feasible for the problem. Furthermore, if the associated cost is minimized, then, it is called optimal.

Proprities of the VNFPRP

In this thesis, we consider the following assumptions:

In NFV environment, several, heterogeneous and similar virtual network functions can be installed at the same node in the graph.

Proprities of the VNFPRP 47

Link capacities are omitted, due to the fact that we are not dealing with the strategic network design phase (which has to be done much earlier in the planning process), but with a tactical planning. The network is designed (by making decisions regarding which links need to be installed, and with which capacity) so that it can accommodate a large amount of trac. This is a long term process involving expensive investments and possible expansion of capacities (if it turns out that some links are heavily loaded). On a shorter time range, various services are managed to make the best possible use of existing resources, but the transmission/transport network is usually not the bottleneck (as the service demands are not of the same order of magnitude as the link capacity). Typically, network capacities are huge and shared by many dierent services. Thus, in this thesis we safely assume that there are enough link capacities in the existing network to handle all demands.

We consider the unicast routing scheme.

All routing paths should be elementary (without circuit). The use of elementary paths is linked to the protocols used, and to a choice of operators to simplify the supervision of ows.

Each commodity requires at least one VNF, otherwise, these commodities can be eliminated (the routing problem can be seen as the shortest path problem).

The Virtual Network Functions Placement and Routing Problem Pairs of nodes (1, 20), [START_REF] Allybokus | Virtual function placement for service chaining with partial orders and anti-anity rules[END_REF][START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF], and (10,[START_REF] Basta | Applying NFV and SDN to LTE mobile core gateways, the functions placement problem[END_REF] represent (source, destination) nodes associated with RA, WV and CM, respectively. From Figure 2.2, the red routing path associated with the service RA, represents an elementary path passing rst through the node 12 on which the VNF Firewall is installed rst, and then through node 16 at which the VNF Parental control is installed. In the same way, the green (resp. yellow) path goes from the source node 11 (resp. 10) to the destination node 19 (resp. 18) and passes through the nodes on which the VNFs are installed in the order satisfying the SF C 2 (resp. SF C 3 ) constraints. Moreover, we can see that services RA and CM share the same VNF Firewall. Also, services RA and WV share the VNF Parental control. Furthermore, from the example, we observe that the VNFs associated to each service cannot be installed on the source node, while these VNFs can be installed on the destination node.

Illustrative example

Complexity analysis

Problem complexity

In this section we show that the problem under consideration is strongly NP-hard, even for the variant in which no conicts between functions and no node and VNFs capacities are imposed.

Theorem 2.1 The Uncapacitated Virtual Network Functions Placement and Routing

Problem is strongly NP-hard, even for a single commodity and without latency, conict, and precedence constraints.

Proof. The proof is done by a polynomial reduction from the Uncapacitated Facility Location Problem (UFLP), which is an NP-hard problem in a strong sense (see, Theorem 3.1 in [START_REF] Cornuéjols | The uncapicitated facility location problem[END_REF]).

Given a set of n facilities (sites) and m customers, let ψ i be the cost of opening the facility i and ψ j i the cost of assigning customer j to facility i. We suppose that all facilities have unlimited capacity. The UFLP consists in nding which of the n facilities to open and how to assign the customers to open facilities so that the facility opening cost plus the assignment cost are minimized.

To reduce the UFLP to the uncapacitated VNFPR problem, we consider an arbitrary connected graph G = (N, A) (G can be a complete graph) with |N | = n nodes (with a

The Virtual Network Functions Placement and Routing Problem one-to-one correspondence between facilities from I and nodes in N ) and a set F of m virtual network functions. We use one commodity |C| = 1, and let s 1 = 1 and d 1 = n be two distinct nodes from N . We set the latency l 1 = +∞ and assume that there are no precedence constraints (i.e., any ordering of functions from F is feasible). The transformation details are given in Table 2 This transformation is polynomial in the number of facilities and customers. An instance of the uncapacitated VNFPR dened in this way consists of installing all virtual network functions from F on a subset of nodes from N , while minimizing costs of installation of VNF and node activation costs. Thus, there is a one-to-one correspondence between an optimal solution of the uncapacitated VNFPR and the optimal solution of the uncapacitated facility location problem with exactly the same solution value. This proves that the Uncapacitated Virtual Network Functions Placement and Routing problem is NP-hard.

.2. UFLP 1-to-1 corresp. Uncapacitated VNFPR - ⇐⇒ |C| = 1 I ⇐⇒ N J ⇐⇒ F 1 Facility i ⇐⇒ Node i Customer j ⇐⇒ Virtual network function j Open facility i ⇐⇒ Activate node i Customer j is supplied by facility i ⇐⇒ Function j is installed at node i ψ i cost of opening the facility i ⇐⇒ ψ i activation cost at node i ψ j i cost of assigning cust. j to facility i ⇐⇒ ψ j i installation cost of function j at i - ⇐⇒ G=(N, A), connected graph - ⇐⇒ s 1 = 1, d 1 = n, l 1 = +∞

Compact MILP formulation

In this section, we propose a compact (i.e., polynomial in size) MILP formulation to model the VNFPRP. We st describe the variables necessary to model it and then the set of constraints.

2.6 Compact MILP formulation 51

Decision variables

The set of variables required in our MILP formulation is described in Table 2.3:

Variables Type

x f k u 1, if

Mathematical model

The VNFPRP can then be modeled as follows:

(P ) : min

u∈N f ∈F ψ f u z f u + u∈N ψ u w u (2.1) (u,v)∈A t k uv - (v,u)∈A t k vu =        -1 if u = d k , 1 if u = s k , 0 otherwise. k ∈ C, u ∈ N (2.2) (u,v)∈A t k uv l uv ≤ l k , k ∈ C (2.3) f ∈F z f u ≤ c u w u , u ∈ N (2.4) k∈C y f k u b k ≤ m f z f u , f ∈ F, u ∈ N (2.5) y f k u + y gk u ≤ 1, k ∈ C, (f, g) ∈ A k , u ∈ N (2.6) 52
The Virtual Network Functions Placement and Routing Problem

(t k uv -1) + (x f k v -x f k u ) ≤ y f k v , k ∈ C, f ∈ F k , (u, v) ∈ A (2.7) x gk u ≤ x f k u , k ∈ C, f, g ∈ F k : f ≺ k g, u ∈ N (2.8) u∈N y f k u = 1, k ∈ C, f ∈ F k (2.9) y f k u ≤ (v,u)∈A t k vu , k ∈ C, f ∈ F k , u ∈ N (2.10) y f k u ≤ x f k u , k ∈ C, f ∈ F k , u ∈ N (2.11) x f k u = 1, if u = d k 0, if u = s k k ∈ C, f ∈ F k (2.12) y f k u , x f k u ∈ {0, 1} u ∈ N, f ∈ F, k ∈ C (2.13) w u ∈ {0, 1} u ∈ N (2.14) t k uv ∈ {0, 1} (u, v) ∈ A, k ∈ C (2.15) z f u ∈ N u ∈ N, f ∈ F (2.16)
The objective function ( the conict constraints and they guarantee that two VNFs in conict are not installed at the same node u ∈ N . Constraints (2.7) are needed to link, the node-installation variables (y), the precedence variables (x) and the arc variables (t): the left-hand-side is forced to 1 (implying that the function f is installed at the node v) if and only if (i) the arc (u, v) is taken in the path associated with the considered commodity k and (ii) the function f is installed at or before the node v and it is not installed at or before the node u. Constraints (2.8) impose the VNFs order for each commodity. Constraints (2.9) guarantee that all required functions for commodity k ∈ C are installed at the graph nodes. Constraints (2.10) ensure that if a VNF f ∈ F k is installed at a node u for a given commodity k, then the associated routing path must enter that node.

Inequalities (2.11) link the precedence and the installation variables, x and y, and express the fact that if VNF f is installed at node u for the commodity k, then f is installed at or before the node u. Finally, constraints (2.12) guarantee that, for each Theorem 2.2 The model (P) is valid. In particular, there always exists an optimal solution of (P) such that for each k ∈ C, the ow variables t k (representing the routing of commodity k) correspond to an elementary path in G.

Proof. Let k ∈ C be a given commodity. Flow constraints (2.2) guarantee that s k and d k are connected. Furthermore, latency limit l k cannot be violated due to the constraint (2.3). It only remains to prove that there exists an optimal solution of the compact model whose routing path for k is elementary, which will imply that precedence constraints (2.8) are satised (due to constraints (2.7)).

Let (x, ŷ, ŵ, t, ẑ) be an optimal solution of the compact model, and let G k be the subgraph of G induced by the arcs (u, v) ∈ A such that tk uv = 1. Assume that G k is not an elementary path connecting s k to d k . We distinguish the following two cases: Case 1: subgraph G k is connected, but it contains cycles (see Figure 2.3). Assume that the routing path is composed by an elementary path and (without loss of generality) a single cycle. The nodes in the elementary path are denoted by u i , ∀i = 1, . . . , q + 1 (where u q+1 = d k ) and all the other nodes are denoted by v i , ∀i = 1, . . . , p. We will show that all VNFs f ∈ F k must be installed along the elementary path, so that the nodes v i , i = 1, . . . , p can be removed from the routing path, without violating feasibility and optimality of the solution. Indeed, in that case, removing the cycle from G k only reduces the latency of the routing path, the ow balance constraints remain valid and the VNFs are properly placed at the remaining nodes.

Let f be an arbitrary VNF from F k and assume the contrary, i.e., that f is installed at a node v and it is not installed at any node u. Therefore:

ŷfk u i = 0, ∀i = 1, . . . , q + 1.
Using the precedence constraints (2.7) tuv -

1 + xfk v -xfk u ≤ ŷfk v we can now calculate the values of xfk u i
, for all i = 1, . . . , q + 1. Let us start with the arc (s k , u 1 ), the constraint (2.7) reads as follows:

tk s k u 1 -1 + xfk u 1 -xfk s k ≤ ŷfk u 1
Recall that VNFs cannot be installed at the source node, i.e., xfk s k = 0 (cf. constraint (2.12)), and that the arc (s k , u 1 ) is taken in the path so tk s k u 1 = 1. After replacing the variables by the values we obtain:

1 -1 + xfk u 1 -0 ≤ 0 ⇒ xfk u 1 ≤ 0.
So, it follows that xfk u 1 = 0. By repeating these steps for all subsequent arcs (u i , u i+1 ), i = 1, . . . , q, we obtain:

xfk u i = 0, ∀i = 2, . . . , q + 1.
Case 2: Subgraph G k is disconnected (see Figure 2.4). In that case, following similar arguments as for the Case 1, we conclude that VNFs from F k can be installed only along the elementary path in G k connecting s k to d k , so that without loss of feasibility and optimality, the remaining nodes from G k can be removed (and the values of t can be appropriately redened).

According to this proof, we conclude that we can always nd an optimal solution of the compact model which is composed by elementary paths only.

Model analysis

In this subsection we analyze the proposed compact formulation in terms of the number of variables and constraints generated by the MILP model.

Number of variables

The total number of variables is:

|C| × |A| + 2 × |N | × |C| × |F k | + |N | + |N | × |F | (2.17) ⇐⇒ O(|C| × |A|) + O(|N | × |C| × |F k |) + O(|N | × |F |) which is equivalent to O(|C| × (|A| + |N | × |F |)

Number of constraints

The total number of constraints is:

|C| × |N | + |C| + |N | + |N | × |F | +|N | × |C| × |S k | + 3 × |C| × |F k | +|C| × |A| × |F k | + |N | × |C| × |F k | 2 +3 × |N | × |C| × |F k | (2.18) which is equivalent to O(|C| × |F | × (|A| + |N | × |F |)) 56
The Virtual Network Functions Placement and Routing Problem

Computational results

In this section, we test the capability and the performance of the compact MILP formulation on a set of realistic instances derived from the SND library in terms of CPU time, and nal gaps. The instances creation is detailed in the next chapter. All the experiments described in this section were made using a computer with Intel(R)Xeon(R) CPU E5-2650 v2 processor clocked at 2.60GHz, 32 cores, two threads per core and 252GB RAM, under Linux operating system. The MILP compact formulation is implemented using the Python API for CPLEX, which is run in single-thread mode and a default memory limited to 20GB. All CPLEX parameters were set to their default values. A default time limit of one hour is set for each tested instance.

In the following, we denote by C, The compact MILP formulation. Outline of the chapter. Section 3.1 is devoted to the path formulation. In Section 3.2, we present computational results and a sensitivity analysis in which we vary the input parameters and discuss the results. In Section 3.3 we provide some concluding remarks.

Path-based MILP Formulation

In this section, we present the path-based MILP formulation to model the VNFPRP.

We rst describe the set of decision variables, and then the constraints dening the model.

Decision variables

The path formulation is characterized by ve families of variables described in Table 3 Routing constraints

p∈P k λ k p = 1, k ∈ C (3.1)
For each commodity k, exactly one s k -d k elementary routing path p ∈ P k satisfying the latency constraints is chosen.

Path-based MILP Formulation

Precedence constraints

( p∈P k (u,v)∈p t pk uv λ k p -1) + (x f k v -x f k u ) ≤ y f k v , k ∈ C, f ∈ F k , (u, v) ∈ A (3.2)
Inequalities (3.2) represent the linking constraints between routing variables (λ), installation variables (y) and precedence variables (x), ensuring that for each commodity k: (i) if the routing path p passes through the arc (u, v), and (ii) the VNF f is installed at or before the node v and (iii) the VNF f is not installed at or before the node u, then the left-hand-side is forced to be equal to 1 (imposing the installation of the VNF f at the node v).

Linking constraints

y f k u ≤ (v,u)∈A p∈P k t pk vu λ k p , k ∈ C, f ∈ F k , u ∈ N (3.3) Constraints (3.
3) impose that the chosen path for each commodity k passes through the nodes where the required functions are installed.

Mathematical model

The Path-based formulation (PF) then reads as follows:

(PF): min

u∈N f ∈F ψ f u z f u + u∈N ψ u w u (λ, x,
λ k p ∈ {0, 1} k ∈ C, p ∈ P k x f k u , y f k u ∈ {0, 1} k ∈ C, u ∈ N, f ∈ F w u ∈ {0, 1} u ∈ N z f u ∈ N u ∈ N, f ∈ F
This model contains a possibly exponential number of variables. We therefore consider a heuristic approach in which only a subset of the most promising feasible paths is selected per each commodity. otherwise it provides a heuristic solution for VNFPRP.

Getting the routing paths

In order to generate routing paths associated with each commodity k, we use Yen's algorithm [START_REF]Finding the k shortest loopless paths in a network[END_REF]. This algorithm aims to nd a limited number of loopless shortest paths between a pair of nodes. Yen's algorithm is based on Dijkstra's algorithm [START_REF] Edsger | A note on two problems in connexion with graphs[END_REF]. For reasonably small instances, instead of recovering a xed number of paths per commodity, we run Yen's algorithm to get all s k -d k paths.

In order to consider the latency constraints 2.3 in the path formulation model using Yen's algorithm, we keep only paths whose length is less or equal to l k , k ∈ C.

Remark: We used NetworkX graph library of Python to get all the paths associated with each commodity. The used Python function is called Shortest_simple_paths, which returns a list of all simple paths between a source node and a destination node, ordered from the shortest to the longest one.

To x the number of generated paths to κ we use the following Python command:

list(islice(nx.shortest_simple_paths(G, source, destination, weight = weight), κ))

We adapted the Python function Shortest_simple_paths to get only the latency constrained paths. In what follows in this manuscript, when we say: we get all the paths using Yen algorithm and we keep only the latency constrained one, we mean that we use this adapted Python function Shortest_simple_paths.

By denition of the sets P k , only elementary paths are considered in this model, and hence, the formulation is correct. Unfortunately, the number of elementary paths can be exponential, which makes this MILP formulation intractable, unless it is embedded within a Branch-and-Price procedure. Nevertheless, this formulation can be eciently used as an MILP heuristic, considering only a subset (polynomial in size) of elementary paths in the model. Proof. Let X = ( λ, x, ŷ, ẑ, ŵ) denotes the optimal solution of (P F ) with λ fractional, i.e., there exists at least one commodity k ∈ C for which the path variable λk p is fractional. The proof aims to show that from the optimal fractional solution, X a totally integer solution can be constructed. Let X = ( λ, x, ŷ, ẑ, ŵ) be the integer solution which is also feasible for (P F ).

Variables ŷ and ŵ appear in the objective function with strictly positive cost and are the same for both solutions X and X; thus, the two solutions will have the same objective value.

Let k ∈ C be the fractional commodity for which the associated vector λ contains fractional values in X. Let P k = {p 1 , p 2 , . . . , p Q } be the set of fractional paths associated with commodity k, satisfying the path constraints (3.1). In that case, by setting an arc upper bounds corresponding to constraints (3.2) and node capacity lower bounds corresponding to constraints (3.3) as follows:

ĉ ap uv := min f ∈F k {ŷ f k v + 1 -xfk v + xfk u } (u, v) ∈ A (3.4) ĉ ap u := max f ∈F k {ŷ f k u } u ∈ N (3.5)
we obtain the following linear inequality system:

0 ≤ p∈P k λ k p t pk vu ≤ ĉ ap uv (u, v) ∈ A (3.6) (v,u)∈A p∈P k λ k p t pk vu ≥ ĉ ap u u ∈ N (3.7)
For each node u ∈ N such that ĉ ap u = 1, by the feasibility of λ, all paths p ∈ P must pass through the node u. Hence, to satisfy constraints (3.6) and (3.7), one could take any of the paths p ∈ P k as a feasible binary solution.

From the above, amid the fractional paths constituting λ, there exists at least one (integer) path which can be used to replace λ. The same procedure can be reiterated for each fractional commodity that admits a fractional λk p in the solution X, keeping the same structure of the binary component (x, ŷ, ẑ, ŵ), due to the fact that the routing paths are chosen separately for each commodity.

Computational results

The purpose of this section is to test the eciency and the sensitivity of the proposed path formulation. We present some computational results to compare the path formulation with the compact MILP formulation proposed in Chapter 2. To compare both models, we focus on the CPU time, the quality of the obtained solutions, and the nal gaps between the global lower and the best known upper bound. We vary the problem input parameters and conduct a sensitivity analysis to determine the most relevant parameters that aect the model's empirical performance.

All the experiments described in this section were made using a computer with Intel(R)Xeon(R) CPU E5-2650 v2 processor clocked at 2.60GHz, 32 cores, 2 threads per core and 252GB RAM, under Linux operating system. The path formulation is implemented using the Python API for CPLEX, which is run in single-thread mode and a default memory limited to 20GB. All CPLEX parameters were set to their default values. A default time limit of one hour is set for each tested instance.

Benchmark instances

Our benchmark instances are generated using the SND library [125] (Survivable Network Design), which is a repository of realistic telecommunication network design instances. The number of nodes, arcs, and demands varies for each instance. SNDlib provides the graph topology, along with the node coordinates and a set of demands with the associated source node, destination node, and a bandwidth. The remaining parameters required for our setting are generated as follows:

To calculate the distance between two nodes u and v in km having only their longitude ϕ u , ϕ v and latitude ς u , ς v , we use the Spherical Law of Cosine [START_REF] Monawar | Anti-theft vehicle tracking and regaining system with automatic police notifying using haversine formula[END_REF]. First we use the following formula to calculate the angular distance in radians:

S uv = arc cos(cos ς u cos ς v + sin ς u sin ς v cos dϕ) with, dϕ = ϕ v -ϕ u .
The distance in km is then obtained as:

d uv = R × S uv ,
where R is earth's radius (R = 6378137km). The ber propagation delay per km is roughly equal to 10µs/km, see, e.g., [START_REF] Chitimalla | 5g fronthaul-latency and jitter studies of cpri over ethernet[END_REF][START_REF] Line Mp Larsen | A survey of the functional splits proposed for 5G mobile crosshaul networks[END_REF]. To dene the latency l uv of an arc (u, v), we multiply the distance between u and v by the ber propagation delay, and we set l uv = l vu . We consider the following set of six Virtual Network Functions, typically employed in service function chaining [START_REF] Huin | Optimal network service chain provisioning[END_REF][START_REF] Savi | Impact of processing costs on service chain placement in network functions virtualization[END_REF] to construct our VNFs set, namely F = {NAT, FW, TM, WOC, IDPS, VOC }. A detailed description of these VNFs is given in Table 3. We divide the set of demands in ve categories: Online Gaming, Video Streaming, Voice over IP, Web Services, and Other Services. Each category is characterized by a latency value and a set of chained service functions, as depicted in Table 3.3.

To dene the latency value l k and the set of chained VNFs associated with each commodity k, the demands are rst assigned to exactly one category. To do so, we calculate the shortest path SP k between s k and d k , which is equal to the sum of the arc latencies composing it. Based on its length, we randomly assign our demands in one of the ve categories, and so we set the value of l k and the set of chained VNFs. For example, if the length of the shortest path is equal to l(SP k ) = 40ms, then we can assign the demand to any of the ve categories, and we randomly choose one. For each commodity k ∈ C there is at most one anti-anity constraint (AAC) between VNFs. We suppose that we cannot install Firewall and Network Address

Latency value

Translator at the same node. The maximum number of AAC is xed to ve per instance because adding multiple anti-anity constraints renders some instances infeasible. To test the eciency of the proposed path formulation, we consider fteen dierent graphs from SNDlib, each of them corresponding to one instance-type in our benchmark. Table 3. 4 summarizes the details about the graphs used in this study. In order to vary the demand types per instance, we generated ten instances for each instance type. Table 3

Models analysis

In this subsection, we analyze the proposed path-based formulation in terms of the number of variables and constraints. We compare them to the ones of the compact MILP formulation introduced in Chapter 2. In the following, we show only the dierence in terms of number of variables and constraints between both models.

Variables: We observe that the number of variables in the path formulation is, on average, more signicant than the number of variables of the compact MILP formulation. This is due to the number of feasible paths associated with each commodity. We can also observe that the number of constraints in the compact MILP formulation is, on average higher than the one in the path formulation this number depends on the number of nodes and demands.

Obtained results

The following two settings are compared in our computational study:

C : The compact MILP formulation proposed in Chapter 2.

PF : The path-based formulation proposed in Section 3.1.

Table 3.7 shows the results obtained by solving 10 instances for each instance type. The CP U time, GAP , and the costs of each instance type are obtained by calculating the average over all ten instances for which a feasible solution was found, otherwise these instances are not considered. We xed the maximum number of latency-constrained paths generated by Yen's algorithm to 5000 paths per commodity. In order to dene the nature of the path formulation, for each instance type, we pre-calculate the maximum number of latency constrained paths over all commodities. This number is shown in column #paths in Table 3.6. Based on that, the second column denoted by PF(E/H), describes the nature of the path formulation: exact E means that the maximum number of feasible paths was below our threshold of 5000 paths per commodity; or heuristic H method, meaning that not all feasible paths were taken into consideration by the path formulation, which is the case only for two instance types: Dfn-bwin and Dfn-gwin.

In Table 3.6, we also report the number of instances solved to optimality by the two models. The path formulation solves eight instances to optimality within one hour, whereas C manages to nd the optimum for only 4 instances. Moreover, there are ve instances for which the compact MILP formulation cannot provide any feasible solution within the given time limit. On the contrary, feasible solutions are found for all instances, using the path formulation. Table 3.7 summarizes the results obtained by solving 15 SNDlib instance types with dierent graph topologies and dierent number of commodities. Columns one and two in the table represent the average of the CPU time calculated from instances solved to optimality by PF and C, respectively. We observe that for instances solved to optimality, PF is outperforming C in terms of CPU time. However, for most of the instances, none of the two methods manages to provide an optimal solution within one hour. Notion T L in the CPU time column means that the time limit was reached. When the optimal solution is not found after reaching the time limit, we report the average nal For all instances for which PF is an exact method, the nal gap provided (CPLEX exit gap) by PF is consistently better than the nal gap provided by C (see Figure 3.2).

Each coordinate (x, y) in Figure 3.3 indicates that for y instance types, the average nal gap provided was bellow x. From Figure 3.3, we can observe that all instance types solved by PF as an exact method (we exclude instance types Dfn-bwin and Dfn-gwin) are solved within a gap less than 55%, while the same instances type are solved within a nal gaps up to 75% by the compact MILP formulation.

Sensitivity analysis

For the sensitivity analysis we vary the instances parameters such as: latency value l k , bandwidth b k , node capacities c u and functions capacities m f and observe how the algorithm will behave, by putting l k = α ×l k , b k = α ×b k , c u = α ×c u and m f = α ×m f with α ∈ {1.0, 1.5}. The aim of these tests is to show the parameters that aect directly the solution, and so the costs. We compare the values of the best solutions found within one hour (averaged over 10 graphs per instance type) for each of these settings.

Instances used in Figure 3.6 need at most 500 paths to be solved to optimality. For these instances, we varied the number of max_paths in {50, 100, 500}. We observe that adding a small number of paths permits PF to converge quickly to a very good solution for Geant and Ta1 for which the number of required paths exceeds 200 paths. Conversely, for instances Atlanta and Nobel-germany, which required only MILP-based Heuristic 53 and 149 paths, respectively, to be solved to optimality, 50 paths were insucient to nd a high-quality solution. This can be explained by the sparsity of their graphs and the number of demands.

Instances used in Figure 3.7 require more than 500 paths to be solved to optimality. We notice that for instances "Janos-us", the PF formulation was struggling to nd a good quality solution. For Géant and Ta1, the increasing of costs is explained by the fact that Cplex reaches the time limit for all instances. Moreover, when there are too many paths (variables), Cplex needs more time to nd a good solution, while it is faster with a fewer number of variables. We observe that in this case there is a trade-o between the quality of heuristic solution and the size of the underlying formulation.

For all other instances we can see that the number of xed paths does not aect a lot the behavior of PF.

Conclusions

In this chapter, we have proposed a path-based MILP formulation for the VNFPRP.

With this formulation we were able to provide optimal solutions for additional instances and to improve the gaps obtained from the compact formulation. We then derived an MILP-based heuristic that we have compared to the compact MILP formulation, in terms of the CPU time and the quality of the obtained solution. We have also tested the path formulation with dierent values of problem parameters and discussed the dierence in terms of costs and (in)feasibility.

The proposed Path-based formulation has shown its eectiveness for nding highquality solutions. In the next chapter, we will apply a column generation algorithm to the path formulation to solve its linear relaxation and derive a Branch-and-Price algorithm to solve it.

Chapter 4

Extended formulations

In this chapter, we present two column generation approaches to tackle the Virtual Network Functions Placement and Routing problem. First, we propose two extended MILP formulations to model it, for which we then derive branch-and-price algorithms. For each formulation, we present the pricing problem, its resolution method, computation of the Lagrangian bound, and the branching scheme. Both approaches are empirically compared, and experiments are conducted on realistic instances to show their eciency.

Previous studies related to the simultaneous placement of ordered VNFs and routing of trac demands mainly rely on compact MILP formulations, whose computational performance is fairly limited. To enhance the capabilities of nding exact solutions for larger instances of practical relevance, we focus in this part of the thesis, on the extended formulations with an exponential number of variables. Several studies and works were proposed to deal with MILP compact formulations performance limitation.

According to Sadykov and Vanderbeck [START_REF] Sadykov | Column generation for extended formulations[END_REF], working in an extended variable space allows one to develop tight reformulations for mixed-integer programs. Ford and Fulkerson [START_REF] Randolph | A suggested computation for maximal multi-commodity network ows[END_REF] are the rst that use an optimization sub-problem to price out an exponential number of non-basic variables, which allow them to add only basic variables and reduce the model size. Bilde and Krarup [START_REF] Bilde | Sharp lower bounds and ecient algorithms for the simple plant location problem[END_REF] show that the extended facility location reformulation for uncapacitated lot-sizing was integral, this means that solving the linear relaxation of the reformulation can give an integer optimal solution without any branching. Chopra et al. [START_REF] Chopra | An extended formulation of the convex recoloring problem on a tree[END_REF] propose an extended formulation for the Convex Recoloring problem on a tree. They show that the LP-relaxation of the extended formulation provides an integer optimum for all considered problem instances, which indicates that the LP-relaxation of the extended formulation is a very good approximation of the integer polytope. Vanderbeck and Wolsey [START_REF] Vanderbeck | Reformulation and decomposition of integer programs[END_REF] survey the main reformulations based on decomposition methods, such as Lagrangian relaxation, Dantzig-Wolfe, and the resulting branch-and-price algorithms and Benders' reformulation, they also discuss in detail extended formulations.

in this part of the thesis, we propose two extended MIP formulations for the VNF-PRP. In both cases, we consider an extended variable space admitting an exponential number of variables which allows us to develop tighter MIP reformulations.

In the rst reformulation, we separate the VNF placement problem, which is treated at the master level, from the routing problem, which is solved separately for each commodity in the pricing problem. The second extended formulation arises from an alternative Dantzig-Wolfe decomposition approach. The problem is decomposed per commodity in such a way that: the master problem ensures that exactly one path with its associated VNF installations is chosen for each commodity, and that node and VNFs capacity constraints are satised. The routing and VNF placement constraints associated for each commodity are managed in the pricing problems.

In order to improve the LP-relaxation bounds of our formulations, new families of valid inequalities are also proposed. All these elements are combined in two ecient Branch-and-Price algorithms. A detailed computational comparison of the Branchand-Price algorithms against the compact MIP formulation and the automatic Benders decomposition approach by the commercial solver Cplex is given.

Outline of the chapter The chapter is organized as follows. In Sections 4.1 and 4.2, we present the two extended formulations. We discuss the associated pricing problems, detail the computation of the Lagrangian bound and present branching schemes.

In Section 4.3, we provide a new set of valid inequalities that aim to strengthen the LP-relaxation bounds. In Section 4.4 we provide implementation details of Branchand-Price (B&P) algorithms for each formulation. In Section 4.5, we provide a theoretical result proving that the Dantzig-Wolfe formulation provides better LP-relaxation bounds than the path formulation. In Section 4.6, we discuss the obtained numerical results, and conclude with some remarks and perspectives in Section 4.7.

First extended formulation: the model PF

In this section, we present the rst extended MILP formulation, denoted by PF (which stands for path-based formulation) to model the VNFPRP. The formulation has been introduced in Chapter 3 to generate heuristic solutions, by considering a compact model obtained from choosing a small but promising number of columns. In this chapter instead we focus on developing an exact method for solving the path formulation, based on a Branch-and-Price procedure. In this model we use latency-constrained elementary path variables associated to each commodity to model routing decisions.

In this section we discuss theoretical properties of this model, along with a derivation of a valid Lagrangian bound, whereas the details related to the B&P implementation are given in Section 4.4. The master problem aims to nd the optimal VNFs installation for a given routing path for each trac request. One pricing problem is dened per commodity k ∈ C, and it consists on determining an s k -d k latency-constrained elementary routing path.

Decision variables

The set of variables required for the path formulation is already dened in Chapter 3, Table 3.1. Similarly, the constraints modeling the path formulation are presented in Chapter 3. We recall the path formulation variables and MILP model in order to derive the associated dual formulation.

Let us denote by P k the set of all latency-constrained elementary paths associated with commodity k ∈ C. We assume that the set P k is given and that, for each chosen path, the arcs composing it are known. Let t pk uv be the parameter that is equal to 1 if arc (u, v) belongs to path p, p ∈ P k , and equal to 0 otherwise.

The master problem formulation

The VNFPRP can be modeled as follows:

(P F ) : min

u∈N f ∈F ψ f u z f u + u∈N ψ u w u (4.1) p∈P k λ k p = 1 k ∈ C (α k ) (4.2) f ∈F z f u ≤ c u w u u ∈ N (β u ) (4.3) k∈C y f k u b k ≤ m f z f u f ∈ F, u ∈ N (γ f u ) (4.4) y f k u + y gk u ≤ 1 k ∈ C, (f, g) ∈ A k , u ∈ N (δ (f,g)k u ) (4.5) ( p∈P k (u,v)∈p t pk uv λ k p -1) + (x f k v -x f k u ) ≤ y f k v k ∈ C, f ∈ F k , (u, v) ∈ A (η f k uv ) (4.6) 100
Extended formulations ) admits an exponential number of path variables; thus, a column generation (CG) procedure is needed to solve its continuous relaxation. In the following, we describe the pricing problem and discuss three possible procedures for its resolution.

x gk u ≤ x f k u k ∈ C, f, g ∈ F k : f ≺ k g, u ∈ N (ϕ (f,g)k u ) (4.7) y f k u ≤ x f k u k ∈ C, f ∈ F k , u ∈ N (θ f k u ) (4.8) y f k u ≤ (v,u)∈A p∈P k t pk vu λ k p k ∈ C, f ∈ F k , u ∈ N (π f k u ) (4.9) u∈N y f k u ≥ 1 k ∈ C, f ∈ F k (ϑ f k ) (4.10) x f k u = 0, u = s k 1, u = d k k ∈ C, f ∈ F k (σ f k d k ) (4.
The LP-relaxation of this model is obtained by replacing constraints (4.12)-(4.13) with λ ≥ 0, 0 ≤ x ≤ 1, y ≥ 0, 0 ≤ w ≤ 1 and z ≥ 0. The dual of the master problem is given by the following linear program:

(DP F ) : max k∈C α k - k∈C u∈N (f,g)∈S k δ (f,g)k u - k∈C f ∈F k (u,v)∈A η f k uv + k∈C f ∈F k ϑ f k + k∈C f ∈F k σ f k d k - k∈C f ∈F k u∈N ∆ f k u - u∈N Ω u α k - (u,v)∈A f ∈F k t pk uv η f k uv + u∈N v∈Γ -(u) f ∈F k t pk vu π f k u ≤ 0 k ∈ C, p ∈ P k c u β u -Ω u ≤ ψ u u ∈ N -β u + m f γ f u ≤ ψ f u f ∈ F, u ∈ N -b k γ f u - g∈F k (f,g)∈S k δ (f,g)k u - g∈F k (g,f )∈S k δ (g,f )k u v∈Γ + (u) η f k uv -θ f k u -π f k u + ϑ f k ≤ 0 k ∈ C, f ∈ F k , u ∈ N - v∈Γ -(u) η f k vu + v∈Γ + (u) η f k uv - g∈F k (f,g)∈S k ϕ (f,g)k u + g∈F k (g,f )∈A k ϕ (g,f )k u + θ f k u -∆ f k u ≤ 0 k ∈ C, f ∈ F k , u ∈ N \ {s k , d k } - v∈Γ -(s k ) η f k vs k + v∈Γ + (s k ) η f k s k v - g∈F k (f,g)∈S k ϕ (f,g)k u + g∈F k (g,f )∈A k ϕ (g,f )k s k + θ f k s k + σ f k s k -∆ f k s k ≤ 0 k ∈ C, f ∈ F k 102 Extended formulations - v∈Γ -(d k ) η f k vd k + v∈Γ + (d k ) η f k d k v - g∈F k (f,g)∈S k ϕ (f,g)k d k + g∈F k (g,f )∈A k ϕ (g,f )k d k + θ f k d k + σ f k d k -∆ f k d k ≤ 0 k ∈ C, f ∈ F k β, γ, δ, η, ϕ, θ, π, ϑ, σ, ∆, Ω ≥ 0, α k ∈ R |C|
Then, for each k ∈ C and p ∈ P k , the dual constraint associated with the path variable λ k p is given as follows:

α k - f ∈F k [ (u,v)∈A t pk uv η f k uv + u∈N v∈Γ -(u) t pk vu π f k u ] ≤ 0, k ∈ C, p ∈ P k (4.14) ⇐⇒ α k + (u,v)∈A f ∈F k (π f k v -η f k uv )t pk uv ≤ 0, k ∈ C, p ∈ P k (4.15)

The pricing problem

As customary in column generation, the master problem is initialized with a subset of λ variables (resulting in the so-called restricted master problem ), and then the additional variables necessary to solve the LP-relaxation of the model are generated on the y by separating the associated dual constraints (4.15). The pricing problem then consists of nding for each commodity k, a path p ∈ P k with negative reduced costs, i.e., a path p such that:

α * k + (u,v)∈p f ∈F k (π * f k v -η * f k uv ) > 0, (4.16)
where (α * , η * , π * ) refers to a sub-vector of an optimal dual solution of the restricted master problem. This dual solution will be used for dening the pricing problems and computing the Lagrangian bound (cf. Section 4.1.5).

Thus, for each commodity k ∈ C, a separate pricing problem is dened in order to nd an s k -d k latency-constrained elementary path of minimum cost. Cost per each arc is dened as

cuv = f ∈F k (η * f k uv -π * f k v ), (u, v) ∈ A.
(4.17 4.1 First extended formulation: the model PF 103

If we nd a path p ∈ P k such that (u,v)∈p cuv -α * k < 0, the associated variable λ k p will be inserted in the restricted master problem. Based on the dual solution, the values dened by (4.17) may be negative. Therefore the pricing problem consists of nding an elementary shortest path satisfying latency constraints on a graph that may contain negative cycles. Hence, the pricing problem is strongly NP-hard [START_REF] Dror | Note on the complexity of the shortest path models for column generation in VRPTW[END_REF]. Three dierent methods are proposed and computationally evaluated for solving this pricing problem (cf. Sections 4.4.3 and 4.6, respectively).

Lagrangian bound

In the following, we describe how to derive the Lagrangian bound of the model PF which can be used to minimize the number of iterations in the column generation procedure. Let DRPF denote the LP-dual of the restricted master problem of the PF model, let Υ * be its optimal solution, Z DRP F the associated objective value, and let DPF denote the dual of the model PF (including all columns). We will slightly abuse the notation and focus only on the (α, π, η) components of the vector Υ.

Υ * satises all constraints of DRPF, but not necessarily all constraints of DPF, i.e., Υ * is optimal for DRPF, but not necessarily feasible for DPF, as the constraints associated with the columns that are left out from the master are missing.

In what follows, we will show how to construct a feasible solution for DPF during the column generation phase. Let Z * k be dened as:

Z * k = min p∈P k [-α * k - (u,v)∈p f ∈F k (π * f v -η * f uv )] = -α * k -max p∈P k (u,v)∈p f ∈F k (π * f v -η * f uv ).
So, Z * k is the minimum reduced cost associated with commodity k (over all paths p ∈ P k ). Therefore, for each commodity k ∈ C and each path p ∈ P k , it holds:

Z * k ≤ -α * k - (u,v)∈p f ∈F k (π * f v -η * f uv )
Consequently, constraints (4.15) in the DPF can be written as:

Z * k + α * k + (u,v)∈p f ∈F k (π * f v -η * f uv ) ≤ 0 k ∈ C, p ∈ P k . (4.18)
The feasible solution Υ D for the DPF is build based on Υ * , and by using the following variable change: ᾱk = α * k + Z * k , k ∈ C. It follows that for all k ∈ C and for all p ∈ P k :

(4.18) ⇐⇒ ᾱk + (u,v)∈p f ∈F k (π * f v -η * f uv ) ≤ 0
Thus, Υ D = (ᾱ * , η * , π * ) is a feasible solution for DPF. Correspondingly, LB = Z DRP F + k∈C Z * k is the objective value of the DPF dening a valid lower bound for the model PF.

The path formulation can be obtained by applying Dantzig-Wolfe decomposition on the compact formulation proposed in Chapter 2. The pricing problem could generate cyclic graphs (columns). Therefore, these columns will not appear in any feasible solution of the master problem because of precedence constraints. In the path formulation, considered in this thesis, we have added sub-tour elimination cuts in the pricing problem in such a way that it generates only elementary paths. The later formulation cannot be obtained from Dantzig-Wolfe decomposition of the compact model given in Chapter 2. This is due to the fact that sub-tour elimination cuts are needed in the ILP pricing problem.

However, the model PF can be seen as a Dantzig-Wolfe problem reformulation obtained from the compact model studied in Chapter 3 after enhancing the latter with subtour elimination constraints imposed on the ow variables. The lower bound derived in this section corresponds to the Lagrangian bound of this Dantzig-Wolfe problem reformulation (we leave out further details and refer the reader to e.g., [START_REF] Marco | Selected topics in column generation[END_REF], for the general theory of Dantzig-Wolfe decomposition). Corollary 3.3 dened in Chapter 3 indicates that, when implementing a B&P procedure to solve the model PF, we can sidestep branching on the exponential variables λ and apply the regular branching scheme dened only for x, y, z, and w variables. In particular, this means that (apart from the change of the coecients in the objective function), the pricing problem will not be aected by branching decisions.

Second extended formulation: the model DW

In this section, we present an alternative extended formulation for the VNFPRP, to which we refer as the model DW (where DW stands for Dantzig-Wolfe). First, we introduce the master problem and then, based on the master problem's dual solution, 4.2 Second extended formulation: the model DW 105 we describe the pricing problem and show how to calculate the value of the Lagrangian bound. At the end of this section, a branching scheme is proposed.

In the following, we use the term path-installation to indicate a path p with preinstalled VNFs satisfying latency, conict, and precedence constraints. The master problem aims to choose one path-installation per each commodity k while respecting node and function capacity constraints.

Decision variables

Let us denote by T k the set of all path-installations associated with commodity k. To create our master problem, we need three families of variables, as described in Table 4 The set T k associated with each commodity k containing all feasible path-installations is supposed to be known. Thus, the placement of each VNF for each path-installation at network nodes is uniquely dened. Let us denote by a f pk u the parameter that is equal to 1 if the VNF f is used at node u for the path-installation p associated with commodity k; and that is equal to 0 otherwise.

a f pk u =        1, if VNF
f is used at node u for the path-installation p associated with commodity k, 0, otherwise.

ILP formulation

The model DW is then given as: [START_REF] Monawar | Anti-theft vehicle tracking and regaining system with automatic police notifying using haversine formula[END_REF] Extended formulations

(DW ) : min u∈N f ∈F ψ f u z f u + u∈N ψ u w u (4.19) p∈T k τ k p = 1 k ∈ C (α k ) (4.20) f ∈F z f u ≤ c u w u u ∈ N (β u ) (4.21) k∈C p∈T k a f pk u τ k p b k ≤ m f z f u f ∈ F, u ∈ N (γ f u ) (4.22) τ k p ∈ {0, 1} k ∈ C, p ∈ T k (4.23) w u ∈ {0, 1} u ∈ N (η u ) (4.24) z f u ∈ N u ∈ N, f ∈ F

The dual of the master problem

The dual of the master problem is given by the following linear program:

max k∈C α k - u∈N η u α k - u∈N f ∈F k a f pk u b k γ f u ≤ 0 k ∈ C, p ∈ T k c u β u -η u ≤ ψ u u ∈ N -β u + m f γ f u ≤ ψ f u f ∈ F, u ∈ N β, γ, η ≥ 0, α ∈ R |C|
α k - u∈N f ∈F k a f pk u b k γ f u ≤ 0 k ∈ C, p ∈ T k , ( 4
α * k - u∈N f ∈F k a f pk u b k γ * f u > 0.

The pricing problem

For each commodity k, we have one pricing problem that aims to nd a path-installation.

The left-hand-side in inequalities (4.26) characterizes the objective function of the pricing problem. The set of variables required in the pricing problem is described in Table 4.3.

Variables Type Proof. Observe that variables n do not appear in the objective function, and that the location variables h (which basically determine the value of the solution) remain binary. Υ * is optimal for DRDW but not necessarily feasible for the DDW; this means that there exists at least one constraint (4.26) in DDW which is violated by this solution.

d f u 1, if
max α * k - u∈N f ∈F k h f u b k γ * f u 108 Extended formulations (u,v)∈A n uv - (v,u)∈A n vu =        -1 if u = d k , 1 if u = s k , 0 otherwise. u ∈ N (4.27a) (u,v)∈A n uv l uv ≤ l k (4.27b) h f u + h g u ≤ 1 (f, g) ∈ A k , u ∈ N (4.27c) (n uv -1) + (d f v -d f u ) ≤ h f v f ∈ F k , (u, v) ∈ A (4.27d) d g u ≤ d f u f, g ∈ F k : f ≺ k g, u ∈ N (4.27e) h f u ≤ d f u f ∈ F k , u ∈ N (4.27f ) h f u ≤ (v,u)∈A n vu f ∈ F k , u ∈ N (4.27g) u∈N h f u ≥ 1 f ∈ F k (4.27h) d f s k = 0 f ∈ F k (4.27i) d f d k = 1 f ∈ F k (4.27j) (d,
For a xed k ∈ C we have:

Z * k = -α * k + min p∈T k { u∈N f ∈F k a f pk u b k γ * f u }
As Z * k represents the minimum reduced cost (over all path installations p ∈ T k ), the dual constraint (4.26) in DDW can be written as follows:

Z * k ≤ -α * k + u∈N f ∈F k a f pk u b k γ * f u ⇐⇒ Z * k + α * k - u∈N f ∈F k a f pk u b k γ * f u ≤ 0

Branching on τ variables

The LP-relaxation of the Dantzig-Wolfe formulation solved by column generation procedure is not necessarily integral. Furthermore, applying the Branch-and-Bound algorithm on the restricted master problem with only the generated columns at the root node will not guarantee a feasible solution and so an optimal solution. Moreover, at each branching node, there may exist new columns with a negative reduced cost which 4.2 Second extended formulation: the model DW 111 should be added to the master problem. Therefore, in order to nd an optimal integer solution, we should generate columns at each branching node.

Various branching schemes, specic (like the one proposed below) or generic (see e.g., [START_REF] Vanderbeck | Branching in Branch-and-Price: a generic scheme[END_REF]), can be used to generate integer solutions using column generation procedure embedded within the Branch-and-Bound algorithm. The resulting algorithm is called Branch-and-Price.

In the following, a commodity k ∈ C is called fractional if it admits a fractional τ variable. For a path-installation p ∈ T k , we use the notation u ∈ p to indicate that the path-installation p passes through the node u. Furthermore, the notation a f pk u = 1 is used to indicate that the path-installation p passes through the node u on which the VNF f is installed for commodity k.

0 < p∈T k a f pk u =1 τ k p < 1 (4.30)
Proof. Let us suppose that k is a fractional commodity but neither Case 1 nor Case 2 holds. Hence, we have:

p∈T k u∈p,v∈p τ k p ∈ {0, 1} u, v ∈ N (4.31) p∈T k a f pk u =1 τ k p ∈ {0, 1} u ∈ N, f ∈ F k . (4.32)
From (4.32) we can distinguish two cases: (a) there exists a node u ∈ N and a function f ∈ F k such that

p∈T k a f pk u =1 τ k p = 1, or (b) p∈T k a f pk u =1
τ k p = 0, for each u ∈ N and f ∈ F k .

(a) By constraints (4.20), and because k is a fractional commodity, all path-installations of T k in the solution are fractional. Let p 1 ∈ T k be a fractional path-installation passing through node u on which VNF f is installed (i.e., 0 < τ k p 1 < 1). As

p∈T k a f pk u =1 τ k p =
p∈T k a gpk v =1 τ k p ∈ {0, 1} must hold. Given that τ k p 1 > 0, then (1) p∈T k a gpk v =1 τ k p > 0. Thus p∈T k a gpk v =1 τ k p = 0 cannot hold. Accordingly, p∈T k a gpk v =1
τ k p = 1. We know that 0 < τ k p 2 < 1 and that g is not installed on v for p 2 , we will have: ( 2)

p∈T k a gpk v =1 τ k p < 1. Therefore, (1) and
(2) contradict hypothesis (4.32).

(I) p 1 and p 2 pass through at least one dierent node, i.e., there should exist another node v = u belonging to p 2 and not to p 1 . We notice that VNF f is installed only on node u for both p 1 and p 2 , and that another VNF g = f can or not be installed on node v belonging to the path-installation p 2 . Since p 2 contains u and v and τ k p 2 > 0, this implies: (i)

p∈T k u∈p,v∈p τ k p > 0.
Moreover, as p∈T k a f pk u =1

τ k p = 1 and τ k p 1 > 0 and we know that path-installation p 1 does not pass through node v, then the value τ k p 1 can be deleted from the following sum: τ k p = 0).

p∈T k u∈p,v∈p τ k p , which implies that (ii) p∈T k u∈p,v∈p τ k p < 1.
Therefore, the result holds. A branching scheme is said to be complete, if it can generate any feasible solution.

From Proposition 4.2 we conclude that our branching scheme proposed for the model DW is complete, as at least one of the two cases should hold for any fractional commodity.

Strengthening inequalities

In this section, we derive several families of valid inequalities that can strengthen the LP-bounds of both proposed extended formulations. We rst present inequalities that can be used to directly enhance the model PF. We then present inequalities that are valid for both models and that can be exploited in case a function's capacity is smaller than the respective trac demand. We close this section by explaining how some of inequalities proposed for the model PF can be used within the Dantzig-Wolfe decomposition to strengthen the model DW. Given a commodity k ∈ C and a node u ∈ N , if there is a unique path p going from s k to u in G, let A p be the arcs belonging to the path p.

y f k u ≤ w u , u ∈ N, k ∈ C, f ∈ F k . ( 4
y f k u + y gk u ≤ w u , u ∈ N, k ∈ C, (f, g) ∈ A k .
f ∈Q y f k u ≤ 1, u ∈ N, k ∈ C, Q ∈ D k . ( 4
y f k u ≤ w u , u ∈ N, k ∈ C, Q ∈ D k .
Proposition 4.7 For a given commodity k ∈ C, and a node u ∈ N , if there exists a unique path from s k to u in G, then inequalities (4.37) are valid for the VNFPRP:

f ∈Q x f k u ≤ |A p |, u ∈ N, k ∈ C, Q ∈ D k : |A p | < |Q|. (4.37)
Proof. The proof is given for a xed commodity k and is valid for all commodities.

Let u be the node for which we have a unique path p going from s k to u, and let Q be a clique in the graph D k , such that |A p | < Q. The number of VNFs in conict installed at or before node u should be less than or equal the number of arcs in the path p; otherwise, two or more functions in conict need to be installed at the same node, which leads to an infeasible solution. Proof. Recall that each commodity requires at least one VNF. As all required VNFs should be installed at graph nodes, at least one node in the graph must be activated.

Furthermore, if there is a conict between VNFs associated with one commodity k, then the number of activated nodes should be at least equal to the maximum number of VNFs in conict, which is the clique number of D k .

Proposition 4.9 Inequalities (4.39) are valid for the VNFPRP.

u∈N f ∈Q y f k u ≥ |Q|, k ∈ C, Q ∈ D k .

Extended formulations

Let C u ⊆ C be a subset of commodities for which there exists at least one latencyconstrained path visiting node u, (i.e., if k / ∈ C u , this means that all paths associated with k do not enter the node u). Let N k be the set of nodes belonging to at least one latency-constrained path associated with commodity k ∈ C (this can be checked in polynomial time using a min-cost ow algorithm for example).

Proposition 4.10 The following inequalities are valid for the VNFPRP.

k∈C\C u f ∈F k y f k u = 0, u ∈ N, (4.40
)

k∈C\C u f ∈F k x f k u = 0, u ∈ N. (4.41)
Moreover, inequalities (4.42) are valid and dominate inequalities (4.10); Proposition 4.12 Inequalities (4.46) are valid for the VNFPRP: x f k u ≤ 1 -

u∈N k y f k u ≥ 1, k ∈ C, f ∈ F k , (4.42) f ∈F z f u = 0, u ∈ N \ {∪ k∈C N k }, (4.43
z f u ≥ b k m f y f k u , k ∈ C, f ∈ F k , u ∈ N.
x f k u ≤ 1 -y f k d k , k ∈ C, u ∈ Γ -(d k ), f ∈ F k , ( 4
v∈N d k y f k v , k ∈ C, u ∈ N s k , f ∈ F k (4.47)
Proof.

x f k u ≤ 1 - v∈N d k y f k v + v∈S k x f k v , k ∈ C, u ∈ N s k , f ∈ F k (4.48)
Proof.

c u > k∈C f ∈F k b k m f .
1) Therefore, inequalities (4.49) are valid for the VNFPRP and dominate inequalities (4.3). 3) In addition if there exists a conict between functions in F k for a given commodity

f ∈F z f u ≤ k∈C f ∈F k b k m f w u , u ∈ N. ( 4
f ∈F z f u ≤ k∈C u f ∈F k b k m f w u , u ∈ N.
k ∈ C, with m f 1 ≤ m f 2 ≤ • • • ≤ m f |Q| , Q ∈ D k and c u ≥ k∈C f ∈F k b k m f , then
inequalities (4.50) are dominated by the following inequalities.

f ∈F z f u ≤ k∈C Q∈D k [( f ∈F k b k m f ) - |Q| i=2 f i ∈Q b k m f i ]w u , u ∈ N. (4.51)
Proof.

1) If there exists a node u having enough capacity to install VNFs required for all commodities in C, then the number of functions needed to treat all commodities bandwidth is bounded by

k∈C f ∈F k b k m f
, which represents the maximum number of VNFs necessary to handle all demands.

2) Only VNFs associated with commodities having at least one path passing through a node u can be installed at node u.

3) The number of VNFs installed at node u when the conict constraints are considered is bounded by the maximum number of copies needed to install the VNFs with the smallest capacity for each commodity (i.e., in the worst case we will keep the VNFs with the maximum instantiation installed at node u and install other VNFs at the other nodes).

p∈T k a f pk u τ k p ≤ w u , k ∈ C, u ∈ N, f ∈ F k (4.52) p∈T k (a f pk u + a gpk u )τ k p ≤ w u , k ∈ C, u ∈ N, (f, g) ∈ A k (4.53) z f u ≥ b k m f p∈T k a f pk u τ k p , k ∈ C, f ∈ F k , u ∈ N (4.54)
Adding these valid inequalities generates new (non-negative) dual variables in the dual of the master program, that we denote by δ, η and ϕ, respectively. Thus, dual constraints (4.26) need to be replaced by inequalities (4.55)

α k - u∈N   f ∈F k (b k γ f u + b k m f ϕ f k u + δ f k u )a f pk u + (f,g)∈A k η f gk u (a f pk u + a gpk u )   ≤ 0, k ∈ C, p ∈ P k (4.55)

Branch-and-Price algorithms

In this section, we present two B&P algorithms implemented for the models PF and DW, respectively.

Generic column generation framework

Initialization The restricted master problem of both models is initialized by a subset of columns building a heuristic solution which is obtained in the initialization phase of the algorithm (see Section 4.4.4). If no solution has been found during a time-limit, the CG framework is initialized with an articial column whose cost is set to a very large number.

Bounding At each iteration of the column generation procedure, the restricted master problem is solved, and a dual solution is provided. Accordingly, the objective function of the pricing problem for each commodity k is updated, and the pricing problem is solved. Depending on the pricing strategy (see Section 4.4.3) multiple columns per commodity having negative reduced costs (or at most one) are added to the restricted master problem. During this process, to reduce the number of CG iterations, we keep track of the Lagrangian bound and compare it to the objective value of the current restricted master problem. If the dierence between the two values is smaller than ε, the column generation procedure is stopped and we resort to branching.

Branching

At the end of the column generation phase, the integrality of the solution of the relaxed master problem is veried. If the current solution is not integer, we branch on the most fractional variable, applying the BFS (Breadth-First Search) based branching strategy. Specically, we explore all the nodes of the same level in the branching tree before moving to the next level. In our implementation the algorithm explores all nodes admitting a fractional feasible solution of the same level by applying the respective branching scheme described below. For each branching node with a fractional solution, two children nodes are created and saved in a queue. The nodes in the queue are explored using the FIFO (First In First Out) method. A global lower bound is calculated at each level. In our preliminary experiments, we also tried the diving strategy as an alternative to the BFS-based branching. Whereas diving is very useful when searching for feasible solutions (see e.g., [START_REF] Furini | A column generation heuristic for the two-dimensional two-staged guillotine cutting stock problem with multiple stock size[END_REF][START_REF] Gérard | Column generation based approaches for a tour scheduling problem with a multi-skill heterogeneous workforce[END_REF][START_REF] Joncour | Column generation based primal heuristics[END_REF]), in our case this strategy did not prove useful, because a high-quality feasible solution is used to initialize the CG procedure (see above).

Branching scheme for the model PF In Corollary 3.3 we showed that the binary constraints on λ variables in the PF model can be relaxed to λ ≥ 0. Hence, in our BP implementation of the PF model, we branch only on the (x, y, z, w) variables.

Dierent branching schemes have been tested for the model PF; The one outperforming all others is to branch rst on the most fractional w variables (by imposing either w ≥ 1 or w ≤ 0), secondly on z by setting either z ≥ z * or z ≤ z * , thirdly on y which are forced to be y ≥ 1 or y ≤ 0 and nally on x variables using x ≥ 1 or x ≤ 0.

Branch-and-Price algorithms 121

Branching scheme for the model DW In our branching scheme for the model DW we start branching on z variables by setting z ≥ z , or z ≤ z and then on w variables by setting w ≤ 0, or w ≥ 1. When z and w variables are integers, we continue by branching on the τ variables. Given that the path-installation variables are generated as and when they are needed, we follow the specic branching scheme for τ variables proposed in Proposition 4.2. Specically, if we nd a pair of two distinct nodes u, v ∈ N such that

p∈T k :u∈p,v∈p τ k
p is fractional we create two branches by imposing constraints that limit the respective sum to 0, or 1, respectively. If none such pair can be found, we search for a node u ∈ N and a function f ∈ F k such that

p∈T k :a f pk u =1 τ k p is
fractional, and create two branches correspondingly.

Extended formulations

Pricing strategy

In the sequel, we rst propose three dierent pricing strategies for the model PF.

These strategies are later computationally evaluated in Section 4.6. For the model DW, the underlying pricing problem exhibits a rich combinatorial structure. We therefore employ an o-the-shelf MIP solver to price in the columns. 1) We utilize dynamic programming; 2) Based on Yen's algorithm [START_REF]Finding the k shortest loopless paths in a network[END_REF], we derive a reduced cost method, and use it in two dierent ways; and 3) We model the pricing problem as an ILP that we solve using an o-the-shelf solver.

Dynamic programming Let k ∈ C be a given commodity and let, for each arc (u, v) ∈ A, cuv given by (4.17) be the corresponding arc cost. We denote by F v (S, l v )

the minimal cost of the partial path going from source node s k to node v visiting all nodes in S exactly once and ready to leave node v with latency value equal to l v . The following equations (see also [START_REF] Righini | New dynamic programming algorithms for the resource constrained elementary shortest path problem[END_REF]) illustrate the recursive function used for the pricing problem associated with commodity k:

DP :    F s k (∅, 0) = 0 F v (S, l v ) = min v∈Γ + (u) {F u (S -{v}, l u ) + cuv | (u ,v )∈A S l u v ≤ l k }
where S is the set of nodes in the path, A S is the set of arcs connecting nodes in S, l u (resp. l v ) represents the sum of the latency of arcs in A S from source node s k until node u (resp. v). We utilize the dynamic programming algorithm proposed in [START_REF] Boland | Accelerated label setting algorithms for the elementary resource constrained shortest path problem[END_REF] (cf.

the set P k . However, multiple columns with negative reduced cost can be added for each commodity. We consider two settings: In the rst one (that we refer to as Red Cost 1), we add at most 10 columns with negative reduced cost per commodity; In the second one (that we refer to as Red Cost 2), all columns from P k with negative reduced cost are added in each iteration. We tested both approaches, and the results are shown in Section 4.6.

At the end of the column generation phase, and before resorting to branching, we also try to x some columns based on their reduced cost. Using an upper bound U B (obtained from the incumbent solution), and the value Z of the current restricted master problem, if there exists a commodity k ∈ C and a path p ∈ P k such that RC(p) + Z > U B, this means that the variable λ k p can be xed to zero. Therefore, in relaxed pricing problem we introduce binary variables t uv which are set to one if arc (u, v) ∈ A belongs to the path associated with the current commodity; and to 0, otherwise. The RPP associated with commodity k ∈ C is equivalent to the following integer linear program, in which constraints (4.56) ensure that the ow goes from the source node s k to the destination node d k , whereas (4.57) represents the latency constraint:

(RP P ) : min v∈Γ + (u) t uv ≤ 1, which aims to force the path to pass only one time through a node u, with degree(u) ≥ 3.

(u,v)∈A cuv t uv (u,v)∈A t uv - (v,u)∈A t vu =        -1 if u = d k , 1 if u = s k , 0 otherwise. u ∈ N (4.56) (u,v)∈A t uv l uv ≤ l k (4.57) t uv ∈ {0, 1} (u, v) ∈ A Let C = (N C , A C ) be
Disconnected paths like the one shown in Figure 4.5(b) are eliminated by the inequality

(u,v)∈A C t uv ≤ |N C | -1.
We point out that the eciency of this procedure highly depends on the underlying cost function c. As we will see in Section 4.6, when there are very few negative cycles in the graph, this procedure may outperform the other alternatives, including dynamic programming. In order to price the columns associated with the DW formulation, the MIP model presented in Section 4.2.3, is solved for each commodity k ∈ C, using an o-the-shelf solver. Thus, at each iteration of the column generation procedure at most one column per commodity with negative reduced cost is generated by the pricing problems. As in the case of the model PF, the column generation phase terminates when no more columns with a negative reduced cost are found, or when the gap between the current value of the restricted master problem and the Lagrangian bound is smaller than ε.

Heuristics

Before entering the B&P phase, we generate a heuristic solution which provides a high-quality upper bound and a promising set of columns that we use to initialize the CG procedure. For the model PF, we employ the MIP-based heuristic presented in Chapter 3. The heuristic solves a compact model derived from the formulation PF in which only a small subset of latency-constrained paths is considered. To obtain this subset, we run Yen's algorithm [START_REF]Finding the k shortest loopless paths in a network[END_REF] which provides κ-elementary paths between two nodes, sorted from the shortest to the longest one. The number of generated paths per commodity is capped by κ, and we let κ ∈ {10, 15, 20, . . . , 50}. As soon as the underlying MIP-based heuristic nds a feasible solution for a xed κ value, we stop. However, for some instances, even for κ = 50 we fail to nd a feasible solution.

For the model DW, we start with an articial column, price-in the columns with negative reduced cost, and then convert the obtained linear program into a MIP.

In both cases, there is a time limit after which this MIP-based heuristic initialization is aborted.

Comparing the LP-relaxations

In this section, we compare the LP-relaxation values of the path formulation (PF) and Dantzig-Wolfe (DW) formulation. First, a graphic that shows that the LP-relaxation bounds provided by the DW formulation are better than those given by the PF is presented, followed by the LP-relaxation comparison proof. Proof. We prove that from any feasible solution of DW, we can construct a feasible solution for PF. Let X DW = (τ DW , z DW , w DW ) be the feasible solution of the relaxed DW associated with Z DW . From the solution X DW , a feasible solution X P F = (λ P F , x P F , y P F , z P F , w P F ) for PF is built as follows :

Variables τ DW in Dantzig-Wolfe formulation represent the path-installation chosen for each commodity k ∈ C. From τ DW variables, the information concerning the arcs belonging to each path and the associated VNFs installation and their order are extracted. Therefore, from the values of τ DW in the solution X DW , the value of λ P F variables in path formulation can be obtained. Furthermore, the set of arcs belonging to each path are provided.

Moreover, in order to get the value of variables x P F and y P F in the path formulation we set x P F = p∈T k τ DW d DW and y P F = p∈T k τ DW h DW , with d DW , h DW ∈ {0, 1} are parameters extracted from the variables τ DW (obtained from the pricing problem associated with commodity k in DW formulation).

Computational results 131

Finally, Variables z DW and w DW are exactly the same for both formulations, therefore, their values are obtained by setting z P F = z DW and w P F = w DW Now, we show that the solution X P F built from X DW is feasible for the path formulation, i.e., satises all the constraints in the PF. From the DW formulation we can see clearly that path constraints (4.2), node-capacity constraints (4.3) and VNF-capacity constraints (4.4) are satised; otherwise d DW is not feasible for DW. Furthermore, each path-installation in DW formulation represents an s k -d k elementary path satisfying ow, latency, conict, installation and precedence constraints. Conict constraints (4.5) are satised by construction, we know from path constraints (4.20) that the sum of τ DW variables is equal to 1, as the parameter h DW ∈ {0, 1}, the conict constraints are satised in DW formulation, we conclude that variables y P F satisfy the conict constraints in the PF. Similarly, as the parameter d DW ∈ {0, 1} and from (4.20), we can see clearly that precedence constraints (4.6), (4.7) and (4.8) are satised. Each pathinstallation passes through the VNFs installations, which means that constraints (4.9) are satised. Each path-installation guarantees that all VNFs associated with a given commodity are installed at nodes, accordingly, and by transformation, y P F variables in the path formulation satisfy the installation constraints (4.10). Therefore, from any solution of the Dantzig-Wolfe formulation, a feasible solution for the path formulation can be derived with Z P F ≤ Z DW . We refer the reader to our computational study, where we report on the instances from which the LP-bound Z DW is strictly stronger than Z P F .

Computational results

In this section we analyze the scalability and the eciency of the two proposed B&P algorithms and show the benets of the valid inequalities dened in Section 4.3. The B&P algorithms are compared to two other exact methods using the commercial solver CPLEX: the rst one is a compact MIP formulation (denoted by C) proposed in Chapter 2 and the second one is the Automatic Benders approach by CPLEX [START_REF] Bonami | Implementing Automatic Benders Decomposition in a Modern MIP Solver[END_REF] applied to the model C in which a family of ow variables is linearly relaxed. Our experiments are designed so as to evaluate the eectiveness and the performance of the proposed extended formulations in terms of CPU time, quality of bounds and nal gaps. Eventually we also measure the advantage of the proposed valid inequalities in improving the LP-bounds and reducing the nal gaps.

All the experiments described in this section were made using a computer with Intel(R)Xeon(R) CPU E5-2650 v2 processor clocked at 2.60GHz, 32 cores, 2 threads per core and 252GB RAM, under Linux operating system. All methods are implemented using the Python API for CPLEX, which is run in single-thread mode with a default memory limited to 20GB. All CPLEX parameters were set to their default values. A default time limit of one hour is set for each tested instance. For the initialization heuristic used within the B&P algorithm (cf. Section 4.4.4), the time limit is xed to AB: The Automatic Benders decomposition available in Cplex [START_REF] Bonami | Implementing Automatic Benders Decomposition in a Modern MIP Solver[END_REF], applied to the setting C in which binary ow variables are linearly relaxed.

Computational results

133

Benchmark instances

In order to perform our experiments, we have generated a set of instances derived from the SNDlib library [125] of telecommunication networks. The instances creation is detailed in Chapter 3. The set of instance-type used to test our algorithms is dened as follows: {Abilene, Atlanta, Di-yuan, France, Geant, Newyork, Nobel-eu, Nobel-germany, Nobel-us, Pdh, Polska, Ta1}.

Obtained results

In the sequel we summarize the major results obtained by our computational study. Tables with more detailed information provided per each instance can be found in Section 4.6.2. This subsection is divided into two parts. The rst part is devoted to comparing the LP-relaxation bounds generated by applying the column generation algorithm to the models PF and DW, respectively. We start by comparing four dierent pricing methods proposed for the model PF (namely, DP, ILP, Red Cost 1, Red Cost 2) in terms of CPU time, number of added columns and generated iterations. Next, after determining the best conguration (i.e., the best pricing method) for the model PF, we focus on the improvement of the LP-gap, with respect to the LP-relaxation bounds provided by the compact formulation (C). We compare the two extended formulations, with and without adding valid inequalities.

In the second part of this subsection, we compare the two proposed B&P algorithms (PF and DW) with two others alternative MIP approaches, namely the Compact formulation (C) and the Automatic Benders of Cplex (AB). We report the overall CPU time in seconds, and compare the quality of lower bounds after reaching the time limit. We also report the number of added columns and generated iterations during the B&P algorithms. From Figure 4.7 we observe that all 120 instances are solved by the Red Cost 2 (resp. ILP) setting at the root node within less than 1000s (resp. 1534s). On the contrary, the Red Cost 1 solves 108 instances with CPU time up to 2700s and the DP setting solves only 71 instances without exceeding the time limit. The fact that the DP consumes more time to nd an elementary resource constrained shortest path, compared to the ILP method, can be explained as follows. In the setting ILP, the sub-tour elimination cuts are added to the pricing problem on the y, only if they are needed. Furthermore, from one CG iteration to the other, we only modify the objective function of the underlying ILP, which allows Cplex to heavily exploit warm-starting techniques.

In order to compare the number of added columns and iterations generated during the CG procedure, we focus on 71 instances for which all four methods were able to solve the problem at the root node without exceeding the time limit of one hour. Red Cost 2 during the column generation procedure. We notice that settings DP and ILP behave in the same way; this is explained by the fact that at each CG iteration, at most one column (the most violated) is added per commodity. On the other hand, for Red Cost 1 (resp. Red Cost 2) several (resp. all) columns having negative reduced costs are added. This explains the huge number of added columns using Red Cost 2. We observe that, as all columns with negative reduced cost are added at each CG iteration, the Red Cost 2 needs less than 5 iterations to nd a solution at the root node. In comparison, the setting DP (resp. ILP) needs up to 20 (resp. [START_REF] Bhamare | Optimal virtual network function placement in multi-cloud service function chaining architecture[END_REF] iterations for all instances to nd the LP-bound. Finally, we observe that Red Cost 1 needs up to 25 CG iterations to solve 65 instances, and sometimes even more than 50 iterations to solve the LP-relaxation.

This can be accounted to the fact that we are adding a xed number of columns, which are not necessarily the most violated ones.

Based on the comparison of CPU times between the four pricing methods, in the following, we will consider Red Cost 2 as our default pricing method for the model PF. This setting which will be denoted by PF in the remainder of this section.

Comparison of LP-relaxation bounds

We now turn our attention to the comparison of the quality of LP-relaxation bounds and the CPU time required for solving LP-relaxation of the formulations PF, DW and C. We also show the relative improvement of lower bounds, with respect to the LPbounds provided by C, by both extended formulations with and without adding valid inequalities. Figure 4.9 provides a cumulative chart in which we report the CPU time for all 120 instances from our test bed. We observe that in less 31 seconds the LP-solution for any of the considered 120 instances can be found by the compact formulation C. The CPU time consumed by the model PF is below 1000 seconds. Finally, the model DW can solve only 91 instances at the root node without exceeding the time limit. In the following, we focus on 52 instances whose LP-relaxation could be solved at the root node by PF, PF+VI, DW, DW+VI without exceeding the time limit (2 among 54 instances solved by DW+VI are not solved by PF+VI). We compare the relative improvement of lower bounds with respect to the LP-bounds obtained by the compact formulation C. Figure 4.12 illustrates the relative improvement of lower bounds, calculated as ((LB a -LB C ) / LB C ) * 100, with a ∈ {PF, PF+VI, DW, DW+VI}. We observe that the lower bound at the root node could be improved by between 1% to 14% (resp. by between 3% to 29%) with PF (resp. DW) formulation for all considered instances without adding valid inequalities. Moreover, the potential benets of adding the valid inequalities to both extended formulations are shown. We notice that adding the valid inequalities further improves the quality of their bounds at the root node.

These improvements are signicant for all instances. The LP-relaxation bounds of the C formulation are improved from 6% to 90% for 75% of instances solved by the model PF. This improvement ranges between 14% and 100% for 75% of instances solved by the model DW. of added columns during the B&P algorithm, the number of branching nodes and the number of generated iterations for PF+VI and DW+VI, respectively. We observe that the setting PF+VI needs more columns and more iterations and also branches more than DW+VI. The vast number of added columns for the setting PF+VI is explained by the fact that we are adding all columns with negative reduced cost at each iteration of the CG procedure. This is in contrast to DW+VI where we are using LP-based pricing method to add only the most violated ones. Furthermore, the number of variables that are required to be integer in the model PF is much larger compared to the model DW, which also explains the larger number of branching nodes for the setting PF+VI. In this chapter we assume that VNF capacity is suciently large, and that the only capacity we have to consider at the nodes is the number of VNFs that we can install.

We call this problem variant the node-capacitated VNFPRP . To improve the capa- In this chapter, we consider two variants for the problem, in the rst variant, we relax only the VNF-capacity constraints, and we call the problem the node-capacitated and conict-constrained VNFPRP. In the second variant, we relax VNF-capacity, nodecapacity and conict constraints, and we call it the Uncapacitated VNFPRP.

Adapted compact MILP formulation

In this section, we rst adapt the compact (i.e., polynomial in size) MILP formulation already introduced in Chapter 2 to model the node-capacitated and conict constrained variant and the uncapacitated variant and the of the VNFPRP. We also show some favorable theoretical properties of this model that allow us to decompose it using Benders decomposition. We conclude this section by proposing three additional families of strengthening inequalities.

Since from constraints (5.11) we have pq∈P r q = 1, it follows that: pq∈P r q l(p q ) > l, which contradicts the fact that the system (5.11)-(5.14) admits a fractional feasible solution.

Proof. [Theorem 5.1] Let X = (x, ŷ, ŵ, t) be an optimal solution of (P ) such that there exists at least one commodity k ∈ C for which the ow tk is fractional. We will show how to construct a purely integer solution X = (x, ŷ, ŵ, t) which is also feasible for (P ). As the components ŷ and ŵ that appear in the objective function with nonzero coecients are the same for both X and X, it follows that the two solutions will have the same objective value.

Let k ∈ C be a commodity such that in the solution X, the associated vector tk contains fractional values. In that case, by setting l := l k , s := s k , d := d k and ĉuv := min

f ∈F k {ŷ f k v + 1 -xfk v + xfk u } ∀(u, v) ∈ A (5.15)
ĉu := max

f ∈F k {ŷ f k u } ∀u ∈ N (5.16)
we obtain the linear inequality system (5.11)- (5.14).

By Lemma 5.2, among the fractional paths composing the ow tk , there exists at least one (integer) path which can be used to substitute tk . This procedure can be repeated for all commodities k ∈ C that admit a fractional ow solution tk in X, without changing the structure of the binary component (x, ŷ, ŵ), due to the fact that the ows are routed separately for each commodity. Corollary 5.3 Without loss of generality, constraints t k uv ∈ {0, 1}, for all (u, v) ∈ A, k ∈ C can be relaxed into t k uv ≥ 0.

Strengthening inequalities

In the following, we derive three families of valid inequalities that can strengthen the quality of the LP-bounds of the formulation (P). The two rst families exploit the precedence constraints imposed for a given commodity, the last one exploit the nodecapacity constraints imposed for each node.

max α s k -α d k - f ∈F k (u,v)∈A β f uv (1 -xf v + xf u + ŷf v ) + u∈N f ∈F k γ f u ŷf u -δl k α u -α v - f ∈F k β f uv + f ∈F k γ f v -δl uv ≤ 0 (u, v) ∈ A (5.25)
α free, (β, γ, δ) ≥ 0

αs k -αd k - f ∈F k (u,v)∈A βf uv (1 -x f v + x f u + y f v ) + u∈N f ∈F k γf u y f u -δl k ≤ 0 (5.27)
These Benders cuts are actively separated at every node of the branch-and-bound tree, and added to the relaxed master problem. The resulting strategy is commonly denoted as Branch-and-Benders-Cut (see, e.g., [START_REF] Cordeau | Benders decomposition for very large scale partial set covering and maximal covering location problems[END_REF]). Eective separation of Benders cuts involves advanced stabilization techniques and proper selection of extreme rays.

To this end, we have decided to use the implementation of Benders feasibility cuts, in an annotated Benders setting of CPLEX [START_REF]Annotating a model for CPLEX[END_REF], which utilizes the major stabilization and implementation techniques from [START_REF] Fischetti | Benders decomposition without separability: A computational study for capacitated facility location problems[END_REF][START_REF] Fischetti | Redesigning Benders decomposition for large-scale facility location[END_REF][START_REF] Fischetti | A note on the selection of Benders' cuts[END_REF].

Remark: This Benders reformulation can be applied on the node-capacitated and the conict constrained version of the problem. This can be done by adding constraints (2.4) and (2.6), respectively, to the Benders master problem. Furthermore, this Benders reformulation can be applied on the compact MILP formulation presented in chapter 2

i.e., VNF-capacitated version of the problem. As the t variables do not appear in the VNF capacity constraints.

MILP-based Heuristic

When models with an exponential number of constraints are implemented within a Branch-and-Cut scheme (as this is the case with our Branch-and-Benders-Cut approach), general purpose solvers in general struggle with nding feasible solutions. This is because they are missing an important information concerning the structure of feasible solutions, as majority of constraints are left out from the model and are separated on the y during the branching process. Hence, only after enumerating many branching nodes, much deeper in the branching tree, a larger number of cuts becomes available to the solver so that more useful feasible solutions can be found.

It is therefore crucial for the performance of Branch-and-Benders-Cut algorithms that

Obtained results

Major results are summarized in this section, whereas more detailed results, given per each instance, can be found in Section 5.4.3. For each of the dierent settings listed above, we report the overall CPU time in seconds and the percentage gap after reaching the time or memory limit (in case the optimal solution has not been found). To nal gap is calculated as GAP = ((U B -LB)/LB) * 100%, where U B denotes the best feasible solution, and LB the global lower bound found in each run.

We start by comparing the four Branch-and-Benders-cut congurations (B, B+VI, B+PH, B+VI+PH) on the set of Erd®s-Rényi graphs. After determining the best conguration, we compare it with three other alternative MILP approaches (namely, C, AB, and RC). In Section 5.4.2.2 we then study the instances derived from SNDlib and compare the best Branch-and-Benders-cut conguration against the available alternatives, and we also study the sensitivity of their performance with respect to the introduction of node-capacity and conict constraints. For example, the number of instances with |N | ∈ {25, 50}, d = 0.5 that can be solved to optimality in less than ten minutes is equal to 20 for B+VI+PH, whereas it is only 14 [START_REF] Alves | A stabilized branch-and-price-and-cut algorithm for the multiple length cutting stock problem[END_REF][START_REF] Allybokus | Virtual function placement for service chaining with partial orders and anti-anity rules[END_REF] for B+VI (B+PH, B, respectively). All 50 instances are solved with a GAP smaller than or equal to 15% by B+VI+PH and B+PH, whereas, this is true for only 41 and 42 instances when settings B, respectively B+VI, are employed. A similar behavior can be observed for instances with |N | ∈ {25, 50} and d = 0.9. These results clearly exhibit the benets of adding valid inequalities and path-based heuristic to the basic Benders decomposition setting. Consequently, in the following we will consider B+VI+PH as our default Benders setting and compare it with alternative MILP approaches. (GAPs). We notice that all models have more diculties in nding optimal solutions Benders reformulation for the node-capacitated VNFPRP for sparser graphs. Nevertheless, the setting B+VI+PH remains the best performing one, followed by the compact formulation C. We observe that 20 instances were solved with Branch-and-Benders-Cut in less than 500 seconds whereas the same number of instances was solved by the compact formulation in less than 2000s. On the contrary, the relaxed compact formulation RC (resp. the Automatic Benders) can solve only 19 (resp. 18) instances within the given time limit of three hours. In terms of the nal gap, Branch-and-Benders-Cut algorithm solves all Benders reformulation for the node-capacitated VNFPRP instances within the nal gap which is below 15%, whereas the same instances were solved with a gap as high as 40% for the compact formulation, 60% for the relaxed compact model, and 75% for the automatic Benders.

Overall, we observe that RC and AB methods provide solutions of poor quality compared to the Branch-and-Benders-Cut and the Compact model. Therefore, for larger instances with 100 nodes, we compare only the B+VI+PH and C settings. Figures 5.6 and 5.7 summarize the results obtained by solving instances with 100 nodes and density d = 0.9 and d = 0.5, respectively. We observe that in the few cases when the number of commodities is rather small, the compact MILP formulation is outperforming (in terms of the CPU time) B+VI+PH. However, by increasing the number of commodities, the compact MILP model struggles with the size of the underlying LP-formulation.

On the contrary, the Branch-and-Benders-Cut algorithm provides relatively small nal gaps. For example, for d = 0.9 nal gaps are below 3% for B+VI+PH, whereas for the compact MILP formulation the nal gaps can be as large as 16% (cf. solved by the Branch-and-Benders-Cut algorithm without exceeding the time limit is 18 (7 are solved to optimality and for 11, the memory limit was exceeded). From Quality of solutions found by the MILP-based heuristic In the following we analyze the quality of solutions found by the MILP-based heuristic (PH) and its overall eciency. Figure 5.8 provides boxplots in which we are comparing the quality of the solutions and the CPU times of the MILP-based heuristic (PH) and the Compact formulation (C). The latter one is considered as an o-the-shelf alternative for nding feasible solutions. In the provided boxplots, all instances derived from Erd®s-Rényi graphs are taken into account. The relative gaps are calculated with respect to the optimal solution or, alternatively, the best-known upper bound found by any of the four exact methods (contrary to Tables 5.2-5.6 from Section 5.4.3 where CPLEX exit gaps are reported). We observe that in most of the cases the heuristic provides feasible solutions of almost the same (or even better) quality compared to those found by the Compact formulation, but within a signicantly less CPU time. By comparing the solution gaps, we notice that the heuristic performs particularly well on sparser graphs. This can be explained by the fact that, by decreasing the graph density, the number of feasible routing paths decreases too, hence the path-based MILP-formulation (which is used within the heuristic) often captures the optimal routing paths. By enlarging the number of nodes or increasing the density, this benet of the MILP-heuristic might be lost, as it can be observed in Figure 5. 8 (a), where e.g., for instances with 100 nodes and d = 0.9 the solutions found by the Compact formulation are of better quality than those found by the MILP-heuristic.

Computational results

185

Eect of adding valid inequalities Finally, we test the potential benets of valid inequalities proposed in Section 5.1.4. Figure 5.9(a) compares the size of the branching tree for the basic setting B and the setting B+VI in which valid inequalities are added to the model. We notice that adding the valid inequalities denitely helps our Benders approach to reduce the number of branching nodes. These reductions can be signicant (a few orders of magnitude). For 50% of the solved instances the number of branching nodes was reduced from ≈ 10 5 to ≈ 10 4 . In Figure 5.9(b) we also show the relative improvement of lower bounds, calculated as (LB B+VI -LB B ) / LB B * 100.

We observe that the lower bound at the root node could be improved by 1% to 3%

for majority of considered instances. constraints. We notice that for 96 out of 100 instances, B+VI+PH nds a feasible solution, and the worst obtained gap is below 47%. On the contrary, less than half of the instances could be solved within the same gap for the remaining three approaches (AB, C and RC).

When it comes to SNDlib instances without node-capacity and conict constraints, For each of the dierent settings listed in Section 5.4, we report the overall CPU time in seconds (or T L/if the time/memory limit is reached, respectively), and the nal gap in percent after reaching the time or memory limit (in case the optimal solution has not been found). To nal gap is reported as Gap = (U B -LB)/LB * 100%, where U B denotes the best feasible solution, and LB the global lower bound found in each run.

In addition, for our MILP-heuristic, in the columns labeled by PH, we report the CPU time of the heuristic (T L means that the time limit of 900 seconds was reached without solving the model to optimality), and the gap with respect to the optimal solution (or, the best known upper bound, in case the latter is not available). Tables 5.7-5.9 (resp. Tables 5.10-5.12) show detailed results of 100 instances derived from SNDlib with (resp. without) node-capacity and conict constraints. The CPU time of all solved instances exceeds the time limit (this explains the removal of t[s] column in tables). Signin Gap column in result tables indicates that no feasible solution is found.

Conclusion

First, it would be interesting to formulate the problem using bi-level optimization.

The leader problem takes care of the VNFs-installation part, and the follower problems routes the trac taking into account these installations. Also, we can tackle the problem using two-stage optimization either to x the VNFs-installation on nodes, or to x the routing paths. This approach would allow a heuristic to obtain a larger instances, high-quality solutions.

Moreover, more ecient separation heuristics and more sophisticated preprocessing methods can be developed in order to improve the resolution of the problem. Furthermore, some meta-heuristics can be used, such as the ant colony optimization algorithm (ACO) or the genetic algorithm. Also, investigating new valid inequalities for the problem would be of interest.

It would be also interesting to further investigate dierent column generation strategies and heuristics to solve the pricing problems, test dierent branching schemes and dene a gap value for the Lagrangian bound to improve the convergence of our Branchand-Price algorithms.

Finally, we can use robust optimization to deal with the uncertainty of some parameters, for example, the number of trac requests in a given period (hour, day, week, ...), in order to minimize the costs generated by the unplanned demands.

Contents 1. 1

 1 Polyhedra and Integer Linear Programming methods . . . 1.1.1 Elements of polyhedral theory . . . . . . . . . . . . . . . . . . 1.1.2 Cutting plane method . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Branch-and-Cut algorithm . . . . . . . . . . . . . . . . . . . . 1.1.4 Benders decomposition . . . . . . . . . . . . . . . . . . . . . . 1.1.5 Dantzig-Wolfe decomposition . . . . . . . . . . . . . . . . . . 1.1.6 Column generation procedure . . . . . . . . . . . . . . . . . . 1.1.7 Branch-and-Price algorithm . . . . . . . . . . . . . . . . . . . 1.2 Graph theory . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Undirected graphs . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Directed graphs . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Shortest path algorithms . . . . . . . . . . . . . . . . . . . . . 1.3 Telecommunication Networks . . . . . . . . . . . . . . . . . 1.3.1 Network Structure . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Network devices . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Network Function (NF) . . . . . . . . . . . . . . . . . . . . . 1.3.4 Routing schemes . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.5 Network Function Virtualization (NFV) . . . . . . . . . . . . 1.3.6 Software Dened Networking (SDN) . . . . . . . . . . . . . . 1.3.7 Service Functions Chaining (SFC) . . . . . . . . . . . . . . . 1.4 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Overview of the related works . . . . . . . . . . . . . . . . . . 1.1 Polyhedra and Integer Linear Programming methods 7

Figure 1 . 1 :

 11 Figure 1.1: Example of a convex hull

Figure 1 . 2 :

 12 Figure 1.2: Example of valid and non valid inequalities and extreme points

Figure 1 . 3 :

 13 Figure 1.3: A hyper-plane separating x * and P

  ȳ∈Y {dȳ + min x≥0 {cx : Ax ≥ b -B ȳ}}, such that the Benders subproblem is dened by φ(ȳ) as follows: (SP ) : φ(ȳ) = min x≥0 {cx : Ax ≥ b -B ȳ}, The Benders subproblem could be feasible, infeasible, or unbounded for a given value of ȳ. In order to prevent infeasibility, we use the LP-duality theory. Let us denote by (D) the dual formulation associated with (SP ) such that: (D) : φ D (ȳ) = max u≥0 {(b -B ȳ) t u : A t u ≤ c} 1.1 Polyhedra and Integer Linear Programming methods 13

  min dy + w s.t. w ≥ (b -By) t u p , p = 1, . . . , P (OptCut) 0 ≥ (b -By) t u r , r = 1, . . . , Q (F easCut) y ∈ Y, w ≥ 0 where OptCut represent the Benders optimality cuts and FeasCut represent the Benders feasibility cuts.

  Dantzig-Wolfe decomposition developed by George Dantzig and Philip Wolfe in 1960. This decomposition aims to reformulate an original problem into a master program and one or multiple pricing problems. Let us consider the following linear program representing a compact formulation (original problem): (C) : min c T x s.t. Ax ≥ b Dx ≥ d x ≥ 0

1. 2

 2 Graph theory 21 cycles. The complexity of Moore-Bellman-Ford's algorithm is in O(|E| × |V |) Algorithm 4: Moore-Bellman-Ford's algorithm(Graph, source): for each vertex v in the graph G \ {s} do dist[v] := +∞ predesessor[v] := ∅ end dist[source] := 0

  algorithm is in O(kn(m + nlogn)) Algorithm 5: Yen's algorithm(Graph, source, destination, K): SP = Get the shortest path using Dijkstra algorithm(Graph, source, destination). Initialize the set of solutions with the shortest path, Sol = [SP ]. Initialize the set of potential paths, P otentialP aths = [ ].

  Some denitions in this section have been collected from the thesis of A. Tomassilli[START_REF] Tomassilli | Towards Next Generation Networks with SDN and NFV[END_REF].

Figure 1 . 4 :

 14 Figure 1.4: Network structure (xed and wireless access networks, core networks and data centers) [84].

FirewallsFigure 1 . 5 :

 15 Figure 1.5: Network devices [3]

Figure 1 . 6 :

 16 Figure 1.6: Routing Schemes [96]

[ 71 ] 29 Figure 1 . 7 :

 712917 Figure 1.7: NFV substitute the network functions by network applications.

Figure 1 . 8 :

 18 Figure 1.8: NFV Architecture [4]

Figure 1 . 9 :

 19 Figure1.9: Traditional network Vs. SDN network[START_REF] Maleki | An SDN Perspective to Mitigate the Energy Consumption of Core NetworksGÉANT2[END_REF]

Figure 1 .

 1 Figure 1.10: SDN architecture

Figure 1. 11 :

 11 Figure 1.11: Example of Service Function Chain [71]

Figure 2 .

 2 Figure 2.1: VNFPRP is composed of multiple subproblems and constraints.

:

  Set of all commodities (trac requests).

F:

  Set of all Virtual Network Functions. (F k , ≺ k ) : Ordered set of VNFs associated with commodity k, k ∈ C. A k : Set of pairs of VNFs which are in conict for commodity k, k ∈ C.

2. 2 . 2

 22 Problem denition Denition 1 (VNFPRP) The Virtual Network Functions Placement and Routing Problem consists of nding for each commodity k ∈ C, the placement of VNFs f ∈ F k at nodes, and the routing paths so that the sum of the VNF installation costs ψ f u plus The Virtual Network Functions Placement and Routing Problem the sum of node activation costs ψ u is minimized. Moreover, the following constraints have to be satised: Node capacity constraints each node u ∈ N has an installation capacity c u , which means that the number of VNFs installed at u should not exceed c u . VNF capacity constraints each function f ∈ F has a capacity m f to manage the amount of data. The sum of bandwidths of all commodities handled by f should be below m f . Routing constraints the s k -d k routing path associated with each commodity k ∈ C should be elementary.

Figure 2 .

 2 Figure 2.2: Example

  2.1) aims to minimize the sum of the VNF installation costs and the node activation costs. Constraints (2.2) are the standard ow conservation constraints which ensure that one unit of ow is routed from s k to d k for each commodity k ∈ C. Constraints (2.3) are the latency constraints: the sum of the arc latency values along the routing path must be less or equal to l k , for each commodity k ∈ C. Constraints (2.4) represent the node capacity constraints; they guarantee that the number of VNFs installed at each node u ∈ N is bounded enough by its capacity c u . Constraints (2.5) are the VNF capacity constraints; they ensure that the volume of data treated by each function f ∈ F does not exceed its capacity m f . Constraints (2.6) are

Figure 2 . 3 :

 23 Figure 2.3: Connected subgraph.

Figure 2 . 4 :

 24 Figure 2.4: Disconnected subgraph.

Figure 2 .

 2 Figure 2.5 shows that the number of variables increases with the number of commodity, which this is caused by the variables indexed by k (x, y and t variables). The number of constraints increases accordingly with the density of the graph: for example, instances P dh having 24 commodities and 68 bidirected arcs, the number of generated variables is small compared to the number of variables. We can also observe that the number of constraints generated by the compact formulation is in general very large compared to the number of variables.

Figure 2 . 5 :

 25 Figure 2.5: Number of variables and constraints generated by the compact formulation.

Theorem 3 .

 3 1 If Yen's algorithm generates all elementary latency-constrained paths for each commodity k ∈ C, then the path-based formulation gives optimal solutions;

3. 1

 1 Path-based MILP Formulation 71 3.1.5 Linear relaxation of path variables λ Constraints λ k p ∈ {0, 1} are the integrality constraints guaranteeing that the latencyconstrained path cannot be split. Together with constraints (3.1), they ensure that there is exactly one path used to route the ow for each commodity. In the following, we show that λ k p ∈ {0, 1} can be relaxed. Let PF' denote the model PF for which the integrality constraints associated with variables λ are replaced by: λ k p ≥ 0, k ∈ C, p ∈ P k . Proposition 3.2 If the relaxed path formulation PF' has an optimal solution with fractional λ values, then it must necessarily admit an integer solution with the same objective value.

Corollary 3 . 3

 33 Without loss of generality, constraints λ k p ∈ {0, 1}, for all k ∈ C, p ∈ P k , can be replaced by λ k p ≥ 0 in the Path-based formulation.

2 . 2 :

 22 The maximum capacity rate associated with each VNF is given as an integer number in[a, b], where a (resp. b) is the minimum (resp. maximum) bandwidth of all demands in the same set of demand (i.e. a = min k∈C b k and b = max k∈C b k ), the capacity is measured on M bits/s. The functions installation cost ψ f u ∈ N is chosen randomly as integer in [50, 1000], and is expressed in dollars $. Virtual Network Functions The total number of demands varies for each SNDlib instance. For each commodity k ∈ C the source node s k , destination node d k and the bandwidth b k in M bits/s are given.

  Service SFC ≤ 60 ms Online Gaming (O-G) NAT-FW-TM-WOC-IDPS ≤ 100 ms Video Streaming (V-S) NAT-FW-TM-VOC-IDPS VoIP NAT-FW-TM-FW-NAT ≤ 500 ms Web Services (W-S) NAT-FW-TM-WOC-IDPS ≤ ∞ ms Other services (O-S) NAT-FW-TM-WOC-VOC Table 3.3: Five services and their respective SFC and latency value. In order to generate feasible instances, the node capacity value is obtained based on the number of demands and the number of VNFs per commodity, which is equal to 5, ∀k ∈ C. Thus, the node capacity is an integer chosen randomly c u ∈ [ |C| × 5 |N | , 2 × |C| × 5 |N | ], ∀u ∈ N . Node activation cost ψ u ∈ N is chosen uniformly at random from [3000, 5000] [65].

Figure 3 . 1 :

 31 Figure 3.1: Comparison between the path-based MILP formulation and the compact MILP formulation with respect to the number of variables and constraints.

Figure 3 .

 3 Figure 3.1 shows the average number of variables and constraints per instance-type.

  gaps, calculated as GAP (%) = ((U B -LB)/LB) * 100, where U B represents the best found feasible solution and the LB represent the global lower bound. The global upper bound is shown in column Costs($), and the value of the LP-relaxation of the two models is given in the last two columns. Detailed results of this study are shown in Tables 3.9-3.13 in Subsection 3.2.4.

Figure 3 . 2 :

 32 Figure 3.2: GAP comparison between PF and C.

Figure 3 . 3 :

 33 Figure 3.3: GAP comparison between PF and C.

Figure 3 . 4 :

 34 Figure 3.4: Costs comparison between PF and C.

Figure 3 .

 3 Figure 3.5 shows the relative improvement of the LP-relaxation value provided byPF with respect to C. We notice that the path formulation provides a much better relaxation bound compared with the compact MILP formulation. In Figure3.5 we observe that for all solved instances the value of the relaxation at the root node provided by the PF formulation is always greater than or equal to the given LP-relaxation of the compact MILP formulation, so closer to the optimal solution.

Figure 3 . 5 :

 35 Figure 3.5: LP-relaxation improvement by the path formulation.

Figures 3 .

 3 Figures 3.6 and 3.7 show the obtained results by varying the number of paths generated by Yen algorithm, with max_paths ∈ {5000, 500, 100, 50} and α = 1.0. Fixing the number of paths for Yen's algorithm means that the demands latency are tightened.

Figure 3 . 6 :

 36 Figure 3.6: Costs comparison between path formulation with 500, 100 and 50 paths.

Figure 3 . 7 :

 37 Figure 3.7: Costs comparison between path formulation with 5000, 500, 100 and 50 paths.

Figure 3 .

 3 Figure 3.8 shows that increasing the function capacities allows for a signicant cost reduction for all instances by -25,96 %. A similar eect can be achieved by increasing node capacities up to -2,78 %.

Figure 3 . 8 :

 38 Figure 3.8:Costs comparison between path formulations with 50 path and with/without increasing node and functions capacities.

3. 2 results 87 Figure 3 . 9 :

 28739 Figure 3.9: Costs comparison between the path formulations with 50 path and with increasing function capacities, bandwidth and latency

4. 1

 1 First extended formulation: the model PF 101 4.1.3 The dual of the master problem Let DPF denote the dual formulation of the LP-relaxation of the model PF. The number of routing paths associated with each commodity can be exponential. Thus, the number of constraints associated to the path variables is exponential in the dual. The associate dual variables are shown next to each constraint from the formulation given in Subsection 4.1.2, in particular, we associate α, η, and π to constraints (4.2), (4.6) and (4.9), respectively.

  )

(4. 25 )

 25 Constraints (4.20) represent the routing constraints ensuring that one path-installation is chosen for each commodity k ∈ C. Inequalities (4.21) represent the nodes capacity constraints. Constraints (4.22) are the VNF-capacity constraints. The LP-relaxation of this model is obtained by replacing (4.23) by τ k p ≥ 0, (4.24) by 0 ≤ w u ≤ 1, and (4.25) by z f u ≥ 0

Proposition 4 . 1

 41 The binary constraints imposed on the arc variables n in the pricing problem (4.27a)-(4.27k) can be relaxed and replaced by n uv ≥ 0, for all (u, v) ∈ A.

4. 2

 2 Second extended formulation: the model DW 109 Therefore, if there exists a feasible solution of the pricing problem with fractional n values, there also exists a latency constrained path (corresponding to binary n values) which satises all the constraints (4.27a)-(4.27k). To show the latter result, one has to follow similar arguments as those given in the proof of Proposition 3.2.4.2.4 Lagrangian boundIn this subsection, we provide details on how we calculate a Lagrangian bound for the DW formulation. Let DDW be the dual of the model DW and let DRDW be the restricted master problem of the model DW. The Lagrangian bound computation is made in the same way as for the model PF, i.e., by constructing a feasible solution for the DDW from the optimal solution of the DRDW. In the following, we show how to construct this feasible solution during the column generation procedure. Let Υ * = (α * , γ * ) be the relevant component of the dual vector obtained by solving DRDW, Z DRDW the associated objective value, and Z * k the optimal value of the pricing problem (see previous section) associated with commodity k with respect to Υ * .

(4. 28 ) 0

 280 To construct a feasible solution Υ D for the DDW, we can use the following transformation: ᾱk = α * k + Z * k , for all k ∈ C. It follows that for all k ∈ C and for all p ∈ T k : (4.28) ⇐⇒ ᾱk -u∈N f ∈F k a f pk u b k γ * f u ≤ Thus, Υ D = (ᾱ * , γ * ) is a feasible solution for DDW and the associated (Lagrangian) bound is calculated as LB = Z DRDW + k∈C Z * k . 110 Extended formulations Given that, any value of LB during the column generation procedure represents a lower bound for the relaxed master problem, in our code, at each iteration of the column generation procedure, the value of the LB is saved and the maximum overall its values is compared to the objective value of the restricted master problem, which we call Primal bound in Figures 4.1 and 4.2. The CG procedure is stopped if the dierence between both values is smaller than ε = 10 -4 , i.e., if |LB -Z DRDW | < ε holds.

Figure 4 .Figure 4 . 1 :Figure 4 . 2 :

 44142 Figure 4.1: Example of the Lagrangian bound evolution during the column generation procedure for the Dantzig-Wolfe formulation on P dh 1 instance without heuristic solution and with valid inequalities.

Proposition 4 . 2 Case 2 .

 422 For any given LP-solution of the (restricted) master problem with fractional τ variables, at least one of the following cases is valid for each fractional commodity k ∈ C: Case 1. There exist two nodes u, v ∈ N \ {s k } satisfying: There exist a function f ∈ F k and a node u ∈ N \ {s k } satisfying:

Figure 4 . 3 :

 43 Figure 4.3: Branching scheme for Dantzig-Wolfe formulation.

(4. 34 )

 34 Proof. If two VNFs f and g are in conict for a commodity k ∈ C, then, both functions cannot be installed at the same node u, if u is activated.Let D k = (F k , E) bethe conict graph associated with commodity k, k ∈ C, where nodes in D k represent the VNFs f ∈ F k . An edge e ∈ E between two nodes f and g in D k represents the fact that f and g are in conict. Let D k denote the set of all maximal cliques in D k and let ω(D k ) be the clique-number (i.e., the size of the maximum clique) in the conict graph. Proposition 4.5 Inequalities (4.35) are valid for the VNFPRP:

Proposition 4 . 8

 48 Inequalities(4.38) are valid for the VNFPRP. u∈N w u ≥ max{1, max k∈C ω(D k )}.

(4. 39 )

 39 Proof. From inequalities (4.10), all VNFs required for each commodity k ∈ C should be installed at graph nodes, thus the number of nodes necessary to install VNFs in conict should be at least equal to the number of VNFs in conict, which is equal to the size of the cliques from D k , k ∈ C.

(4. 45 )

 45 Proof. The number of VNFs installed at node u is at least equal to the number of VNFs needed to handle one commodity k ∈ C.

Figure 4 . 4 :

 44 Figure 4.4: Branch-and-Price algorithm.

4. 4

 4 Branch-and-Price algorithms 123 4.4.3.1 Pricing the columns for the model PF We consider three methods for solving the pricing problem for the formulation PF:

  the support graph obtained from the binary solution of the RPP model. Since negative arc costs are allowed, C may contain subtours. Subtour elimination cuts are added at each iteration of the column generation procedure if C represents one of the two congurations illustrated in Figure 4.5. The connected path shown in Figure 4.5(a) is eliminated by adding the inequality:

Figure 4 . 5 : 2

 452 Figure 4.5: Example of two solutions generated by the relaxed pricing problem.

Proposition 4 .Figure 4 . 6 :

 446 Figure 4.6: Relative improvement of LP-relaxation bounds of the path formulation by using Dantzig-Wolfe formulation on Pdh instances.

  900 seconds. The following settings represent all tested methods in our computational experiments: DP: The model PF (introduced in Section 4.1) with the dynamic programming algorithm proposed in Section 4.4.3.1 as pricing method; Red Cost 1: The model PF with the reduced-cost-based pricing proposed in Section 4.4.3.1, in which a subset of up to 10 paths per commodity with negative reduced cost is added at each iteration of the CG procedure; Red Cost 2: The model PF with the reduced-cost-based pricing proposed in Section 4.4.3.1, in which all paths with negative reduced cost are added at each iteration of the CG procedure; ILP: The model PF with the pricing method based on solving the ILP given in Section 4.4.3.1; PF: The model PF with the best-performing pricing method; PF+VI: The setting PF in which Valid Inequalities (4.33), (4.34), (4.37) -(4.46) and (4.49)-(4.51) presented in Section 4.3 are additionally used to initialize the model; DW: The model DW presented in Section 4.2; DW+VI: The setting DW in which Valid Inequalities (4.38), (4.43), (4.44), (4.49)-(4.51) and (4.52) -(4.54) from Section 4.3 are added to the model; C: The compact MIP formulation proposed in Chapter 2;

4. 6

 6 .1.1 Comparison between dierent pricing methods for solving the LPrelaxation of the model PF Graphical summary of the obtained results comparing the four pricing methods proposed for the model PF is given in Figures 4.7and 4.8. Cumulative charts representing the CPU times (Figure 4.7), the number of added columns (Figure 4.8(a)) and generated iterations (Figure 4.8(b)) during the column generation procedure are provided. A point with coordinates (x, y) in Figure 4.7 indicates that for y instances, the CPU Extended formulations time needed to solve the LP-relaxation was x seconds or less. A similar representation of the number of added columns and generated iterations is shown in Figures 4.8

Figure 4 . 7 :

 47 Figure 4.7: CPU time comparison between dierent pricing methods proposed for the relaxed path formulation.

Figure 4 . 8 :

 48 Figure 4.8: Comparison between four dierent pricing methods proposed for the relaxed path formulation with respect to the number of added columns and generated iterations.

Figure 4 .

 4 Figure 4.8(a) shows the number of added columns for DP, ILP, Red Cost 1 and

Figure 4 .

 4 Figure 4.8(b) allows to compare the number of generated iterations between DP, ILP, Red Cost 1 and Red Cost 2 during the CG procedure. We observe that, as all

Figures 4. 10

 10 Figures 4.10and 4.11 depict the CPU time consumed by PF, PF+VI, DW and DW+VI

Figure 4 . 10 :Figure 4 . 11 :

 410411 Figures 4.10and 4.11 depict the CPU time consumed by PF, PF+VI, DW and DW+VIsettings at the root node of the branching tree. We observe that adding valid inequalities to both formulations increases the CPU time. Valid inequalities slow down the resolution at the root node for the model PF; 75% of instances were solved with CPU time below 300 seconds, whereas after adding valid inequalities, the same percentage of instances needs up to 1025 seconds. Moreover, the overall number of solved instances decreases from 120 to 109. Similarly, for the model DW, after adding valid inequalities the number of instances for which LP-relaxation can be solved drops from 92 to 66 instances.

Figure 4 . 12 :

 412 Figure 4.12: Relative improvement of LP-relaxation bounds of C formulation by using PF and DW formulation with and without valid inequalities.

Figure 4 . 13 :Figure 4 . 14 :Figure 4 . 15 :Figure 4 . 16 :

 413414415416 Figure 4.13: Number of added columns by PF and DW formulation with and without valid inequalities.

Figure 4 . 17 :Figure 4 . 18 :

 417418 Figure 4.17: Gap improvement comparison between the path formulation, Dantzig-Wolfe formulation, the compact formulation and the Automatic Benders of Cplex.

, 4 .

 4 20 and 4.21 compare the overall number

Figure 4 . 19 :Figure 4 . 20 :Figure 4 . 21 :

 419420421 Figure 4.19: Number of added columns comparison between the path formulation and Dantzig-Wolfe formulation.

Table 4 . 6 :

 46 Results for SNDlib instances with Path formulation in which dierent pricing problems are tested.

-eu 4 28 82 378 6 3 . 5 28

 435 30 746.16 79113 4 23.28 1222.18 88498 2 4347.38 1281.69 147 27 -T L --N obel -eu

5

 5 we have proposed two extended formulations for solving the Virtual Network Functions Placement and Routing problem. The variables of the rst formulation (denoted by PF) are latency-constrained paths, whereas the variables of the second formulation (denoted by DW) are latency-constrained paths that also embed the information regarding the function installations at their nodes. In order to strengthen the LP-bounds, we have proposed several families of valid inequalities for both formulations. Their benets have been computationally demonstrated on a set of instances derived from telecommunication networks. We have presented a branching scheme for each formulation and have developed and implemented the associated Branch-and-Price algorithms. The latter ones are computationally compared with the compact MIP formulation presented in Chapter 2 and the automatic Benders decomposition applied to it. The obtained results have shown that there is a trade-o between the quality of global lower bounds and the time needed to solve the LP-relaxations. Whereas the overall best global lower bounds can be obtained by the model DW, its CPU time is sacriced by the expensive MIP-based pricing procedure. Our results show that the full potential of extended formulations is still to be exploited. This can be done by developing problem-tailored pricing procedures, alternative branching schemes, or more advanced heuristics that can be used to initialize the upper bounds, or for solving the pricing problem. Chapter Benders reformulation for the node-capacitated VNFPRP In this chapter, we study the Uncapacitated Virtual Network Functions Placement and Routing problem, for which node-capacity, VNF-capacity, and conict constraints are relaxed. We propose some valid inequalities and a Benders decomposition scheme, allowing us to solve the problem, within a Branch-and-Benders-Cut algorithm. Finally, we provide computational results to compare some variants of the proposed formulations and show signicant improvements over an approach based on automatic Benders cuts (generated with Cplex).

  bilities of nding exact solutions for larger instances of practical relevance, we exploit some theoretical properties of a compact model and develop a Benders decomposition approach, in which the VNF placement problem is treated at the master level and the routing problem (which becomes decomposable per commodity) is solved at the subproblem level. We also propose several new families of valid inequalities and use the path-based MILP formulation to provide heuristic solutions. These elements are combined into an ecient Branch-and-Benders-Cut framework which is capable of beating an o-the-shelf solver (in terms of the CPU time and the overall solution quality) on a set of realistic and random benchmark instances considered in our computational study.Outline of the chapter The chapter is organized as follows. In Section 5.1, we propose a compact MILP formulation based on the model proposed in Chapter 2, study its theoretical properties and introduce new valid inequalities. Section 5.2 is devoted to the Benders reformulation, whereas the Uncapacitated MILP-based heuristic, already dened in Chapter 3, is given in Section 5.3. In Section 5.4, detailed computational results are provided, and some concluding remarks are derived in Section 5.5.

(5. 26 )

 26 Hence, any (extreme) ray (α, β, γ, δ) of the above dual induces a valid Benders feasibility cut:

5. 4 . 2 . 1

 421 Results obtained on Erd®s-Rényi graphs Graphical summary of the obtained results is given in Figures 5.2-5.7, where we provide cumulative charts representing the CPU times and the GAPs for instances with |N | ∈ {25, 50, 100}. Tables 5.1-5.6 in Section 5.4.3 provide detailed results obtained for this set of instances. Comparing four Branch-and-Benders-cut congurations Figures 5.2(a) and 5.3(a) allow to compare the CPU times of the four Branch-and-Benders-Cut settings for the instances with |N | ∈ {25, 50} and d ∈ {0.5, 0.9}. A point with coordinates (x, y) in this chart indicates that for y instances, the CPU time needed to prove optimality was below x seconds. A similar representation of the nal gaps for the same group of instances is shown in Figures 5.2(b) and 5.3(b). Figures 5.2 and 5.3 show that the B+VI+PH setting is outperforming the remaining three Benders settings (B, B+VI, B+PH).

  (a) CPU time (s) (b) GAP (%)

Figure 5 . 2 :

 52 Figure 5.2: CPU time and GAP comparison between Branch-and-Benders-Cut algorithms with and without valid inequalities and with and without MILP-heuristic (|N | ∈ {25, 50} and d = 0.5).

Figure 5 . 3 :

 53 Figure 5.3: CPU time and GAP comparison between Branch-and-Benders-Cut algorithms with and without valid inequalities and with and without MILP-heuristic (|N | ∈ {25, 50} and d = 0.9).

Figure 5 .

 5 5 indicates that for sparser graphs, for 90% (respectively 100%) of the instances the nal gaps obtained by the Branch-and-Benders-Cut remain below 10% (respectively 15%), whereas the corresponding gures for the nal gaps achieved by the compact model are 20% (respectively 40%).

  (a) d = 0.9 (b) d = 0.5

Figure 5 . 4 :

 54 Figure 5.4: CPU time comparison between Branch-and-Benders-Cut algorithm, compact and relaxed compact formulation and Automatic Benders for graphs with 25 and 50 nodes and density 0.9 and 0.5.

  (a) d = 0.9 (b) d = 0.5

Figure 5 . 5 :

 55 Figure 5.5: GAP comparison between Branch-and-Benders-Cut algorithm, compact and relaxed compact formulation and Automatic Benders for graphs with 25 and 50 nodes and density 0.9 and 0.5.

  Figure 5.6(b)).Also, from Figure5.6(a) we notice that 16 out of 22 instances with d = 0.9 are solved within less than 2500 seconds by B+VI+PH, while this is true for only 9 instances when solved using the Compact formulation.

Figure 5 . 6 :

 56 Figure 5.6: CPU time and GAP comparison between Branch-and-Benders-Cut algorithm with valid inequalities and MILP-heuristic and compact formulation (|N | = 100 and d = 0.9).

Figure 5 .Figure 5 . 7 :

 557 Figure 5.7(b) we observe that the largest GAP for B+VI+PH was 17%, whereas the same (or smaller) GAP is achieved for only 23 out of 27 instances using the Compact

Figure 5 . 8 :

 58 Figure 5.8: GAP and CPU time comparison between compact MILP formulation and MILP-based Heuristic.

  (a) The number of branching nodes for settings B and B+VI . (b) Relative improvement of lower bounds at the root node after adding valid inequalities.

Figure 5 . 9 :Figure 5 .

 595 Figure 5.9: Number of branching nodes and bounds improvement of the Branch-and-Benders-Cut approach with and without valid inequalities for graphs with |N | = 100 and d = 0.9.

Figure 5 .

 5 Figure 5.10(a) shows results for SNDlib instances with node-capacity and conict

Figure 5 .

 5 Figure 5.10(b) shows that for all instances B+VI+PH nds a feasible solution and the worst obtained gap is below 27%. On the contrary, only 86%, 55% and 36% of these instances can be solved within the same gap using the Compact formulation C, the Relaxed Compact formulation RC and the Automatic Benders AB, respectively.

Figure 5 . 11 :

 511 Figure 5.11: Average GAPs per instance-type for SNDlib instances with node-capacities and conict constraints.

Figure 5 .

 5 Figure 5.12 shows a similar comparison for the SNDlib instances without nodecapacity and conict constraints. We observe that the Branch-and-Benders-cut algorithm (B+VI+PH) provides small nal gaps compared to those obtained by the Compact formulation for Atlanta, Newyork, Nobel-germany, Nobel-us and Polska. For

Figure 5 . 12 :

 512 Figure 5.12: Average GAPs per instance-type for SNDlib instances without nodecapacities and conict constraints.

  Tables 5.1, 5.2 and 5.3 summarize the results obtained for the instances derived from Erd®s-Rényi graphs with d = 0.5 and with 25, 50 and 100 nodes, respectively. Tables 5.4, 5.5 and 5.6 summarize the results obtained for the instances derived from Erd®s-Rényi graphs with d = 0.9 and with 25, 50 and 100 nodes, respectively. In the last two columns of Tables 5.2-5.6 we report the CPU time per instance required by the MILP-based heuristic (PH) and the relative gap of the obtained solution with respect to the optimal solution (or, alternatively, the best-known upper bound found by any of the four exact methods).

3 : 5 .- 9 - 9 . 5 - 9 50 2192 50 8 - 8 50 2238 50 10 - 9 . 3 - 6 . 7 T

 35995988109367 Results for graphs with |N | = 100 and d = 0.25.3 T L 32.8.1 T L 17.8 T L 54.3 T L 32.6 182 8.1.7 T L 1.1 T L 1.5 T L 1.8 T L 3.1 T L 27.2 T L 6.6 96 1.L 9.3 T L 6.3 T L 5.7 T L 15.9 T L 12.0 348 7.5 Table 5.5: Results for graphs with |N | = 50 and d = 0.9.

Table 5 . 6 :

 56 |C| |F | t[s] Gap t[s] Gap t[s] Gap t[s] Gap t[s] Gap t[s] Gap t[s] Gap t[sResults for graphs with |N | = 100 and d = 0.9.

  SpurN ode = current node i. RootP ath = subpath in Sol[nbP ath -1] from the source node to i. for each path p in Sol do if the same nodes are in path p and the RootP ath then Remove the edge (i, i + 1) from the graph end end for each node j in the RootP ath, with j = SpurN ode do Remove the node j from the graph end SpurP ath = Dijkstra algorithm(Graph, SpurN ode, destination). Create the new path T otalP ath by concatenating the RootP ath and

the SpurP ath: T otalP ath = RootP ath ∪ SpurP ath if T otalP ath / ∈ P otentialP aths then Add T otalP ath to P otentialP aths end Restore edges to the graph. Restore nodes in RootP ath to Graph. end if P otentialP aths is empty then Break end Sort paths in P otentialP aths from the smaller to the bigger cost. Add the lowest path to the set of solutions Sol. Delete the path from P otentialP aths.

Table 2

 2

.1 summarizes our notation.

Table 2 .

 2 1: Main notation, parameters and sets.

Table 2 .

 2

2: Transformation details

  the virtual network function f is installed at or before node u for commodity k; 0, otherwise.

			Binary
	y f k u	1, if virtual network function f is installed at node u for	Binary
		commodity k; 0, otherwise.
	w u	1, if node u is activated; 0, otherwise.

Binary t k uv 1, if arc (u, v) is taken in the path associated with commodity k; 0, otherwise. Binary z f u number of VNF f installed at node u. Integer Table 2.3: Decision variables of the compact MILP formulation

Table 2 .

 2 4 summarizes the more detailed results which are given in Tables 2.6-2.9 (see Subsection 2.7.1). Tables show the obtained results by solving 15 SNDlib instance types with dierent graph topologies and dierent number of commodities, namely, Abilene, Atlanta, Dfn-bwin, Dfn-gwin, Di-yuan, France, Geant, Janos-us, Newyork, Nobel-eu, Nobel-germany, Nobel-us, Pdh, Polska and Ta1. For each instance type, ten instances are generated randomly. From the tables, we observe that only two Abilene instances and two Pdh instances are solved to optimality. Also, we remark that no feasible solution is provided by the model for four Nobel-eu instances and for one instance Ta1.

	2.7 Computational results

Table 2 .

 2

	Instance_name CPU_time(s) GAP(%)	Costs($) Relaxation value
	Abilene	1688.92	4.85	79990.30	58929.67
	Atlanta	T L	58.37	278836.30	92640.77
	Dfn-bwin	T L	1.58	70656.50	48731.37
	Dfn-gwin	T L	8.79	131783.70	100959.17
	Di-yuan	T L	20.49	35104.00	25435.27
	France	T L	67.97	651469.10	187780.30
	Geant	T L	45.80	451038.30	109636.63
	Janos-us	T L	4.10 13775446.70	13204901.84
	Newyork	T L	50.11	342506.20	141291.17
	Nobel-eu	T L	74.45	600379.83	139001.24
	Nobel-germany	T L	45.11	124062.90	58908.65
	Nobel-us	T L	38.93	118642.50	63906.55
	Pdh	1587.34	1.20	54561.60	38874.36
	Polska	T L	31.30	137285.30	88343.00
	Ta1	T L	58.86	376160.33	78467.77

4 represents the average value for the CPU time, nal Gap, cost, and relaxation value of the ten instances associated with each instance type. A value in the column CPU_time(s) represents the average of the CPU time of the instances solved to optimality without exceeding the time limit, whereas T L in this column illustrates the fact that the compact formulation needs more than 1 hour to nd an optimal solution for all instances. From the table, we observe that for some instance types, nding a good feasible solution is easier than for others, for example, for instance types Abilene, Dfn-bwin, Dfn-gwin, Janos-us and Pdh the nal gaps are below 10%, while, for Atlanta, France, Geant, Newyork, Nobel-eu, Nobel-germany, Nobel-us, Polska and Ta1 the nal gaps are between 20% and 74%.

Table 2 .

 2

	58	The Virtual Network Functions Placement and Routing Problem
	2.7.1 Detailed results
	In this subsection, we report detailed results obtained by solving the SNDlib instances
	using the compact MILP formulation.
	Abbreviations Description
	CPU_C	CPU time in seconds.
	RC	LP-relaxation value.
	LB_C	Best known lower bound, provided by Cplex.
	Gap_C	Relative gap, provided by Cplex.
	UB_C	Best known global upper bound, provided by cplex.

4: Average of CPU time, Gap, Cost and relaxation value.

Table 2 . 5

 25 denotes the best feasible solution, and LB the global lower bound found in each run.

	Instance	CPU_C	RC	LB_C Gap_C	UB_C
		T L 96349.60 121732	11.05 108280.35
		T L 70104.38 102227	17.51	84325.84
		T L 50815.42 68999	1.29	68106.14
		T L 48503.40 68227	0.34	67997.68
		1533.09	40452.10	55712	0.00	55712.00
	Abilene	T L 57305.83 74926	2.58	72989.17
		T L 89779.21 111873	7.24 103773.29
		T L 50181.90 67942	4.64	64787.33
		1844.75	40260.73	58900	0.00	58900.00
		T L 45544.12 69365	3.86	66684.41
		T L 187685.36 455555	55.87 201047.36
		T L 65595.14 171546	53.71	79415.00
		T L 75171.64 270340	67.91	86748.61
		T L 65463.21 140690	42.68	80637.43
	Atlanta	T L 79399.97 244878 T L 60813.42 152840	61.30 53.12	94777.48 71651.54
		T L 98220.61 473321	76.86 109548.74
		T L 84582.53 361504	73.80	94728.71
		T L 64238.76 128957	39.33	78239.26
		T L 145237.07 388732	59.10 158977.90
		T L 48543.97 72194	1.71	70960.00
		T L 43190.17 66269	2.36	64703.00
		T L 51086.81 73695	1.26	72765.00
		T L 42190.82 62340	2.76	60618.34
	Dfn-bwin	T L 59611.15 82931 T L 42820.20 65148	2.45 0.60	80896.00 64758.00
		T L 43048.82 63340	2.22	61934.45
		T L 60745.40 84250	0.41	83908.00
		T L 52554.58 73317	0.50	72951.62
		T L 43521.76 63081	1.50	62137.64
		T L 93821.72 129047	11.44 114280.35
		T L 118367.32 151463	11.40 134192.26
		T L 77911.29 113279	11.66 100072.27
		T L 72339.77 105326	11.79	92910.44
	Dfn-gwin	T L 105110.36 132226 T L 101920.72 130431	6.90 123096.30 6.19 122360.65
		T L 105783.60 133469	6.79 124406.07
		T L 123237.29 152082	4.48 145263.03
		T L 86216.46 116757	10.54 104454.14
		T L 124883.21 153757	6.69 143469.36

: Description of abbreviations used on Tables 2.6-2.9.

The nal gap is calculated as GAP_C = (U B -LB)/LB * 100%, where U B

Table 2 .

 2

	60	The Virtual Network Functions Placement and Routing Problem
	Instance	CPU_C	RC	LB_C Gap_C	UB_C
		T L	23367.81	30114	17.55	24829.23
		T L	32295.12	43227	22.46	33518.40
		T L	24334.12	33484	19.14	27076.03
		T L	23803.91	31854	18.16	26069.27
	Di-yuan	T L T L	22808.12 25128.40	34687 37195	29.56 25.33	24431.85 27773.67
		T L	29156.42	40179	16.12	33700.80
		T L	25619.04	39170	27.21	28510.55
		T L	25561.37	32548	13.88	28031.58
		T L	22278.40	28582	15.51	24148.23
		T L	177774.80	673370	69.42	205891.57
		T L	129220.07	646723	77.15	147748.06
		T L	174563.06	704671	72.75	191988.23
		T L	181826.41	637928	68.07	203678.70
	France	T L T L	139979.59 235734.13	690638 743322	76.90 66.02	159531.84 252556.79
		T L	224804.23	693888	64.66	245227.88
		T L	230949.85	609195	59.66	245733.54
		T L	118290.50	463475	68.98	143772.28
		T L	264660.39	651481	56.08	286114.53
		T L	77238.42	147568	26.64	108253.33
		T L	90112.74	167340	25.19	125190.60
		T L	188864.38	1012275	78.06	222085.94
		T L	86437.31	157715	23.96	119922.21
	Geant	T L T L	111260.20 123382.34	191529 1090981	25.49 85.85	142712.27 154365.54
		T L	100105.48	676904	81.02	128477.31
		T L	147328.25	778519	76.54	182635.56
		T L	78866.47	142442	22.11	110949.97
		T L	92770.69	145110	13.17	125993.68
		T L 10887359.23 11382718	4.28 10895750.90
		T L 11342337.18 13070200	13.13 11353640.23
		T L 14005192.67 15379304	8.88 14014082.08
		T L 11873470.30 13822127	14.00 11887417.16
	Janos-us	T L 15534868.92 15570677 T L 14555386.83 14580575	0.16 15545076.64 0.10 14565999.76
		T L 12625202.51 12648938	0.10 12636608.73
		T L 12848111.78 12870288	0.09 12858854.98
		T L 12537209.83 12570924	0.16 12550414.28
		T L 15839879.13 15858716	0.06 15849306.31

6: Obtained results for the compact MILP formulation.

Table 2 .

 2

	2.7 Computational results

7: Obtained results for the compact MILP formulation.

Table 2 .

 2 Decision variables . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.1.2 Linear constraints . . . . . . . . . . . . . . . . . . . . . . . . 68 3.1.3 Mathematical model . . . . . . . . . . . . . . . . . . . . . . . 69 3.1.4 Getting the routing paths . . . . . . . . . . . . . . . . . . . . 70 3.1.5 Linear relaxation of path variables λ . . . . . . . . . . . . . . 71 3.2 Computational results . . . . . . . . . . . . . . . . . . . . . . 72To tackle this, some papers divide the problem into two parts and treat the VNFs placement and the routing problems separately. A large body of the literature deals with heuristics, in order to handle each part of the problem, or for their relaxed versions.

	62	The Virtual Network Functions Placement and Routing Problem
	Instance CPU_C	RC	LB_C Gap_C	UB_C
		2141.91	38805.94	56518	0	56518
		T L 45729.75 66862	0.49	66537.25
		T L 38932.42 51989	1.05	51442.24
		T L 36016.09 49517	1.5	48772.25
	Pdh	T L 30850.73 44116 T L 32338.46 47046	3.33 1.6	42647.65 46291.68
		1032.78	39393.92	55876	0.0	55876
		T L 37821.69 50606	1.82	49684.82
		T L	38734.6	57071	1.64	56135.3
		T L 50119.97 66015	0.54	65660
		T L 82914.01 137978	37.29	86522.85
		T L 106159.63 142278	21.88 111142.83
		T L 85717.67 132934	31.35	91255.38
		T L 85426.18 146856	35.6	94570.84
	Polska	T L 81055.23 130412 T L 94797.85 143218	32.72 28.64 102198.21 87744.08
		T L 95199.63 135530	26.05 100218.98
		T L 91039.84 148651	34.62	97190.71
		T L 92574.06 132398	26.11	97832.97
		T L 68545.87 122598	38.7	75147.37
		T L 82715.84	-	-	-
		T L 73074.28 148201	38.87	90592.97
		T L 63090.38 126016	29.69	88602.56
		T L 65060.97 207975	57.13	89154.76
	Ta1	T L 71865.36 294566 T L 87691.33 306170	66.63 63.5 111764.07 98300.08
		T L 100768.98 925598	85.98 129798.82
		T L 96286.92 960568	87.46	120413.5
		T L 80692.01 250241	56.22 109551.25
		T L 67679.75 166108	44.24	92618.86
	Table 2.9: Obtained results for the compact MILP formulation.

8: Obtained results for the compact MILP formulation.

In this chapter, we propose a path-based MILP formulation to model the VNFPR problem. We also demonstrate how to eciently use it to derive high-quality heuristic solutions in a reasonable computational time. The study is conducted on a set of realistic telecommunication instances derived from the SND library. To test the eciency of our approach, we also compare the obtained results with the compact MILP formulation introduced in Chapter 2. We vary the problem's parameters such as the node and VNF capacities, the commodities bandwidth and latency, and we discuss the trade-os between saved costs and (in)feasibility. The purpose of this chapter is to provide an empirical study and evaluate the viability of the MIP-based heuristic, as an alternative to the compact MILP formulation presented on Chapter 2. 3.2.1 Benchmark instances . . . . . . . . . . . . . . . . . . . . . . . 72 3.2.2 Models analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.2.3 Obtained results . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.2.4 Detailed results . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 In this chapter, we propose a computationally eective path-based MILP formulation to model the Virtual Network Functions Placement and Routing Problem. The Pathbased Formulation (PF) is established using Yen's algorithm, which aims to nd a xed number of elementary shortest paths between two nodes. Based on that, the latency-constrained routing paths associated with each commodity are obtained.

In the case where all feasible paths associated with each commodity are taken into account in the model, the path formulation represents an exact method; thus, it provides optimal solutions. Alternatively, if only a subset of feasible paths is included in the model, the path-based formulation provides heuristic solutions. We demonstrate the path-based formulation's eectiveness by comparing it with the MILP compact formulation presented in Chapter 2 on a set of realistic benchmark instances derived from the SNDLib. In addition, we vary the problem's parameters, such as node and VNF capacities and commodities latency and bandwidth, and analyze the computational behavior and the cost saving achieved by enlarging the capacities.

  .1.We keep the variables x, y, z and w already dened for the compact MILP formulation in Chapter 2 and we dene a new family of variables representing the routing variables.For each commodity k ∈ C, let P k denote the set of all shortest s k -d k -paths whose total length does not exceed l k , where the arc latency is used as the length measure. For each commodity k, we suppose that all feasible paths in P k are given. Let t pk uv be the parameter that is equal to 1 if arc (u, v) belongs to path p for commodity k; and equal to 0 otherwise.

	Variables	Type

λ k p 1, if path p ∈ P k associated with commodity k is chosen; 0, otherwise. Binary x f k u 1, if the virtual network function f is installed at or before node u for commodity k; 0, otherwise. Binary y f k u 1, if virtual network function f is installed at node u for commodity k; 0, otherwise. Binary w u 1, if node u is activated; 0, otherwise. Binary z f u number of VNF f installed at node u. Integer Table 3.1: Decision variables of the path-based formulation

3.1.2 Linear constraints

The constraints containing x, y, z and w variables are the same for the path-based MILP formulation as in the compact formulation. As we replace the arc variables t in the compact formulation, by the routing variables λ in the path formulation, the ow (2.2) constraints, the precedence constraints (2.7), and the linking constraints (2.10), will be replaced by the following constraints:

Table 3 .

 3

	3.2 Computational results

4: Details about the instances. |N |: the number of nodes, |A|: the number of arcs, |C|: the number of demands, |F |: the number of functions, # AAC: the number of anti anity constraints.

Table 3 .

 3 .5 summarizes the average demand per service for each instance type.

	Instance_type	O-G	V-S VOIP	W-S	O-S
	Abilene	4.2	6.2	7.4	57.6	55.6
	Atlanta	1.7	6.1	4.9	98.7	97.6
	Dfn-bwin	0.0 22.1	21.4	21.9	23.6
	Dfn-gwin	0.0 23.6	23.5	31.0	30.9
	Di-yuan	0.6	1.6	1.2	8.1	9.5
	France	2.5 11.0	8.8 140.1 136.6
	Geant	0.8 20.9	17.5 208.1 213.7
	Janos-us	10.5 24.2	22.5 302.9 288.9
	Newyork	4.5 12.8	11.3 105.6 104.8
	Nobel-eu	0.3 10.8	9.4 168.9 187.6
	Nobel-Germany	0.0	6.6	5.7	52.0	55.7
	Nobel-us	1.6	5.3	6.2	39.0	37.9
	Pdh	0.0	6.3	4.9	5.5	6.3
	Polska	0.0	5.0	4.1	27.5	28.4
	Ta1	6.5 20.5	20.1 169.2 178.7

5: Average demands per service for each instance type.

  Both formulations have practically the same set of variables, except for the routing part. For the path formulation there are |C| × |P k | path variables, for the compact MILP formulation there are |C| × |A| arc variables. The set P k may contain an exponential number of elements. Constraints: The number of constraints in the path formulation is smaller than the one in the compact MILP formulation. Accordingly, the compact MILP formulation has |C| × |N | more constraints.

	# variables	PF C	# constraints	10 5 10 6	PF C
				10 4
	Abilene Atlanta Dfn-bwin Dfn-gwin Di-yuan France Geant Janos-us Newyork Nobel-eu Nobel-ge Nobel-us Pdh Polska Ta1 Instances		Abilene Atlanta Dfn-bwin Dfn-gwin Di-yuan France Geant Janos-us Newyork Nobel-eu Nobel-ge Nobel-us Pdh Polska Ta1 Instances
	(a) # Variables			(b) # Constraints

Table 3 .

 3

	3.2 Computational results

6: Nature of the path formulation (exact or heuristic) and a comparison between both formulations with respect to the number of instances solved to optimality (#optimal) and the number of instances for which no feasible solution is found.

Table 3 .

 3 Table 3.7 no nal Gaps are shown for PF for instance types Dfn-bwin and Dfngwin because PF is a heuristic for them, we are not generating all feasible paths for these instances.

7: Average CPU time, Gap, Cost and relaxation value comparison between PF and C.

In

Table 3 .

 3

	3.2.4 Detailed results
	Abbreviations Description
	CPU_H	Heuristic CPU time in seconds.
	RH	Heuristic LP-relaxation value.
	LB_H	Heuristic best known lower bound, provided by Cplex.
	Gap_H	Heuristic relative gap, provided by Cplex.
	UB_C	Heuristic best known global upper bound, provided by Cplex.
	CPU_C	Compact formulation CPU time in seconds.
	RC	LP-relaxation value of the compact formulation.
	LB_C	Best known lower bound, provided by Cplex for the compact for-
		mulation.
	Gap_C	Relative gap, provided by Cplex for the compact formulation.
	UB_C	Best known global upper bound, provided by Cplex for the compact
		formulation.

8: Description of Tables 3.9-3.13 abbreviations Notation `T L in tables means that the time limit is exceeded and the signmeans that no feasible solution is found by the model.

Table 3 .

 3 9: Results comparison between MILP-Based Heuristic and compact MILP formulation.

	90

Table 3 .

 3 10: Results comparison between MILP-Based Heuristic and compact MILP formulation.

	3.2 Computational results

Table 3 .

 3 11: Results comparison between MILP-Based Heuristic and compact MILP formulation.

	92

Table 3 .

 3 12: Results comparison between MILP-Based Heuristic and compact MILP formulation.

	3.2 Computational results

Table 3 .

 3 13: Results comparison between MILP-Based Heuristic and compact MILP formulation.

	94

  and express the fact that if VNF f is installed at node u for the commodity k, then f is installed at or before the node u. Constraints (4.9) ensure that if a VNF f ∈ F k is installed at a node u for a given commodity k, then the associated routing path p must enter that node. Constraints (4.10) guarantee that all required functions for commodity k ∈ C are installed at the graph nodes. Finally, constraints (4.11) guarantee that, for each commodity k ∈ C, no VNF is installed at or before the source node s k and all VNFs are installed at or before the destination node d k . Model (4.1)-(4.13

		11)
	(λ, x, y, w) are binary	(4.12)
	z is integer	(4.13)

Constraints (4.2) are the path constraints which ensure that exactly one elementary latency-constrained path p ∈ P k is chosen for each commodity k ∈ C. Constraints (4.3) represent the node capacity constraints, which guarantee that the number of VNFs installed at each node u ∈ N is bounded by its capacity c u . Constraints

(4.4)

are the VNF capacity constraints. They ensure that the volume of data treated by each function f ∈ F should not exceed its capacity m f . Constraints (4.5) are the conict constraints and they guarantee that two VNFs in conict are not installed at the same node u ∈ N . Constraints (4.6) are needed to link node installation variables (y), precedence variables (x) and path variables (λ): the left-hand-side is forced to 1 (implying that the function f is installed at the node v) if and only if (i) the path p passing through the arc (u, v) is chosen for the considered commodity k and (ii) the function f is installed at or before the node v and it is not installed at or before the node u. Constraints (4.7) impose the VNFs order for each commodity. Inequalities (4.8) link the precedence and the installation variables, x and y,

  .2.

	Variables	Type

τ k p 1, if path-installation p associated with commodity k Binary is chosen; 0, otherwise. w u 1, if node u is activated; 0, otherwise. Binary z f u number of VNF f installed at node u. Integer Table 4.2: Decision variables of the Dantzig-Wolfe formulation

  .[START_REF] Bonami | Implementing Automatic Benders Decomposition in a Modern MIP Solver[END_REF] where we associate variables α and γ to constraints (4.20) and(4.22), respectively. The separation of constraints (4.26) represents the pricing problem. Let (α * , γ * ) be components of the dual solution of the restricted master problem, the pricing problem consists of nding a commodity k and a path p ∈ T k such that:

Table 4 .

 4 virtual network function f is installed at or before node u; 3: Decision variables of the pricing problem for the model DW The MILP formulation of the pricing problem is given as follows:

	Binary

h f u 1, if virtual network function f is installed at node u; 0, otherwise. Binary n uv 1, if arc (u, v) belongs to the routing path; 0, otherwise. Binary

  constraints between variables d and h, the right-hand-site is forced to 1 in order to ensure that if VNF f is used at node u, then it is used at or before u. Inequalities (4.27g) link variables h and n, they guarantee that the routing path enter all nodes at which VNFs are installed. Finally, constraints (4.27h) guarantee that all required VNFs for the current commodity are installed at nodes.

	Constraints (4.27d),
	(4.27e), (4.27i) and (4.27j) represent the precedence constraints. Constraints (4.27f)
	are the linking

h, n) is binary (4.27k) Constraints (4.27a) are the ow-preservation constraints ensuring that the path goes from the source node s k to the destination node d k . Inequalities (4.27b) represent the latency constraints. Inequalities (4.27c) are the anti-anity constraints which ensure that two VNFs f and g in conict are not installed at the same node u.

  1, there must exist another fractional path-installation p 2 ∈ T k \ {p 1 } passing through node u on which VNF f is installed. Since p 1 = p 2 , we have two cases: (I) p 1 and p 2 pass through the same nodes but with dierent function installations, i.e., there exists at least one VNF g ∈ F k \ {f } installed on a dierent node. Let denote by v (resp. w, v = w) the node belonging to the pathinstallation p 1 (resp. p 2 ) on which g is installed. As the hypothesis (4.32) is valid for any VNF in F k and any node in N \ {s k },

  .[START_REF] Chitimalla | 5g fronthaul-latency and jitter studies of cpri over ethernet[END_REF] Proof. For a given commodity k ∈ C, if the VNF f ∈ F k is installed at node u, then node u should be activated.Proposition 4.4 Inequalities (4.34) are valid for the VNFPRP:

  .[START_REF] Mosharaf | A survey of network virtualization[END_REF] Proof. For a given commodity k ∈ C, nodes in Q represent the set of VNFs in conict, i.e., they cannot be installed at the same node. Therefore, only one VNF in Q can be installed at node u as otherwise the conict constraints (4.5) are violated.

	Linear inequalities (4.34) and (4.35) can be combined and generalized for each clique
	Q in D k .
	Proposition 4.6 Inequalities (4.36) are valid for the VNFPRP, and they also domi-
	nate inequalities (4.35).

f ∈Q

  If a VNF f associated with commodity k ∈ C is installed at the destination node, then, this function cannot be installed at or before any predecessor of d k . Inequalities (4.46) can be generalized by the following inequalities (4.47). Let us consider a cut separating the source nodes s k from the destination node d k in G. We denote by N s k (resp. N d k ) the component containing s k (resp. d k ) in G. Suppose that arcs in the cut go from N s k to N d k . The following inequalities are valid for the VNFPRP.

	4.3 Strengthening inequalities	117
	Proof.

.46)

where Γ -(d k ) denotes all incoming neighbors of d k .

  For a given commodity k ∈ C, if a VNF f is installed at nodes in N d k and in addition we have no arc going from N d k to N s k , then f is not installed at or before nodes in N s k .

Proposition 4.13 Inequalities (4.46) can be generalized by inequalities (4.48) which are valid for the VNFPRP. Let S k be the node separator disconnecting N s k from N d k ,

  For a given commodity k ∈ C, if the VNF f is installed at nodes in N d k and it is not installed at or before nodes in the separator S k , then, the chosen path does not go from N d k to N s k ; This means that the function f cannot be installed at or before nodes of N s k . Note that arcs from N d k to N s k may exist. 4.14 Node capacity constraints(4.3) do not dene facets of the polyhedron of the VNFPRP if there exists a node u, such that

	4.3.2 Strengthening inequalities for both models
	Besides inequalities (4.38), (4.43) and (4.44) which are also valid for the model DW, in
	the following we propose additional inequalities that involve only z and w variables, and
	are therefore valid for both formulations studied in this chapter. With valid inequalities
	given in Proposition 4.11 we address the setting in which the capacity of a function is

not sucient to handle the full demand of a given commodity (i.e., multiple copies of the same function need to be installed). Proposition 4.14 provides further generalizations of this setting. Proposition

  .49) 2) Moreover, if |C u | < |C|, then the linear inequalities (4.50) dominate (4.49):

  All proposed methods reach the time limit for the 52 instances for which all methods could solve the LP-relaxation. Thus, in what follows we show only the Gap improvement for all methods, which is calculated as GAP = (GLB -LB)/LB * 100%, where GLB represents the global lower bound found by each setting listed above, (the best known solution provided by Cplex for C and AB settings, and the best global lower bound provided by the Branch-and-Price algorithm based BFS strategy for PF+VI and DW+VI settings). Also, we compare the number of added columns and generated nodes and iterations by both Branch-and-Price algorithms. Figures 4.17and 4.18 display the relative Gap improvement of lower bounds. We notice that the relative Gap is improved by 5% for 75% of solved instances (i.e., 39 instances) by DW+VI, by 0.6% with PF+VI, by 2.5% with the Compact formulation and by 2.1% with the Automatic Benders of Cplex. Moreover, for 50 out of 52 instances, the Gap is improved by 10% with Dantzig-Wolfe formulation, whereas for 50 instances the path formulation improves the gap only by 2.5%, the Compact formulation and the automatic Benders by 5%. Therefore, Dantzig-Wolfe formulation provides a better solution quality.

	142	Extended formulations
	In our Branch-and-Price (B&P) algorithms we initialize UBs by using one of the two
	heuristics dened above. Consequently, we implemented BFS (Breadth-First Search)
	strategy-based branching in which a global lower bound is updated at each level of the
	B&P tree.

Table 4 .

 4 5: Results for SNDlib instances with Path formulation in which dierent pricing problems are tested.

	4.6.2 Detailed results
	Tables provided in this subsection show more detailed results for each setting intro-
	duced in Section 4.6. First, we describe abbreviations shown in each column in the
	tables:

Table 4 .

 4 7: Results for SNDlib instances with Path formulation in which dierent pricing problems are tested.

	N obel -us 10 14	N obel -us 9 14	N obel -us 8 14	N obel -us 7 14	N obel -us 6 14	N obel -us 5 14	N obel -us 4 14	N obel -us 3 14	N obel -us 2 14	N obel -us 1 14	N obel -ger 10 17	N obel -ger 9 17	N obel -ger 8 17	N obel -ger 7 17	N obel -ger 6 17	N obel -ger 5 17	N obel -ger 4 17	N obel -ger 3 17	N obel -ger 2 17	N obel -ger 1 17
	42	42	42	42	42	42	42	42	42	42	52	52	52	52	52	52	52	52	52	52
	91	91	91	91	91	91	91	91	91	91	121	121	121	121	121	121	121	121	121	121
	6	6	6	6	6	6	6	6	6	6	6	6	6	6	6	6	6	6	6	6
	89.50	63.50	86.47	95.96	102.52	90.21	96.13	50.70	106.22	88.14	204.65	129.05	183.60	154.46	192.33	236.01	142.38	181.51	248.75	243.14	T L
	946	895	958	941	982	942	846	924	932	921	1321	1346	1289	1328	1400	1340	1263	1281	1329	1320	-
	7 46.87	6 37.02	7 42.39	8 42.98	8 51.08	7 57.38	8 28.27	5 24.75	7 57.42	7 29.19	9 102.34	8 177.77	8 74.09	8 98.02	9 159.23	9 121.14	7 71.24	7 70.12	9 104.88	11 97.41	-747.60
	941	894	962	949	1005	948	837	929	930	918	1308	1358	1299	1338	1408	1374	1260	1284	1324	1302	6423
	7	6	8	8	8	9	7	5	7	6	10	8	8	8	10	11	7	8	8	10	14
	1.88 2220	2.08 2040	2.72 2320	2.19 2257	4.66 2293	2.27 2123	2.25 1841	1.81 2175	3.14 2247	1.70 1958	15.45 4212	21.26 4976	12.09 3897	23.55 4634	20.40 4889	19.95 4243	15.54 3734	18.03 4369	18.80 4129	14.79 4098	1941.05 34771
	4	4	4	4	5	5	4	3	4	3	8	14	7	14	13	11	9	13	11	8	34
	1.66 2324	1.96 2160	2.59 2483	2.06 2402	5.01 2482	1.75 2271	2.27 1967	1.87 2322	2.73 2385	1.76 2013	10.99 6051	8.42 6026	10.19 5300	11.51 6032	13.77 5809	10.36 5724	8.58 5028	8.82 5783	12.50 5705	10.18 5900	681.54 92449
	3	3	3	3	3	3	4	3	3	3	3	3	3	3	4	3	3	3	3	3	3

Table 4 .

 4 8: Results for SNDlib instances with Path formulation in which dierent pricing problems are tested.

	3
	269.40 44598
	19
	22210
	460.19
	14
	-1534.39 5260
	-
	T L

Table 4 .

 4 9: Results for SNDlib instances with relaxed Path formulation, relaxed Dantzig-

	Wolfe formulation and Relaxed Compact formulation.

Table 4 .

 4 10: Results for SNDlib instances with relaxed Path formulation, relaxed Dantzig-Wolfe formulation and Relaxed Compact formulation.

	154

Table 4 .

 4 11: Results for SNDlib instances with relaxed Path formulation, relaxed Dantzig-Wolfe formulation and Relaxed Compact formulation.

	Instances PF PF+VI DW DW+VI Gap t[s] |N | |A| |C| |F | Col Iter Gap t[s] Col Iter Gap t[s] Col Iter Gap t[s] Col Iter 12 30 132 6 3.32 0.78 1718 3 47.56 3.19 1794 1 7.83 133.49 330 8 48.50 311.32 2097 35 12 30 132 6 5.77 0.57 1718 2 247.07 2.29 1664 1 17.03 55.75 282 6 248.86 153.00 1982 29 12 30 132 6 4.44 0.48 1687 2 89.31 1.86 1764 1 12.43 51.87 218 5 95.19 248.14 1748 35 12 30 132 6 4.60 0.58 1698 2 168.61 1.92 1774 1 8.58 61.98 334 7 172.38 171.88 1809 30 12 30 132 6 13.78 0.49 1683 3 478.48 0.69 1748 1 17.42 72.64 214 5 478.48 36.06 917 10 12 30 132 6 5.12 0.48 1689 2 570.43 0.81 1752 1 11.30 52.01 232 6 570.43 236.55 2576 40 12 30 132 6 5.50 0.42 1680 2 160.11 1.04 1740 1 8.92 57.97 220 7 161.50 201.55 2080 32 12 30 132 6 1.04 0.53 1695 2 51.77 1.65 1770 1 5.05 61.45 437 7 55.91 225.81 2146 34 12 30 132 6 8.70 0.52 1679 2 220.14 2.15 1762 1 19.77 57.34 308 5 222.68 266.56 3031 41 12 30 132 6 3.61 0.65 1717 3 438.28 0.97 1782 1 10.52 43.41 240 6 438.57 188.90 2083 33 15 44 210 6 5.48 11.74 7392 2 15.93 45.32 7401 2 13.69 503.66 501 8 24.71 1845.42 2296 29 15 44 210 6 4.49 10.49 4563 4 27.17 24.47 4696 2 11.39 172.56 505 5 31.90 2303.34 3239 37 15 44 210 6 7.81 6.90 4929 3 112.10 30.09 5018 2 16.68 302.41 206 6 121.42 2422.32 5144 41 15 44 210 6 7.01 6.54 4996 3 79.83 15.64 5307 2 14.13 281.46 207 6 84.59 1779.73 2929 28 15 44 210 6 6.85 5.75 4151 3 99.32 4.08 6540 1 11.09 330.98 417 8 101.44 2782.86 5653 54 15 44 210 6 13.05 10.02 4930 3 93.83 21.75 4825 2 19.76 194.32 358 7 97.71 2833.35 4342 45 15 44 210 6 14.57 7.42 4985 3 113.69 30.85 4911 2 29.05 197.63 208 6 122.65 1383.84 3133 29 15 44 210 6 5.62 8.46 5046 3 49.84 26.84 5115 2 14.34 502.80 341 8 58.50 1602.83 2688 26 15 44 210 6 5.22 5.85 4851 4 70.18 24.97 4869 2 13.54 316.95 446 10 78.37 1409.37 3245 31 15 44 210 6 9.85 7.16 5046 3 105.56 31.01 4958 2 15.48 358.66 370 12 108.92 2025.36 3789 37 11 84 22 6 1.37 13.51 7361 2 5.96 7.42 5386 2 4.98 228.98 117 10 9.37 1037.92 550 39 11 84 22 6 3.43 13.42 8227 2 10.87 9.40 4817 2 12.85 281.78 136 16 20.01 2025.82 551 38 11 84 22 6 0.39 27.97 14277 2 3.96 20.58 10295 2 3.81 384.60 135 11 -T L --11 84 22 6 1.58 22.35 13831 2 18.71 13.99 7343 2 9.36 408.55 99 10 -T L --11 84 22 6 0.97 22.14 13549 2 18.53 14.65 7724 2 5.82 216.27 103 9 -T L --11 84 22 6 1.11 16.21 10228 3 7.58 13.92 7263 2 7.20 355.63 143 15 12.16 2389.28 635 48 11 84 22 6 11 84 22 6 1.74 16.37 12235 2 21.03 14.00 8060 2 4.75 167.02 86 9 11 84 22 6 0.39 34.84 17608 2 6.57 28.17 15888 2 9.88 567.53 110 10 11 84 22 6 1.94 20.52 13190 2 10.41 18.02 9622 2 6.25 151.02 101 8 25 90 300 6 10.61 373.18 55964 3 44.53 990.71 56899 2 -T L --25 90 300 6 5.10 272.30 47514 5 17.22 1640.24 48485 2 -T L --25 90 300 6 6.24 333.01 54355 3 48.88 2769.04 57207 1 14.39 3586.68 500 6 25 90 300 6 10.90 241.07 50397 3 29.38 500.92 52651 2 -T L --25 90 300 6 9.11 251.89 47509 3 32.25 1285.48 50936 2 -T L --25 90 300 6 6.61 326.58 55720 4 45.30 2558.16 57172 2 14.57 3157.17 451 7 25 90 300 6 10.77 301.56 55379 3 21.58 958.75 57837 2 -T L --25 90 300 6 5.86 345.98 46745 4 20.79 420.02 50014 2 -T L --25 90 300 6 5.02 299.81 50328 3 14.56 405.80 49570 2 -T L --25 90 300 6 13.10 296.71 50403 3 74.57 503.32 51435 2 22.46 2700.80 409 5 0.91 19.95 13053 2 14.00 13.41 7036 2 3.86 237.66 114 10 17.26 2652.03 613 46 -T L ---T L ---T L ---T L ---T L ---T L ---T L ---T L ---T L --4402.30 2165.80 1121 98 -T L ---T L -4403.37 3427.42 2103 239 -
	Name	Abilene 1	Abilene 2	Abilene 3	Abilene 4	Abilene 5	Abilene 6	Abilene 7	Abilene 8	Abilene 9	Abilene 10	Atlanta 1	Atlanta 2	Atlanta 3	Atlanta 4	Atlanta 5	Atlanta 6	Atlanta 7	Atlanta 8	Atlanta 9	Atlanta 10	Di -yuan 1	Di -yuan 2	Di -yuan 3	Di -yuan 4	Di -yuan 5	Di -yuan 6	Di -yuan 7	Di -yuan 8	Di -yuan 9	10 Di -yuan	F rance 1	F rance 2	F rance 3	F rance 4	F rance 5	F rance 6	F rance 7	F rance 8	F rance 9	F rance 10

Table 4 .

 4 12: Results for SNDlib instances with relaxed Path formulation, relaxed Dantzig-Wolfe formulation with and without valid inequalities.

	Instances PF PF+VI DW DW+VI Gap t[s] |N | |A| |C| |F | Col Iter Gap t[s] Col Iter Gap t[s] Col Iter Gap t[s] Col Iter 22 72 462 6 2.83 382.87 50877 3 484.94 400.56 41009 2 14.47 1510.34 414 4 -T L --22 72 462 6 6.65 298.15 48652 3 371.93 707.50 45857 2 -T L ---T L --22 72 462 6 7.30 290.98 47741 3 319.00 1077.01 44656 2 16.21 1310.94 430 3 -T L --22 72 462 6 2.59 269.43 45164 4 53.49 723.15 45379 2 6.00 1848.17 1306 7 -T L --22 72 462 6 4.42 303.14 49197 3 431.34 381.11 34961 2 17.45 2429.54 864 5 -T L --22 72 462 6 3.92 271.81 45707 4 157.94 728.84 34729 2 11.27 2346.37 784 5 -T L --22 72 462 6 3.77 283.34 49418 3 130.91 2086.82 49933 2 -T L ---T L --22 72 462 6 7.77 293.47 51436 3 182.56 1647.43 50908 2 14.24 2779.17 436 5 -T L --22 72 462 6 2.58 270.50 46695 4 110.81 2274.09 43065 2 5.29 1838.10 704 6 -T L --22 72 462 6 5.15 286.78 50382 3 457.67 383.16 46688 2 18.57 2312.52 423 4 -T L --16 98 240 6 3.55 823.65 114459 3 20.23 873.57 117283 2 -T L ---T L --16 98 240 6 4.54 616.65 107859 3 36.77 1371.14 113142 2 -T L ---T L --16 98 240 6 2.11 649.80 108805 3 17.36 1407.36 113158 2
	Name	Geant 1	Geant 2	Geant 3	Geant 4	Geant 5	Geant 6	Geant 7	Geant 8	Geant 9	Geant 10	N ewyork 1	N ewyork 2	N ewyork 3

Table 4 .

 4 ger 13: Results for SNDlib instances with relaxed Path formulation, relaxed Dantzig-Wolfe formulation with and without valid inequalities.

	17 52 121 6 3.26 10.18 5900 3 79.45 5.06 10180 2 10.21 443.90 434 8 -T L --17 52 121 6 7.71 12.50 5705 3 129.04 9.63 7903 2 19.61 749.74 459 11 -T L --17 52 121 6 4.44 8.82 5783 3 160.83 9.38 9050 1 10.95 312.16 359 8 -T L --17 52 121 6 4.05 8.58 5028 3 45.33 7.95 8582 2 16.80 642.54 463 8 -T L --17 52 121 6 4.15 10.36 5724 3 57.83 7.49 9388 2 11.12 492.78 572 9 -T L --17 52 121 6 2.82 13.77 5809 4 19.48 39.88 5716 2 7.18 605.83 470 10 -T L --17 52 121 6 6.38 11.51 6032 3 32.35 14.71 8215 2 12.97 416.67 415 7 -T L --17 52 121 6 9.33 10.19 5300 3 105.01 31.42 6854 2 21.36 662.80 318 7 -T L --17 52 121 6 4.45 8.42 6026 3 56.46 4.98 9415 2 16.28 727.86 624 11 -T L --17 52 121 6 4.89 10.99 6051 3 45.50 30.36 6504 2 13.23 703.07 457 11 -T L --
	1	N obel -ger 2	N obel -ger 3	N obel -ger 4	N obel -ger 5	N obel -ger 6	N obel -ger 7	N obel -ger 8	N obel -ger 9	N obel -ger 10

Table 4 .

 4 14: Results for SNDlib instances with relaxed Path formulation, relaxed Dantzig-Wolfe formulation with and without valid inequalities.

Table 5 .

 5 Benders reformulation for the node-capacitated VNFPRP Gap 100 4996 10 6 101.22 0.00 329.27 0.00 88.37 0.00 100 5058 10 8 64.18 0.00 65.08 0.00 35.65 0.00 100 5026 10 10 78.33 0.00 106.68 0.00 43.32 0.00 100 4974 15 6 377.13 0.00 7081.53 0.00 61.22 0.00 100 5000 15 8 622.40 0.00 855.60 0.00 32.49 0.40 100 4966 15 10 T L 2.42 T L 1.14 31.95 0.35 100 5030 20 6 522.76 0.00 T L

	100 4974 50 10 -10.94 T L 15.49 67.43 0.48	100 4954 50 8 -14.33 T L 19.38 106.38 0.00	100 5026 50 6 -16.91 T L 29.20 527.16 0.00	100 5050 45 10 -6.64 T L 11.25 57.13 1.19	100 4954 45 8 T L 0.58 T L 0.06 44.86 0.35	100 4910 45 6 -13.93 T L 29.04 156.51 0.04	100 5060 40 10 T L 5.97 T L 14.85 74.16 0.00	100 4946 40 8 -2.34 T L 10.86 75.34 0.00	100 4986 40 6 -7.52 T L 8.54 125.36 0.35	100 4820 35 10 T L 5.45 T L 9.20 64.12 0.00	100 4878 35 8 -7.83 T L 9.99 69.30 0.58	100 4966 35 6 T L 6.12 T L 16.30 49.54 0.25	100 4866 30 10 -11.35 T L 25.15 73.07 0.67	100 4984 30 8 T L 8.71 T L 12.47 73.05 0.00	100 4896 30 6 T L 4.04 T L 5.88 57.02 0.00	100 4988 25 10 301.03 0.00 6254.61 0.00 41.88 0.00	100 4976 25 8 T L 1.42 1772.46 0.00 25.23 0.00	100 4916 25 6 -10.47 T L 24.34 110.61 0.42	100 4882 20 10 -6.26 T L 9.38 63.97 0.00	100 4918 20 8 T L 2.27 T L 2.77 48.25 2.30	0.40 60.61 1.65	|N | |A| |C| |F | t[s] Gap t[s] Gap t[s]	Instances B+VI+H C PH

This contradicts the constraint (2.12), which states that xfk d k = 1.

Figure 4.9: CPU time needed to solve the LP-relaxation of the PF, C and DW formulation.

Benders reformulation for the node-capacitated VNFPRP high-quality solutions are provided in the initialization phase. To this end, we use the MILP-based heuristic introduced in Chapter 3, whose solution is handed over to the solver immediately before solving the rst LP of the relaxed master problem.

Acknowledgements

I would like to convey my sincere gratitude to my supervisors Ivana Ljubi¢, Nancy Perrot, and Éric Gourdin, for their advice, availability, and for providing me scientic and moral support throughout this work. They were able to transmit to me the knowledge, the rigor, the motivation, and the passion for research in general and combinatorial optimization in particular. I am grateful for all the eorts that have been deployed in order to make this project successful and for all the constructive remarks and criticism that helped me improving my work and acquiring a lot of knowledge. I was very honored that Mrs. Brigitte Jaumard, Professor at the Concordia University, Canada, and Mr. Laurent Alfandari, Professor at the ESSEC Business School of Paris, accepted to review my thesis. My deep gratitude goes to Mr. Ali Ridha Mahjoub, Professor at université Paris-Dauphine, to Mr. Fabio Furini, Researcher at the Italian National Research Council (CNR), and to Mr. Markus Leitner, Associate Professor at Vrije Universiteit Amsterdam, for having accepted to examine my thesis.

I address, likewise, my warmly thanks to all my colleagues in Orange Labs and LAM-SADE for the welcoming and friendly environment, for the lively discussions, for their encouragement, their friendship, and the scientic exchanges that greatly contributed to my scientic and personal development during these years of my thesis.

I would like to especially thank Abdelatif, Quentin, Adrien, Wesley, Fatma, Raja,

Extended formulations

Sections 2.2 and 2.3). Algorithm 6: Labeling algorithm based on dynamic programming -Inspired from Dijkstra algorithm proposed for the SPP.; -The algorithm starts from the source node.;

for each iteration, i.e., going from node i to node j using arc (i, j) do -A label l j is created for node j based on l i , in order to store information on the resources and the predecessor arc (i, j).;

-Only feasible labels are created, i.e., satisfying all resource constraints.;

for each node do -Only the non-dominated labels (i.e., pareto-optimal) are kept according to the following dominance rule: ;

if Both are at the same node then if for each feasible extension of l 2 , there is also do -A feasible extension of l 1 where res l 1 r ≤ res l 2 r with r ∈ R, for exactly the same resources.

end then l 1 dominates l 2 .; end end end end -The algorithm stops when there are no more unprocessed labels for all nodes.; -It then checks whether the destination node could be reached.; -If yes, return all shortest paths from s k to d k satisfying the resource constraints.

Remark: We must have at least one non-decreasing resources (without negative cycles) to use this algorithm.

Reduced cost method For each commodity k ∈ C, we run Yen's algorithm [START_REF]Finding the k shortest loopless paths in a network[END_REF] in order to generate all latency-constrained elementary paths between s k and d k . We keep all these columns in a column pool, denoted by P k . First, the restricted master problem is initialized by a subset of columns from P k , which are then added to P k and removed from P k . At each iteration of the CG procedure, the restricted master problem is solved, and a dual solution is generated. Based on this solution, the reduced cost RC(p) of each path p ∈ P k is calculated. If the path p ∈ P k has a negative reduced cost, then p can be added to the restricted master (i.e., to the set P k ) and deleted from Extended formulations this case the path p is deleted from P k . Sign -" in columns Col", Iter", Gap" and Nodes" in tables indicates that no information can be provided and that no solution is found. Moreover, T L in t[s]" columns illustrates the fact that the time limit is reached. The relative gap improvement is reported as Gap = (GLB -LB)/LB * 100%, where GLB denotes the global lower bound found by each setting and LB the worst known lower bound (the minimum overall global lower bounds). The Uncapacitated VNFPR can then be modeled as follows:

P dh

(P ) : min

(5.10)

Unsplittable routing paths

Integrality constraints t k uv ∈ {0, 1} guarantee that the s k -d k ow cannot be split, i.e., that there is a single path used for routing the ow. In the following, we show that this constraint can be relaxed. Let (P') denote the model (P) in which integrality constraints (5.10) are replaced by

The major result of this section, which will also allow us to apply the Benders decomposition scheme is summarized as follows:

Theorem 5.1 If the compact formulation (P ) has an optimal fractional solution then it must necessarily admit an integer solution with the same objective value.

The proof of this theorem is provided at the end of the section. We start by proving Lemma 5.2, in which we focus on a single commodity, and assume that the capacity upper bounds on arcs corresponding to constraints (2.7) and capacity lower bounds on nodes corresponding to constraints (2.10) are pre-specied.

Lemma 5.2 Let ĉuv ∈ N ∪ {0} be an arbitrary capacity function on arcs (u, v) ∈ A, ĉu ∈ {0, 1} be an arbitrary lower bound capacity on nodes u ∈ N and l ∈ N a owlatency limit for a ow to be sent from a given source s ∈ N to a destination d ∈ N , s = d. The following system of linear inequalities is either infeasible, or it admits a feasible binary solution.

166

Benders reformulation for the node-capacitated VNFPRP

Proof. Assume that the system of linear inequalities is feasible, and let t be its fractional feasible solution. By the ow decomposition theorem, we can decompose the solution t into t1 , t2 , . . . , tQ associated with fractional paths P = {p 1 , p 2 , . . . , p Q }, each of them carrying r q units of ow, such that Q q=1 r q = 1. Let tq,uv be the arc variables that compose tq and the path p q and let tq,uv ∈ {0, 1} be obtained by applying the ceiling function to values of tq,uv , i.e., tq,uv = tq,uv , for all (u, v) ∈ p q , 1 ≤ q ≤ Q.

For each node u ∈ N such that ĉu = 1, by the feasibility of t, all paths p q ∈ P must pass through node u.

The ow preservation constraints (5.11) are clearly satised by each tk q , 1 ≤ q ≤ Q. It is easy to see that the ow preservation constrains remain satised for each individual path dened by tk q , 1 ≤ q ≤ Q. Hence, to satisfy constraints (5.11), (5.13) and (5.14), one could take any of the paths p 1 , . . . , p Q as a feasible binary solution.

What remains to show is that at least one of the paths {p 1 , . . . , p Q } must also satisfy the latency constraint (5.12). We observe that tq,uv has the same value for each arc (u, v) composing the path p q , for each 1 ≤ q ≤ Q, more precisely, this value is equal to r q dened above. As the sum of the ow is integer, at source and destination nodes we must have at least two paths with fractional ows. Suppose that none of the paths from P satisfy the latency constraints, i.e.: l(p) > l ∀p ∈ P , where l(p) = uv∈p l uv denotes the latency of the path p ∈ P . Then:

Node-Precedence inequalities Theorem 5.4 The constraints (5.17) are valid for (P ):

(5.17)

Proof. For given positions i and j in the ordering of F k , with i = 1, 2, . . . , |F k | -2 and j = i+2, . . . , |F k |, if both functions f i and f j are installed at node u for commodity k, then to satisfy the precedence constraints, each function f l between f i and f j , with l = i+1, i+2, . . . , j -1, must be installed at node u as well, otherwise we allow creation of cycles, and hence the violation of precedence constraints (see Figure 5.1).

.1: Explanation of node-precedence inequalities.

We notice that even though these cuts are polynomial in number (there are O(|C||F

of these inequalities), they may impose a signicant burden to the MILP solvers because of their size. We therefore also consider the aggregated version of these constraints.

Aggregated node-precedence inequalities Theorem 5.5 The linear inequalities (5.18) are valid for (P ):

Proof. The constraints are obtained by summing up (5.18) over all l ∈ {i + 1, . . . , j -1}, and hence they are valid for (P).

Problem reformulation using Benders cuts 169

Node-capacity inequalities Theorem 5.6 The linear inequalities (5.19) are valid for (P ):

(5. [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF] Proof. For a given position i in the ordering of

the function f i is installed at node u for the commodity k, then the number of functions that must be installed before and after f i must be less than or equal to the capacity c u of node u, otherwise the node-capacity constraints are violated.

Problem reformulation using Benders cuts

The size of the compact MILP model presented in Section 5.1 easily becomes intractable, due to the fact that the model contains O(|N ||F ||C| + |A||C|) variables and constraints. Hence, to deal with instances of realistic size, it is necessary to rely on a problem decomposition in which smaller and easier-to-solve components are solved separately, and then the gained information is combined, and the process is repeated in an iterative way until convergence.

In this chapter we focus on Benders decomposition approach, which can be eciently exploited thanks to our result given in Section 5.1. The Benders master problem aims at nding an optimal placement of virtual network functions associated to each commodity and deciding in which order these functions should be visited, while guessing that a feasible routing path can be found for such a placement. Hence, in the master problem we decide on the value of (x, y, w) variables, We rst observe that this linear system of inequalities is separable, so that the feasibility check can be performed independently for each commodity k ∈ C. For a given binary vector (x, ŷ), the feasibility subproblem for a commodity k ∈ C can be formulated as the following linear program (we omit the superscript k for the sake of better readability): min 0

Hence, the solution (x, ŷ, ŵ) is feasible, i for each k ∈ C, there are sucient arc and node capacities to route the s k -d k ow and the latency of the routing path is at most l k . According to LP duality, if the dual of the k-th Benders subproblem is unbounded, its primal is infeasible, and correspondingly, the point (x, ŷ, ŵ) has to be cut o.

After associating dual variables α, β, γ and δ to constraints (5.20)-(5.23), respectively, the LP dual of the k-th Benders subproblem reads as follows:

Path-based formulation

The path formulation has been introduced in Chapters 3 and 4. In the Uncapacitated version of the problem the path formulation is a bit dierent from the one already introduced in previous chapters. The objective function is not the same, and the constraints (5.7) should be added to the model. Also, the VNF-capacity, the nodecapacity and the conict constraints are relaxed. Then, the Uncapacitated VNFPR problem can be modeled as the following MILP:

Heuristic

We generate up to κ elementary shortest s k -d k -paths for each commodity k ∈ C, using Yen's algorithm [START_REF]Finding the k shortest loopless paths in a network[END_REF]. The paths whose length is below l k are added to the model in the order from the shortest one to the longest one.

Obviously, if Yen's algorithm would generate all possible s k -d k -paths with latency less or equal to l k and there would be less than κ of them, our resulting compact MILP model would be an exact reformulation of the problem. This however does not happen very often in practice (we usually keep κ ≤ 50). Nevertheless, the quality of solutions obtained by this simple MILP heuristic is quite satisfying, as reported in the following section.

Computational results

The purpose of this section is to test the scalability and the eciency of the proposed Branch-and-Benders-Cut computational framework and to compare it to other exact methods, using an o-the-shelf solver. To this end, we use the commercial solver CPLEX and compare our approach to two available options: rst, solving the problem as a compact formulation (P) introduced in Section 5.1, and second, applying Automatic Benders decomposition [START_REF] Bonami | Implementing Automatic Benders Decomposition in a Modern MIP Solver[END_REF] to the formulation (P').

Our tests are conducted so as to evaluate: the performance of Benders decomposition in terms of computing time; the size limit of realistic instances to be solved to optimality in a reasonable time limit; the quality of the MILP-based heuristic, and, nally, the possible benet of valid inequalities proposed in Section 5.1.4 in improving the computing times or reducing the nal gaps. We also investigate the performance of our models for problem instances with node-capacity and conict constraints.

All the tests reported in this section were made using a machine with Intel(R) Xeon(R) CPU E5-2650 v2 processor clocked at 2.60GHz and 252GB RAM, under Linux operating system. All methods are implemented using the Python API for CPLEX, which is run in single-threaded mode and all CPLEX parameters were set to their default values. A default time limit of three hours is set for each instance from the benchmark set described in Section 5.4.1.1, and of one hour for each instance from the benchmark set described in Section 5.2. The default memory limit was set to 20GB.

The following settings were tested in our computational experiments:

• C: the compact MILP formulation (P) proposed in Section 5.1;

• RC: the relaxed compact MILP formulation (P') in which the ow variables t are relaxed according to Corollary 5.3;

• AB: automatic Benders decomposition available in Cplex [START_REF] Bonami | Implementing Automatic Benders Decomposition in a Modern MIP Solver[END_REF], applied to the model Benders reformulation for the node-capacitated VNFPRP

• B+PH: setting B in which initial solutions are obtained using the MILP-Heuristic described in Section 5.3. The number of generated paths by Yen's algorithm is limited to κ = 50 paths per commodity and the time limit for running the MILP-heuristic is set to 900 seconds;

• B+VI+PH: setting B+PH in which Valid Inequalities from Section 5.1.4 are additionally used to initialize the model.

Benchmark instances

In order to test the performance of the proposed methods, we generate two set of instances. The rst set is generated using Erd®s-Rényi graphs, and the second set is obtained from the well-known library of telecommunication network instances, called SNDlib [125].

Instances derived from Erd®s-Rényi graphs

These instances are generated as Erd®s-Rényi graphs using the Python library Net-workX. We vary the number of nodes |N | ∈ {25, 50, 100} and consider the graph density d ∈ {0.5, 0.9}. For each arc (i, j) ∈ A, the latency l ij ∈ N is chosen uniformly at random from the interval [1, 25], and we set l ij = l ji . For each node u ∈ N , its node activation cost ψ u ∈ N is chosen uniformly at random from [START_REF] Fischetti | Benders decomposition without separability: A computational study for capacitated facility location problems[END_REF]1000]. The total number of VNFs |F | is set to 10, and for each u ∈ N, f ∈ F , the installation cost ψ f u ∈ N is chosen uniformly at random from [10,[START_REF] Mckeown | Openow: enabling innovation in campus networks[END_REF] . We vary the total number of commodities, by considering |C| ∈ {10, 15, . . . , 50}. For each commodity, its distinct source and destination nodes are chosen randomly from N . To dene the set F k of VNFs associated to each commodity, we randomly choose between {6, 8, 10} functions from the set F and order them randomly. The latency l k is set to 1.5 times the length of the shortest path (with respect to the values l ij ). .

Instances derived from the SNDlib

204

Benders reformulation for the node-capacitated VNFPRP

Conclusions

In this chapter we studied two variants of the Virtual Network Functions Placement and Routing problem (VNFPRP). We provided theoretical results that allowed us to reformulate the problem using Benders decomposition. We proposed three families of valid inequalities to strengthen the LP-bounds.

Conclusions

In this dissertation, we have studied the Virtual Network Functions Placement and Routing Problem (VNFPRP), for which the sum of the function installation and node activation costs has to be minimized. The studied problem was considered with routing and latency constraints, node and function capacity constraints, incompatibility and chaining constraints. We have shown that the problem is strongly NP-hard even for its simplest version.

In the rst part of the thesis, we have investigated the basic properties of the problem and proposed a compact MILP formulation. This formulation does not seem to be strong enough for nding a solution using an o-the-shelf solver. To tackle the problem from a computational perspective, we have proposed a path-based heuristic that provides optimal solutions for some realistic instances from the literature. Afterwards, we have proposed two extended formulations to model the problem: path formulation and Dantzig-Wolfe formulation. In order to strengthen the LP-bounds, we have proposed several families of valid inequalities and have demonstrated their benets. We have shown that the LP-bounds of Dantzig-Wolfe formulation are stronger than the LP-bounds of the path formulation. We have presented a branching scheme for each formulation and developed a respective Branch-and-Price algorithm. We have computationally compared both algorithms with the MILP compact formulation and the automatic Benders of Cplex.

In the last part of the dissertation, we have studied a variant of the problem in which the VNFs-capacity and conict constraints are relaxed. We provided theoretical results that allowed us to reformulate the problem using Benders decomposition and three families of valid inequalities to strengthen the LP-bounds. All these ingredients have been combined in a Branch-and-Benders-Cut framework and tested on a set of realistic benchmark instances.

As perspectives, there are diverse directions for which our future research associated with the VNFPRP can be conducted.

MOTS CLÉS

Optimisation combinatoire, Génération de colonnes, Branch-and-Price, inégalités valides, décomposition de Benders, Branch-and-Benders-Cut, Heuristique, Fonctions virtualisées de réseau, Réseaux définis par logiciel, Chaînage de fonctions de service.

RÉSUMÉ

La virtualisation des fonctions réseau et les réseaux définis par logiciel sont deux nouvelles technologies prometteuses qui émergent dans la nouvelle génération des réseaux de télécommunication. Leur introduction permet la minimisation du temps de traitement des services, et un gain en énergie et en coûts. Dans cette thèse, nous étudions un problème confronté par les fournisseurs des services réseau, intitulé le problème de placement des fonctions virtuelles et de routage des services dans les réseaux définis par logiciel.

Nous commençons par prouver que le problème est NP-Difficile au sens fort, même en considérant une seule commodité, et en relâchant les contraintes de capacité sur les noeuds et sur les fonctions virtuelles, et aussi les contraintes de latence et de précédence. Nous proposons ensuite une formulation PLNE (Programmation linéaire en nombres entiers) compacte pour modéliser le problème. Cette formulation ne semble pas être assez forte pour trouver des solutions réalisables en un temps raisonnable à l'aide d'un solveur standard. Afin de remédier à cela, nous fournissons, une heuristique basée sur une formulation PLNE.

Nous proposons également, deux formulations PLNE étendues pour modéliser le problème et nous définissons l'algorithme de Branch-and-Pricer associé. En plus, nous définissons plusieurs familles d'inégalités valides afin de renforcer la relaxation linéaire. Nous comparons les deux algorithmes proposés et nous discutons leur performance.

Ensuite, nous étudions une variante du problème et nous présentons des résultats théoriques qui nous permettent de la reformuler en appliquant la décomposition de Benders. Nous renforçons cette reformulation par un ensemble d'inégalités valides. Tout cela est combiné dans algorithme de Branch-and-Benders-cut que nous avons testé et comparé avec l'algorithme de Benders automatique inclus dans le solveur commercial Cplex.

Les résultats numériques indiquent que nos approches de décomposition et de reformulation sont plus efficaces par rapport aux deux méthodes (la formulation compacte et l'algorithme de Benders automatique) fournies par un solveur standard sur un ensemble d'instances réalistes, à la fois en terme de temps CPU et en terme de qualité de la solution.

Les résultats indiquent également que notre heuristique fournit des solutions de haute qualité. We start proving that the VNFPRP is NP-hard in a strong sense, even for a single commodity and without nodecapacity, VNFs-capacity, latency, and precedence constraints. We provide a compact Mixed Integer Linear Programming (MILP) formulation to model it. This formulation does not seem to be strong enough for finding a solution using an offthe-shelf solver. In order to tackle the problem from a computational perspective, we provide MILP-based heuristic. The obtained results indicate that our MILP-heuristic provides high-quality solutions.

Moreover, we propose two extended formulations for the problem and derive a Branch-and-Price algorithm. We also provide several families of valid inequalities to strengthen the LP-bounds. We present computational results and discuss the performance of each algorithm.

Furthermore, we discuss a variant of the problem, and we provide theoretical results that allow us to derive Benders reformulation of the problem, along with several families of valid inequalities. These ingredients are combined in a Branchand-Benders-Cut framework and computationally tested on a wide range of realistic instances and compared with the Automatic Benders decomposition provided by Cplex.

Computational results indicate that our decomposition and reformulation approaches are more efficient compared to the two methods (the MILP compact formulation and the automatic Benders) provided by the off-the-shelf solver on a set of realistic instances, both in terms of CPU time and overall solution quality.

KEYWORDS

Combinatorial optimization, Column generation, Branch-and-Price, Valid inequalities, Benders decomposition, Branch-and-Benders-cut, Heuristic, Virtual Network Functions, Software-Defined Networking, Service Functions Chaining.