Pedro Liguori

Ridha Mahjoub

Ruslan Sadykov

Eduardo Uchoa

A Branch-Cut-and-Price Algorithm for the Location-Routing Problem

Literature on LRP

A combination of two central OR problems ≈3000 papers in Google Scholar with both "location" and "routing" in the title Important recent works [START_REF] Ii Belenguer | A branch-and-cut method for the capacitated location-routing problem[END_REF] -important valid inequalities & Branch-and-Cut; [START_REF] Baldacci | An exact method for the capacitated location-routing problem[END_REF] -exact "enumeration" & column generation approach [START_REF] Contardo | An exact algorithm based on cut-and-column generation for the capacitated location-routing problem[END_REF] -state-of-the-art exact algorithm [START_REF] Schneider | Large composite neighborhoods for the capacitated location-routing problem[END_REF] -state-of-the-art heuristic [START_REF] Schneider | A survey of the standard location-routing problem[END_REF] -the latest survey on LRP

Our study

Recently, large improvement in exact solution of classic VRP variants [START_REF] Pecin | Improved branch-cut-and-price for capacitated vehicle routing[END_REF]] [Pecin et al., 2017a[START_REF] Vii S | A bucket graph based labeling algorithm with application to vehicle routing[END_REF]] [Pessoa et al., 2018a] A generic Branch-Cut-and-Price VRP solver [START_REF] Vi Pessoa | A generic exact solver for vehicle routing and related problems[END_REF] incorporates all recent advances vrpsolver.math.u-bordeaux.fr

This solver can be applied to the LRP However, problem-specific cuts are necessary for obtaining the state-of-the-art performance

We review existing families of cuts and propose new ones Formulation

λ i r , i ∈ I, r ∈ R i , equals 1 iff route r is used for depot i a r e , e ∈ E ∪ F , r ∈ ∪ i∈I R i , equals 1 iff edge e is used by r y i , i ∈ I, equals 1 iff route depot i is open z ij , i ∈ I, j ∈ J, equals 1 iff client j is assigned to depot i min i∈I f i y i + i∈I r ∈R i e∈E∪F c e a r e λ i r i∈I z ij = 1, ∀ j ∈ J, r ∈R i e∈δ(j) a r e λ i r = 2z ij , ∀ i ∈ I, j ∈ J j∈J d j z ij ≤ w i y i , ∀ i ∈ I, z ij ≤ y i , ∀ i ∈ I, j ∈ J, (z, y , λ) ∈ {0, 1} K Given a subset of clients C ⊂ J, i∈I r ∈R i e∈δ(C) a r e λ i r ≥ 2 • i∈C d i Q .
Separation (embedded in the VRP solver)

CVRPSEP library [Lysgaard et al., 2004] Chvátal-Gomory Rank-1 Cuts [Jepsen et al., 2008] [ [START_REF] Pecin | Limited memory rank-1 cuts for vehicle routing problems[END_REF] Each cut is obtained by a Chvátal-Gomory rounding of a set C ⊆ J of set packing constraints using a vector of multipliers ρ (0

< ρ j < 1, j ∈ C): i∈I r ∈R i     j∈C ρ j e∈δ(j) 1 2 a r e     λ i r ≤     j∈C ρ j    
All best possible vectors ρ of multipliers for |C| ≤ 5 are given in [START_REF] Pecin | Limited memory rank-1 cuts for vehicle routing problems[END_REF].

Non-robust in the terminology of [START_REF] Pessoa | Robust branch-cut-and-price algorithms for vehicle routing problems[END_REF] Separation (embedded in the VRP solver)

A local search for each vector of multipliers.

Depot Capacity Cuts [START_REF] Ii Belenguer | A branch-and-cut method for the capacitated location-routing problem[END_REF] If a subset of clients C ⊂ J cannot be served by a subset of depots S ⊂ I, Separation (in the VRP solver callback)

A heuristic algorithm: combination of GRASP and local search.

Covering inequalities for depot capacities Let W = i∈I w i and D = j∈J d j . We should have

i∈I w i y i ≥ d(J) ⇒ i∈I w i (1 -y i ) ≤ W -D
We can generate any valid inequality for this knapsack. For example covering inequalities: given a subset of depots S ⊂ I,

i∈S w i > W -D, i∈S (1 -y i ) ≤ |S| -1
Separation (in the VRP solver callback)

We optimize an LP which looks for the most violated inequality which is satisfied by all integer solutions of the knapsack.

Route Load Knapsack Cuts (RLKC)

x i q -number of routes with load of exactly q ≤ Q units leaving depot i ∈ I. Then:

Q q=1 qx i q ≤ w i . (1)
Any valid inequality for (1) is valid for the LRP.

Non-robust in the terminology of [START_REF] Pessoa | Robust branch-cut-and-price algorithms for vehicle routing problems[END_REF] First separation algorithm Chvátal-Gomory rounding of (1).

1/k -facets of the master knapsack polytope

Theorem ( [START_REF] Aráoz | Polyhedral neopolarities[END_REF])

The coefficient vectors ξ of the knapsack (non-trivial) facets ξx ≤ 1 of n q=1 qx q = n with ξ 1 = 0, ξ Q = 1 are the extreme points of the following system of linear constraints

ξ 1 = 0, ξ Q = 1, ξ q + ξ Q-q = 1 ∀1 ≤ i ≤ n/2, ξ q + ξ t ≤ ξ q+t whenever q + t < n. Definition A knapsack facet ξx ≤ 1 is called a 1/k -facet if k is the smallest possible integer such that ξ q ∈ {0/k , 1/k , 2/k , . . . , k /k } ∪ {1/2}.
Second separation algorithm 1/6-and 1/8-facets can be efficiently separated using the algorithm by [START_REF] Chopra | A concise characterization of strong knapsack facets[END_REF] Taking into account of RLKCs in the pricing Let μ(q) be the contribution of RLKCs to the reduced cost of a route variable with load q

Pricing problem: Resource Constrained Shortest Path

It is solved by a labelling algorithm, each label L is

(c L + μ(q L ), j L , q L ) Dominance relation L L if cL ≤ cL , j L = j L , q L ≤ q L (2)
is valid, as μ(q) is non-decreasing

Completion bounds can still be efficiently used as μ(q) is super-additive

Other components of the Branch-Cut-and-Price

Bucket graph-based labelling algorithm for the RCSP pricing [START_REF] Righini | Symmetry helps: Bounded bi-directional dynamic programming for the elementary shortest path problem with resource constraints[END_REF]Salani, 2006] [S. et al., 2017] Partially elementary path (ng-path) relaxation [Baldacci et al., 2011a] Automatic dual price smoothing stabilization [Wentges, 1997[START_REF] Pessoa | Automation and combination of linear-programming based stabilization techniques in column generation[END_REF] Reduced cost fixing of (bucket) arcs in the pricing problem [START_REF] Ibaraki | A dynamic programming method for single machine scheduling[END_REF]] [Irnich et al., 2010[START_REF] Vii S | A bucket graph based labeling algorithm with application to vehicle routing[END_REF] Enumeration of elementary routes [START_REF] Baldacci | An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts[END_REF] Multi-phase strong branching [START_REF] Pecin | Improved branch-cut-and-price for capacitated vehicle routing[END_REF] On depot openings (largest priority) On number of vehicles for each depot On number of clients per depot On assignment of clients to depots On edges of the graph

Computational results

Open instances solved to optimality. Could not be solved by the state-of-the-art [START_REF] Contardo | An exact algorithm based on cut-and-column generation for the capacitated location-routing problem[END_REF] in 5-97 hours

  Figure: LRP instance: G = (I ∪ J, E ∪ F )

  one vehicle from a depot i ∈ I \ S should visit C: i∈I\S r ∈R i e∈δ(C)

  Route Load Knapsack Cuts reduce the gap but not yet worth to include in the VRP solver An extension to the Two-Echelon Capacitated Vehicle Routing problem allows us to double the size of instances which can be solved to optimality[START_REF] Iv Marques | An improved branch-cut-and-price algorithm for the two-echelon capacitated vehicle routing problem[END_REF] Perspectives Improve separation of Route Load Knapsack Cuts A polyhedral study is needed for the Multi Capacitated Depot Vehicle Routing Problem.You can use the VRP solver to test new families of cuts for vehicle routing problems within state-of-the-art Branch-Cut-and-Price!

	Conclusions
	A large improvement over the state-of-the-art for the LRP
	Set by applying the VRP solver and providing callbacks for Instance Optimum Time [Prins et al., 2006] 100x5-1b 213568 10m05s problem-specific cuts
	100x10-1a	287661	1h32m
	100x10-1b	230989	1h38m
	100x10-3a	250882	1h17m
	100x10-3b	203114 11h01m
	200x10-1a	474702 20m42s
	200x10-1b	375177	1h55m
	200x10-2a	448005	4h45m
	200x10-2b	373696	5h53m
	[Tuzun and Burke, 1999] P113112	1238.24	2h29m
	P131112 2E-CVRP demo is available on	1892.17 36m52s
	P131212 vrpsolver.math.u-bordeaux.fr 1960.02 34m59s
	Underlined: improved solutions over [Schneider and Löffler, 2019]

instances by[START_REF] Prins | Solving the capacitated location-routing problem by a grasp complemented by a learning process and a path relinking[END_REF] with 5-10 depots and 50-200 clients. Time limit 3 hours


Pascal Bruel

. Mach zero vs low Mach flows

                                         u u u
                    u
Compressible Euler equations for a polytropic ideal gas: dimensionless form [START_REF] Müller | 30th Computational Fluid Dynamics Lecture Series, von Karman Institute for Fluid Dynamics[END_REF]:

Pascal Bruel -NLPDFP Workshop-23 September 2019 -Borovoe (Burabay) -Kazakhstan

From now on, we shall drop the superscript * for denoting the dimensionless variables and it is assumed that the low Mach number asymptotic analysis can be considered as a regular perturbation.

(

) (0) (1) (1) 0
(0) (0) (1) (1) (2) (2)
(0) (0) (1) (1) (2) (2) ˆˆˆ( , ; ) ( ) ( , ) ( ) ( , ) ( ) ( , ) ...... ˆˆˆ( , ; ) ( ) ( , ) ( ) ( , ) ( ) ( , ) ..... ˆˆˆ( , ; ) ( ) ( , ) ( ) ( , ) ( ) ( , p t M M p x t M p x t M p x t u t M M u x t M u x t M u x t t M M x t M x t M x                         x x x ( ) (2) (2)
) ...... the scaling functions are chosen such that:

( ) inject these expansions into the governing equations ând order the powers of (i) Ŝo each independent variables is expanded in terms of a series ( ) ŵhere is the small parameter, for instance:

M M 
The reference Mach number is chosen as the small parameter

                     convective     u u u u u
One-scale expansion at low Mach: the hydrodynamic limit

(0) (0) (0) (0) (0) (0) (0) (2) (0) (0) (0) (0) (0) (0) (0) (0) (0)
The system for O(1) variables reads as:

.( ) 0

.( )

. . 1 ( , ) ( , ) ( )
This represents the Mach zero system

t p t T dp T t dt x t T x t p t                                   u u u u u .
With periodic boundary conditions, the kinetic energy defined by: remains constant in time 2Choose now ( / ( / )) So = is now an time scale, the dimensionless system reads now as:

1 .( ) 0 1 1 .( ) 1 .( ) 0 ˆref ref ref ref ac t L c t t t M p t M M E H t M                       acoustic     u u u u u Pascal Bruel -NLPDFP Workshop -23 September 2019 -Borovoe (Burabay) -Kazakhstan
One-scale expansion at low Mach: the acoustic limit

(1) (0) (0) 1 (0) 1 (0) 0 (1) (0) (0)
Then, the continuity equation reads as (keeping only the leading terms):

.( ) 0 to be satisfied for all when so:

Term 0 Term .( ) 0 co M t t M t M t                            u M M 0 u (0) (0) (1) (1) (0) 2 (0) 2 (0) (0) 0 0
mbined with the energy equation it yields after some algebra the : Low Mach flows: two-time-scale expansion

0 with / p t p c = 0 c p t            
(1) (0) (0) (0) (0) (0) (0)
The first order energy equation can be expressed as:

.( ) .( ) ( 1)( ) Pr Re Differentiating the above equation with respect to and substracting times the dive

p dp p T q dt p                    u 2 2 2 (1) (0) (0) (0) (1) (0) (0) 2 (0)
rgence of the first order momentum equation yields:

( ) ( ) .( ) ( 1)
with ( , ) This is a wave equation and its source is the change of the

p q p t c p c t              over acoustic time    x
See the seminal paper by Klein (1995) for a detailed analysis of the two-space-scale expansion. Examples of results, test cases definitions with a continuous quenched mean reaction rate (Corvellec et al., 1999)

u c Sc A    3 * * 3 2 2 ( / ( / ) ) 1.1886 1.1886 ( - )(1 ) ( - )(1 ) / (1 ) (1 )
The corresponding theoretical turbulent flame velocity is 0.5 m/s (case 1) and 10 m/s (case 2) The density may vary significantly in relation with the exact nature of the flow at hand (contact dicontinuities, reacting flow).

Together with proper initial and boundary conditions

The continuous system of PDE's at hand: the Euler equations with a co-located formulation 2 2

.( ) 0

.( )

.( ) 0 Proposition for a pressure-based algorithm for simulating non reacting flows (Moguen et al., 2019) Predictor step (pressure frozen): solving the continuity and the momentum equations.

1 1 ; 2 2 ( 1) v t p t E H t p E e c T H = E p e                              
Corrector step (density frozen): solving the energy equation to get the pressure correction with input from the momentum equation.

Final conclusion

The versatility of pressure based approaches has been demonstrated: able to cope with zero Mach as well as with low Mach flow systems. The MIAU based scheme works also well for high Mach flows (not shown here).

Future work will concentrate on:

Artificial compressibility: increasing the convergence rate by working on the pseudo waves amplitude system at the boundaries.

MIAU like pressure based approach: extending it to reacting flows.

  Workshop -23 September 2019 -Borovoe (Burabay) -KazakhstanOne-scale expansion at low Mach: the acoustic limit With periodic boundary conditions, the acoustic energy defined by: remains constant in time. Here, the tilde denotes an averaged value over a region of characteristic length of 1/ 𝑀

  Low Mach flows: two-time-scale expansion Basic idea: introducing a finite artificial sound speed in the system Continuity equation equation for the pressure• Proposed byChorin (1967) for inert flows • Unsteady inert flows(Soh et al. 1988, McHugh et al. 2005) • Extended to steady zero Mach number reacting flows byBruel et al. (1996).• Extended to unsteady zero Mach number reacting flows byCorvellec et al. (1999).• Combined with an hybrid-mesh finite volume formulation byDourado et al. (2004) • Combined with discontinuous Galerkin finite element method by[START_REF] Bassi | [END_REF] It is well suited to transform a compressible density based code into a code able to deal with zero Mach inert or reacting flows. speed is given in constant density flow regions by a and the corresponding pseudo-Mach number is M 1 and its value can be controlled through the setting of .

  of low Mach flows: the MIAU pressure based approach Pascal Bruel -NLPDFP Workshop -23 September 2019 -Borovoe (Burabay) -Kazakhstan Important features of low Mach flows Any low Mach flow solution results a priori from the superimposition/interaction of a slow component (with a quasi-like Mach zero behavior) and a fast component (acoustic waves).

Pascal Bruel -NLPDFP Workshop-23 September 2019 -Borovoe (Burabay) -Kazakhstan

	Choose	t	ref	(	ref L	/	u	ref	)
	So	t	ref	=	t	conv	is the	time scale, the dimensionless
	One-scale expansion at low Mach: the hydrodynamic limit system reads now as:
	t						.(		) 0
		t				.(			)	2 .( 1 p M
		t	E						H	) 0

Pascal Bruel -NLPDFP Workshop-23 September 2019 -Borovoe (Burabay) -Kazakhstan

Pascal Bruel -NLPDFP Workshop-23 September 2019 -Borovoe (Burabay) -Kazakhstan

Pascal Bruel -NLPDFP Workshop-23 September 2019 -Borovoe (Burabay) -KazakhstanOne-scale expansion at low Mach: the hydrodynamic limit

Pascal Bruel -NLPDFP Workshop-23 September 2019 -Borovoe (Burabay) -KazakhstanOne-scale expansion at low Mach: the hydrodynamic limit One-scale expansion at low Mach: the acoustic limit

Pascal Bruel -NLPDFP Workshop -23 September 2019 -Borovoe (Burabay) -KazakhstanOne-scale expansion at low Mach: the acoustic limit

With a quenched mean reaction rate, the turbulent flame speed is unique and the flame structure can be determined semi-analytically(Sabel'nikov et al., 1998) Pascal Bruel -NLPDFP Workshop -23 September 2019 -Borovoe (Burabay) -Kazakhstan

Pascal Bruel -NLPDFP Workshop -23 September 2019 -Borovoe (Burabay) -Kazakhstan

Pascal Bruel -NLPDFP Workshop-23 September 2019 -Borovoe (Burabay) -Kazakhstan

Turbulent premixed reaction zone stabilized by a dump [START_REF] Nguyen | [END_REF] Low Mach flows: illustration of a reacting flows featuring turbulence and coherent motion related with thermoacoustic coupling

Low Mach flows simulations: what are the choices ?

Density-based approach: preconditioning techniques, modification of the diffusion matrix of the flux scheme.

Pressure-based approach: this approach originally developed to cope with Mach zero flows must be adapted: the energy equation plays a key role (Klein, 1995) and so will be used for establishing the pressure correction equation.

Updates with A few (~ 5) k-iterations are sufficient to converge the loop.

Low Mach flow configuration: nozzle flow